text:
using Chain: @chain
using Dates: Day, today
using Test

using Stonks.APIClients
using Stonks.Requests
using Stonks: UpdatableSymbol, RequestBuilderError, split_tickers_in_batches, build_fx_pair,  construct_updatable_symbols
using Stonks.Models: AssetPrice

include("test_utils.jl")

@testset "Request utilities" begin
  test_client = APIClients.YahooClient("abc")

  @testset "Build FX pairs" begin
    @test build_fx_pair("EUR/USD") == ("EUR", "USD")
    @test isa(build_fx_pair("EUR&USD"), ArgumentError)
    @test isa(build_fx_pair("EURO/USD"), ArgumentError)
    @test isa(build_fx_pair("EUR/USD/CAD"), ArgumentError)
  end

  @testset "Build UpdatableSymbol using an FX pair" begin
    @test first(construct_updatable_symbols("FOO/BAR")).fx_pair == ("FOO", "BAR")
    @test first(construct_updatable_symbols("FOOOO/BAR")).fx_pair |> ismissing
    @test last(construct_updatable_symbols(["EUR/USD", "USD/CAD"])).fx_pair == ("USD", "CAD")
  end

  @testset "Ticker batching when all dates are missing" begin
    tickers = @chain ["AAPL", "MSFT", "TSLA", "IBM"] map(t -> UpdatableSymbol(t), _)
    batches_max10 = split_tickers_in_batches(tickers, 10)
    batches_max2 = split_tickers_in_batches(tickers, 2)
    @test length(batches_max10) == 1
    @test length(batches_max2) == 2
  end

  # TODO - move it to test_utils.jl
  @testset "Complex ticker batching" begin
    tickers = complex_tickers()
    batches = split_tickers_in_batches(tickers, 2)
    groups = @chain batches map(batch -> map(t -> t.ticker, batch), _)
    @test length(batches) == 4
    @test groups == [["AAPL"], ["MSFT", "TSLA"], ["IBM", "GOOG"], ["NFLX"]]
  end

  @testset "Request building using APIResource" begin
    tickers = complex_tickers()
    resource = test_client.resources["price"]
    resource.max_batch_size = 2
    request_params = Requests.prepare_requests(tickers, resource; interval="1d")
    @test length(request_params) == 4
    @test sum(map(r -> length(r.tickers), request_params)) == length(tickers)
  end

  # @testset "Test request building using APIClient" begin
  #   tickers = complex_tickers()
  #   request_params = Requests.prepare_requests(tickers, test_client, interval = "1d")
  #   #@test isa(request_params, Dict)
  #   #@test keys(request_params) == keys(test_client.resources)
  # end

  @testset "Resolution of url params" begin
    url = "https://example.com/{resource}/{symbol}?lang=en&symbol={symbol}"
    new_url = Requests.resolve_url_params(url; resource="quote", symbol="AAPL", foo="bar")
    @test new_url == "https://example.com/quote/AAPL?lang=en&symbol=AAPL"
  end

  @testset "Resolution of query parameters" begin
    query_params = Dict(
      "interval" => "{interval}",
      "range" => "{range}",
      "symbols" => "{symbols}",
      "content" => "json",
    )
    params_resolved = Requests.resolve_query_params(
      query_params; interval="1d", range="30d", symbols="AAPL,MSFT"
    )
    @test params_resolved["interval"] == "1d"
    @test params_resolved["symbols"] == "AAPL,MSFT"
    @test params_resolved["range"] == "30d"
    @test params_resolved["content"] == "json"
  end

  @testset "Resolution of all request parameters (url, query, others)" begin
    resource = test_client.resources["price"]
    tickers = complex_tickers()
    request_params = Requests.prepare_requests(tickers, resource; interval="1d")
    @test length(request_params) == 4
    @test first(request_params).params["interval"] == "1d"
  end

  @testset "Resolution of all request parameters (url, query, others) for yahoo exchange rates" begin
    resource = test_client.resources["exchange"]
    resource.max_batch_size = 2
    tickers = construct_updatable_symbols(["EUR/USD", "USD/CAD", "FOO/BAR"])
    batches = Requests.split_tickers_in_batches(tickers, resource.max_batch_size)
    request_params = map(b -> Requests.resolve_request_parameters(b, resource; interval="1d"), batches)
    @test length(request_params) == 2
    @test first(request_params).params["symbols"] == "EURUSD=X,USDCAD=X"
  end

  @testset "Resolution of all request parameters (url, query, others) for alphavantage exchange rates" begin
    client = APIClients.AlphavantageJSONClient("abc")
    resource = client.resources["exchange"]
    tickers = construct_updatable_symbols(["EUR/USD", "USD/CAD", "FOO/BAR"])
    batches = Requests.split_tickers_in_batches(tickers, resource.max_batch_size)
    request_params = map(b -> Requests.resolve_request_parameters(b, resource; interval="1d"), batches)
    @test length(request_params) == 3
    @test first(request_params).params["from_symbol"] == "EUR"
    @test first(request_params).params["to_symbol"] == "USD"
  end

  @testset "Resolution of request parameters for alphavantage json client" begin
    client = APIClients.AlphavantageJSONClient("abc")
    resource = client.resources["price"]
    tickers = @chain ["AAPL", "MSFT"] map(t -> UpdatableSymbol(t), _)
    request_params = Requests.prepare_requests(tickers, resource)
    @test first(request_params).params["symbol"] == "AAPL"
  end

  @testset "Resolution of request parameters for exchange rate" begin
    query_params = Dict(
      "function" => "FX_DAILY", "from_symbol" => "{base}", "to_symbol" => "{target}"
    )
    params_resolved = Requests.resolve_query_params(query_params; base="EUR", target="USD")
    @test params_resolved["from_symbol"] == "EUR"
    @test params_resolved["to_symbol"] == "USD"
  end

  @testset "Validation of request parameters" begin
    resource = test_client.resources["price"]
    resource.query_params = Dict(
      "symbols" => "{symbols}", "from_symbol" => "{base}", "foo" => "{target}"
    )
    tickers = @chain ["AAPL", "MSFT"] map(t -> UpdatableSymbol(t), _)
    rp = Requests.resolve_request_parameters(tickers, resource; base="USD")
    @test isa(rp, RequestBuilderError)
  end

  @testset "Optimistic request resolution - success" begin
    params = Requests.RequestParams(;
      url = "http://example.com",
      tickers = map(UpdatableSymbol, ["AAPL", "MSFT"])
    )
    ok_values = [
      (
        params = params,
        value=[
          AssetPrice(; symbol="MSFT", date=Date("2022-02-16"), close=299.5),
          AssetPrice(; symbol="MSFT", date=Date("2022-02-17"), close=299.92)
        ],
        err=nothing
      ),
      (
        params = params,
        value=[AssetPrice(; symbol="MSFT", date=Date("2022-02-17"), close=290.73)],
        err=nothing,
      ),
      (
        params = params,
        value=nothing,
        err=ErrorException("Somwething went wrong")
      )
    ]
    responses = Channel(3)
    foreach(v -> push!(responses, v), ok_values)
    close(responses)
    res = Requests.optimistic_request_resolution(AssetPrice, responses)
    @test isa(res, Vector{AssetPrice})
  end

  @testset "Optimistic request resolution - failure" begin
    params = Requests.RequestParams(;
      url = "http://example.com",
      tickers = map(UpdatableSymbol, ["AAPL", "MSFT"])
    )
    values = [
      (params = params, value=nothing, err = ErrorException("Something went wrong")),
      (params = params, value=nothing, err = ErrorException("Another issue")),
    ]
    responses = Channel(2)
    foreach(v -> push!(responses, v), values)
    close(responses)
    res = Requests.optimistic_request_resolution(AssetPrice, responses)
    @test isa(res, ErrorException)
  end

  @testset "Conversion of date to range expression" begin
    today = Dates.today()
    conv_fn = Requests.convert_date_to_range_expr
    @test conv_fn(today) == "1d"
    @test conv_fn(today - Dates.Day(2)) == "5d"
    @test conv_fn(today - Dates.Day(7)) == "1mo"
    @test conv_fn(today - Dates.Day(28)) == "1mo"
    @test conv_fn(today - Dates.Day(30)) == "3mo"
    @test conv_fn(today - Dates.Day(89)) == "6mo"
    @test conv_fn(today - Dates.Day(179)) == "6mo"
    @test conv_fn(today - Dates.Day(180)) == "1y"
    @test conv_fn(today - Dates.Day(364)) == "1y"
    @test conv_fn(today - Dates.Day(365)) == "5y"
    @test conv_fn(today - Dates.Day(5000)) == "5y"
  end
end


text:
# As you leave the cave and reach open waters, you receive a transmission from
# the Elves back on the ship.
#
# The transmission was sent using the Buoyancy Interchange Transmission System
# (BITS), a method of packing numeric expressions into a binary sequence. Your
# submarine's computer has saved the transmission in hexadecimal (your puzzle
# input).
#
# The first step of decoding the message is to convert the hexadecimal
# representation into binary. Each character of hexadecimal corresponds to four
# bits of binary data:
#
# 0 = 0000
# 1 = 0001
# 2 = 0010
# 3 = 0011
# 4 = 0100
# 5 = 0101
# 6 = 0110
# 7 = 0111
# 8 = 1000
# 9 = 1001
# A = 1010
# B = 1011
# C = 1100
# D = 1101
# E = 1110
# F = 1111
#
# The BITS transmission contains a single packet at its outermost layer which
# itself contains many other packets. The hexadecimal representation of this
# packet might encode a few extra 0 bits at the end; these are not part of the
# transmission and should be ignored.
#
# Every packet begins with a standard header: the first three bits encode the
# packet version, and the next three bits encode the packet type ID. These two
# values are numbers; all numbers encoded in any packet are represented as
# binary with the most significant bit first. For example, a version encoded as
# the binary sequence 100 represents the number 4.
#
# Packets with type ID 4 represent a literal value. Literal value packets
# encode a single binary number. To do this, the binary number is padded with
# leading zeroes until its length is a multiple of four bits, and then it is
# broken into groups of four bits. Each group is prefixed by a 1 bit except the
# last group, which is prefixed by a 0 bit. These groups of five bits
# immediately follow the packet header. For example, the hexadecimal string
# D2FE28 becomes:
#
# 110100101111111000101000
# VVVTTTAAAAABBBBBCCCCC
#
# Below each bit is a label indicating its purpose:
#
# - The three bits labeled V (110) are the packet version, 6.
# - The three bits labeled T (100) are the packet type ID, 4, which means the
#   packet is a literal value.
# - The five bits labeled A (10111) start with a 1 (not the last group, keep
#   reading) and contain the first four bits of the number, 0111.
# - The five bits labeled B (11110) start with a 1 (not the last group, keep
#   reading) and contain four more bits of the number, 1110.
# - The five bits labeled C (00101) start with a 0 (last group, end of packet)
#   and contain the last four bits of the number, 0101.
# - The three unlabeled 0 bits at the end are extra due to the hexadecimal
#   representation and should be ignored.
#
# So, this packet represents a literal value with binary representation
# 011111100101, which is 2021 in decimal.
#
# Every other type of packet (any packet with a type ID other than 4) represent
# an operator that performs some calculation on one or more sub-packets
# contained within. Right now, the specific operations aren't important; focus
# on parsing the hierarchy of sub-packets.
#
# An operator packet contains one or more packets. To indicate which subsequent
# binary data represents its sub-packets, an operator packet can use one of two
# modes indicated by the bit immediately after the packet header; this is
# called the length type ID:
#
# - If the length type ID is 0, then the next 15 bits are a number that
#   represents the total length in bits of the sub-packets contained by this
#   packet.
# - If the length type ID is 1, then the next 11 bits are a number that
#   represents the number of sub-packets immediately contained by this packet.
#
# Finally, after the length type ID bit and the 15-bit or 11-bit field, the
# sub-packets appear.
#
# For example, here is an operator packet (hexadecimal string 38006F45291200)
# with length type ID 0 that contains two sub-packets:
#
# 00111000000000000110111101000101001010010001001000000000
# VVVTTTILLLLLLLLLLLLLLLAAAAAAAAAAABBBBBBBBBBBBBBBB
#
# - The three bits labeled V (001) are the packet version, 1.
# - The three bits labeled T (110) are the packet type ID, 6, which means the
#   packet is an operator.
# - The bit labeled I (0) is the length type ID, which indicates that the
#   length is a 15-bit number representing the number of bits in the sub-packets.
# - The 15 bits labeled L (000000000011011) contain the length of the
#   sub-packets in bits, 27.
# - The 11 bits labeled A contain the first sub-packet, a literal value
#   representing the number 10.
# - The 16 bits labeled B contain the second sub-packet, a literal value
#   representing the number 20.
#
# After reading 11 and 16 bits of sub-packet data, the total length indicated
# in L (27) is reached, and so parsing of this packet stops.
#
# As another example, here is an operator packet (hexadecimal string
# EE00D40C823060) with length type ID 1 that contains three sub-packets:
#
# 11101110000000001101010000001100100000100011000001100000
# VVVTTTILLLLLLLLLLLAAAAAAAAAAABBBBBBBBBBBCCCCCCCCCCC
#
# - The three bits labeled V (111) are the packet version, 7.
# - The three bits labeled T (011) are the packet type ID, 3, which means the
#   packet is an operator.
# - The bit labeled I (1) is the length type ID, which indicates that the
#   length is a 11-bit number representing the number of sub-packets.
# - The 11 bits labeled L (00000000011) contain the number of sub-packets, 3.
# - The 11 bits labeled A contain the first sub-packet, a literal value
#   representing the number 1.
# - The 11 bits labeled B contain the second sub-packet, a literal value
#   representing the number 2.
# - The 11 bits labeled C contain the third sub-packet, a literal value
#   representing the number 3.
#
# After reading 3 complete sub-packets, the number of sub-packets indicated in
# L (3) is reached, and so parsing of this packet stops.
#
# For now, parse the hierarchy of the packets throughout the transmission and
# add up all of the version numbers.
#
# Here are a few more examples of hexadecimal-encoded transmissions:
#
# 8A004A801A8002F478 represents an operator packet (version 4) which contains
# an operator packet (version 1) which contains an operator packet (version 5)
# which contains a literal value (version 6); this packet has a version sum of
# 16.
#
# 620080001611562C8802118E34 represents an operator packet (version 3) which
# contains two sub-packets; each sub-packet is an operator packet that contains
# two literal values. This packet has a version sum of 12.
#
# C0015000016115A2E0802F182340 has the same structure as the previous example,
# but the outermost packet uses a different length type ID. This packet has a
# version sum of 23.
#
# A0016C880162017C3686B18A3D4780 is an operator packet that contains an
# operator packet that contains an operator packet that contains five literal
# values; it has a version sum of 31.
#
# Decode the structure of your hexadecimal-encoded BITS transmission; what do
# you get if you add up the version numbers in all packets?

# read input
line = readline("input.txt")
bits = parse(BigInt, line; base=16)
pos = ndigits(bits; base=2)

# utility function to read n bits
function read_bits(n)
  global pos -= n
  pos < 0 && throw("pos is less than zero")
  mask = (1 << n) - 1
  (bits >> pos) & mask
end

# utility function - returns true if we're done reading bits
function iseob()
  # ensure 1 << pos will fit in an int (assuming 64-bit) - most likely, if
  # there are at least 63 bits remaining, we probably aren't done anyway
  pos > 63 && return false
  mask = (1 << pos) - 1
  iszero(bits & mask)
end

# do eet
total_versions = 0
while !iseob()
  global total_versions += read_bits(3)
  id = read_bits(3)
  if id == 4
    # literal value
    v = read_bits(5)
    while (v & 16) > 0
      v = read_bits(5)
    end
  else
    lentype = read_bits(1)
    if lentype == 0
      # Next 15 bits are the number of bits that make up sub packets, but we
      # don't really care about that at the moment... just consume the 15 bits
      read_bits(15)
    else
      # Next 11 bits are the number of sub packets. But we don't really care
      # about that at the moment, so, just consume the 11 bits
      read_bits(11)
    end
  end
end

println("Result: ", total_versions)

ext: jl


text:
#warlpiri-properties.jl - uses warlpiri numbers to test valid properties.

@testset "warlpiri-properties" begin
    @test isempty.(wmatrix) == map(x -> (x == (WP0001 → WP0000)) || (x == (WP0000 → WP1111)), wmatrix)

    const warlpiri_allreals =
        [WP0010 → WP0001, WP0011 → WP0010, WP0100 → WP0011, WP0101 → WP0100,
         WP0110 → WP0101, WP0111 → WP0110, WP1000 → WP0111, WP1001 → WP1000,
         WP1010 → WP1001, WP1011 → WP1010, WP1100 → WP1011, WP1101 → WP1100,
         WP1110 → WP1101, WP1111 → WP1110]

    @test SigmoidNumbers.isallreals.(wmatrix) == map((x) -> x in warlpiri_allreals, wmatrix)

    #roundsinf test          ->  0000  0001  0010  0011  0100  0101  0110  0111  1000  1001  1010  1011  1100  1101  1110  1111
    const warlpiri_roundsinf = [false false false false false false false false  true  true  true  true  true  true  true false  # 0000 ->
                                false false false false false false false false  true  true  true  true  true  true  true  true  # 0001 ->
                                 true  true false false false false false false  true  true  true  true  true  true  true  true  # 0010 ->
                                 true  true  true false false false false false  true  true  true  true  true  true  true  true  # 0011 ->
                                 true  true  true  true false false false false  true  true  true  true  true  true  true  true  # 0100 ->
                                 true  true  true  true  true false false false  true  true  true  true  true  true  true  true  # 0101 ->
                                 true  true  true  true  true  true false false  true  true  true  true  true  true  true  true  # 0110 ->
                                 true  true  true  true  true  true  true false  true  true  true  true  true  true  true  true  # 0111 ->
                                 true  true  true  true  true  true  true  true  true  true  true  true  true  true  true  true  # 1000 ->
                                false false false false false false false false  true false false false false false false false  # 1001 ->
                                false false false false false false false false  true  true false false false false false false  # 1010 ->
                                false false false false false false false false  true  true  true false false false false false  # 1011 ->
                                false false false false false false false false  true  true  true  true false false false false  # 1100 ->
                                false false false false false false false false  true  true  true  true  true false false false  # 1101 ->
                                false false false false false false false false  true  true  true  true  true  true false false  # 1110 ->
                                false false false false false false false false  true  true  true  true  true  true  true false  # 1111 ->
    ]

    @test SigmoidNumbers.roundsinf.(wmatrix) == warlpiri_roundsinf

    #containszero test       ->  0000  0001  0010  0011  0100  0101  0110  0111  1000  1001  1010  1011  1100  1101  1110  1111
    const warlpiri_containszero = [ true  true  true  true  true  true  true  true  true  true  true  true  true  true  true false  # 0000 ->
                                   false false false false false false false false false false false false false false false false  # 0001 ->
                                    true  true false false false false false false false false false false false false false false  # 0010 ->
                                    true  true  true false false false false false false false false false false false false false  # 0011 ->
                                    true  true  true  true false false false false false false false false false false false false  # 0100 ->
                                    true  true  true  true  true false false false false false false false false false false false  # 0101 ->
                                    true  true  true  true  true  true false false false false false false false false false false  # 0110 ->
                                    true  true  true  true  true  true  true false false false false false false false false false  # 0111 ->
                                    true  true  true  true  true  true  true  true false false false false false false false false  # 1000 ->
                                    true  true  true  true  true  true  true  true  true false false false false false false false  # 1001 ->
                                    true  true  true  true  true  true  true  true  true  true false false false false false false  # 1010 ->
                                    true  true  true  true  true  true  true  true  true  true  true false false false false false  # 1011 ->
                                    true  true  true  true  true  true  true  true  true  true  true  true false false false false  # 1100 ->
                                    true  true  true  true  true  true  true  true  true  true  true  true  true false false false  # 1101 ->
                                    true  true  true  true  true  true  true  true  true  true  true  true  true  true false false  # 1110 ->
                                    true  true  true  true  true  true  true  true  true  true  true  true  true  true  true false  # 1111 ->
    ]

    @test SigmoidNumbers.containszero.(wmatrix) == warlpiri_containszero

    #ispositive test          ->  0000  0001  0010  0011  0100  0101  0110  0111  1000  1001  1010  1011  1100  1101  1110  1111
    const warlpiri_ispositive = [false false false false false false false false false false false false false false false false  # 0000 ->
                                 false  true  true  true  true  true  true  true false false false false false false false false  # 0001 ->
                                 false false  true  true  true  true  true  true false false false false false false false false  # 0010 ->
                                 false false false  true  true  true  true  true false false false false false false false false  # 0011 ->
                                 false false false false  true  true  true  true false false false false false false false false  # 0100 ->
                                 false false false false false  true  true  true false false false false false false false false  # 0101 ->
                                 false false false false false false  true  true false false false false false false false false  # 0110 ->
                                 false false false false false false false  true false false false false false false false false  # 0111 ->
                                 false false false false false false false false false false false false false false false false  # 1000 ->
                                 false false false false false false false false false false false false false false false false  # 1001 ->
                                 false false false false false false false false false false false false false false false false  # 1010 ->
                                 false false false false false false false false false false false false false false false false  # 1011 ->
                                 false false false false false false false false false false false false false false false false  # 1100 ->
                                 false false false false false false false false false false false false false false false false  # 1101 ->
                                 false false false false false false false false false false false false false false false false  # 1110 ->
                                 false false false false false false false false false false false false false false false false  # 1111 ->
    ]

    @test SigmoidNumbers.ispositive.(wmatrix) == warlpiri_ispositive

    #isnegative test          ->  0000  0001  0010  0011  0100  0101  0110  0111  1000  1001  1010  1011  1100  1101  1110  1111
    const warlpiri_isnegative = [false false false false false false false false false false false false false false false false  # 0000 ->
                                 false false false false false false false false false false false false false false false false  # 0001 ->
                                 false false false false false false false false false false false false false false false false  # 0010 ->
                                 false false false false false false false false false false false false false false false false  # 0011 ->
                                 false false false false false false false false false false false false false false false false  # 0100 ->
                                 false false false false false false false false false false false false false false false false  # 0101 ->
                                 false false false false false false false false false false false false false false false false  # 0110 ->
                                 false false false false false false false false false false false false false false false false  # 0111 ->
                                 false false false false false false false false false false false false false false false false  # 1000 ->
                                 false false false false false false false false false  true  true  true  true  true  true  true  # 1001 ->
                                 false false false false false false false false false false  true  true  true  true  true  true  # 1010 ->
                                 false false false false false false false false false false false  true  true  true  true  true  # 1011 ->
                                 false false false false false false false false false false false false  true  true  true  true  # 1100 ->
                                 false false false false false false false false false false false false false  true  true  true  # 1101 ->
                                 false false false false false false false false false false false false false false  true  true  # 1110 ->
                                 false false false false false false false false false false false false false false false  true  # 1111 ->
    ]

    @test SigmoidNumbers.isnegative.(wmatrix) == warlpiri_isnegative

    ################################################################################
    #rounds_positive test     ->  0000  0001  0010  0011  0100  0101  0110  0111  1000  1001  1010  1011  1100  1101  1110  1111
    const warlpiri_rounds_pos = [false false false false false false false false  true  true  true  true  true  true  true  true  # 0000 ->
                                 false false false false false false false false false false false false false false false false  # 0001 ->
                                 false false false false false false false false false false false false false false false false  # 0010 ->
                                 false false false false false false false false false false false false false false false false  # 0011 ->
                                 false false false false false false false false false false false false false false false false  # 0100 ->
                                 false false false false false false false false false false false false false false false false  # 0101 ->
                                 false false false false false false false false false false false false false false false false  # 0110 ->
                                 false false false false false false false false false false false false false false false false  # 0111 ->
                                 false false false false false false false false false false false false false false false false  # 1000 ->
                                 false false false false false false false false  true false false false false false false false  # 1001 ->
                                 false false false false false false false false  true  true false false false false false false  # 1010 ->
                                 false false false false false false false false  true  true  true false false false false false  # 1011 ->
                                 false false false false false false false false  true  true  true  true false false false false  # 1100 ->
                                 false false false false false false false false  true  true  true  true  true false false false  # 1101 ->
                                 false false false false false false false false  true  true  true  true  true  true false false  # 1110 ->
                                 false false false false false false false false  true  true  true  true  true  true  true false  # 1111 ->
                                ]

    @test SigmoidNumbers.rounds_positive.(wmatrix) == warlpiri_rounds_pos

    #test equality of all allreals.
    WP_z = zero(WP)
    WP_allreals = [next(a) → a for a in WP if a != WP_z && next(a) != WP_z]
    for a in WP_allreals, b in WP_allreals
        @test a == b
    end

    #test equality of the two null values.
    @test (next(WP_z) → WP_z) == (WP_z → prev(WP_z))
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
function test_kernels()

    N   = rand(50:200)
    M   = rand(50:200)
    dim = rand(5:10)
    Xtrain = rand(N, dim)
    ytrain = rand(N) .>= 0.75
    Xtest  = rand(M, dim)
    ytest = rand(N) .>= 0.75

    kernels = [LinearKernel(),
               SquaredExponentialKernel(),
               RationalQuadraticKernel()]

    @testset "PatMat with $(nameof(typeof(kernel))) kernel" for kernel in kernels
        test_kernelmatrix(PatMat(0.9), kernel, Xtrain, ytrain, Xtest, ytest)
    end

    @testset "PatMatNP with $(nameof(typeof(kernel))) kernel" for kernel in kernels
        test_kernelmatrix(PatMatNP(0.9), kernel, Xtrain, ytrain, Xtest, ytest)
    end

    @testset "TopPushK with $(nameof(typeof(kernel))) kernel" for kernel in kernels
        test_kernelmatrix(TopPushK(5), kernel, Xtrain, ytrain, Xtest, ytest)
    end

    @testset "TopPush with $(nameof(typeof(kernel))) kernel" for kernel in kernels
        test_kernelmatrix(TopPush(), kernel, Xtrain, ytrain, Xtest, ytest)
    end
end

function getmatrix(model::PatMat, Xtrain, ytrain, kernel)
    pos = findall(ytrain)
    nα  = length(pos)
    K   = MLKernels.kernelmatrix(Val(:row), kernel, vcat(Xtrain[pos,:], Xtrain))
    K[1:nα, nα+1:end] .*= -1
    K[nα+1:end, 1:nα] .*= -1
    return K
end

function getmatrix(model::PatMat, Xtrain, ytrain, Xtest, ytest, kernel)
    return getmatrix(model, Xtrain, ytrain, Xtest, kernel)
end

function getmatrix(model::PatMat, Xtrain, ytrain, Xtest, kernel)
    pos = findall(ytrain)
    nα  = length(pos)
    K   = MLKernels.kernelmatrix(Val(:row), kernel, vcat(Xtrain[pos,:], Xtrain), Xtest)
    K[nα+1:end, :] .*= -1
    return K
end

function getmatrix(model::Union{AbstractTopPushK, PatMatNP}, Xtrain, ytrain, kernel)
    pos = findall(ytrain)
    neg = findall(.~ytrain)
    nα  = length(pos)
    K   = MLKernels.kernelmatrix(Val(:row), kernel, vcat(Xtrain[pos,:], Xtrain[neg,:]))
    K[1:nα, nα+1:end] .*= -1
    K[nα+1:end, 1:nα] .*= -1
    return K
end

function getmatrix(model::Union{AbstractTopPushK, PatMatNP}, Xtrain, ytrain, Xtest, ytest, kernel)
    return getmatrix(model, Xtrain, ytrain, Xtest, kernel)
end

function getmatrix(model::Union{AbstractTopPushK, PatMatNP}, Xtrain, ytrain, Xtest, kernel)
    pos = findall(ytrain)
    neg = findall(.~ytrain)
    nα  = length(pos)
    K   = MLKernels.kernelmatrix(kernel, vcat(Xtrain[pos,:], Xtrain[neg,:]), Xtest)
    K[nα+1:end, :] .*= -1
    return K
end

function test_kernelmatrix(model, kernel, Xtrain, ytrain, Xtest, ytest; atol::Real = 1e-10)

    AccuracyAtTopKernels.save_kernelmatrix(model, "train.bin", Xtrain, ytrain; kernel = kernel, T = Float64)
    AccuracyAtTopKernels.save_kernelmatrix(model, "valid.bin", Xtrain, ytrain, Xtest, ytest; kernel = kernel, T = Float64)
    AccuracyAtTopKernels.save_kernelmatrix(model, "test.bin", Xtrain, ytrain, Xtest; kernel = kernel, T = Float64)

    K1  = getmatrix(model, Xtrain, ytrain, kernel)
    K2, = AccuracyAtTopKernels.kernelmatrix(model, Xtrain, ytrain; kernel = kernel)
    t3, io3, out3 = AccuracyAtTopKernels.load_kernelmatrix("train.bin"; T = Float64)
    K3 = out3[1]

    K4  = getmatrix(model, Xtrain, ytrain, Xtest, ytest, kernel)
    K5, = AccuracyAtTopKernels.kernelmatrix(model, Xtrain, ytrain, Xtest, ytest; kernel = kernel)
    t6, io6, out6 = AccuracyAtTopKernels.load_kernelmatrix("valid.bin"; T = Float64)
    K6 = out6[1]

    K7  = getmatrix(model, Xtrain, ytrain, Xtest, kernel)
    K8, = AccuracyAtTopKernels.kernelmatrix(model, Xtrain, ytrain, Xtest; kernel = kernel)
    t9, io9, out9 = AccuracyAtTopKernels.load_kernelmatrix("test.bin"; T = Float64)
    K9 = out9[1]

    @testset "train kernel matrix" begin
        @test K1 ≈ K2 atol = atol
        @test K1 ≈ K3 atol = atol
    end

    @testset "validation kernel matrix" begin
        @test K4 ≈ K5 atol = atol
        @test K4 ≈ K6 atol = atol
    end

    @testset "test kernel matrix" begin
        @test K7 ≈ K8 atol = atol
        @test K7 ≈ K9 atol = atol
    end

    close(io3)
    close(io6)
    close(io9)
    rm("train.bin")
    rm("valid.bin")
    rm("test.bin")
end


text:
using Test, WhittleLikelihoodInference, FiniteDifferences, StaticArrays, FFTW
using ChainRulesTestUtils
import WhittleLikelihoodInference.ChainRulesCore: NoTangent
import WhittleLikelihoodInference:
    nlowertriangle_dimension, nlowertriangle_parameter,
    indexLT, nalias, minbins,
    Acv2EIStorage, Sdf2AcvStorage, Sdf2EIStorage, SdfStorage,
    Acv2EIStorageUni, Sdf2AcvStorageUni, Sdf2EIStorageUni, SdfStorageUni,
    EIstorage_function, EIstorage_gradient, EIstorage_hessian,
    sdfstorage_function, sdfstorage_gradient, sdfstorage_hessian,
    PreallocatedEI, PreallocatedEIGradient, PreallocatedEIHessian,
    PreallocatedSdf, PreallocatedSdfGradient, PreallocatedSdfHessian,
    allocate_memory_EI_F, allocate_memory_EI_FG, allocate_memory_EI_FGH,
    allocate_memory_sdf_F, allocate_memory_sdf_FG, allocate_memory_sdf_FGH,
    AdditiveStorage, extract_asdf,
    CorrelatedOUUnknown, OUUnknown,
    FreqAcvEst, LagsEI, encodetimescale,
    grad_sdf, hess_sdf, grad_acv, hess_acv, EI!, grad_EI!, hess_EI!

## define useful functions
# Function for approximating the gradient using finite differences.
approx_gradient(func, x) = reinterpret(ComplexF64,jacobian(central_fdm(5, 1), func, x)[1])
approx_gradient_uni(func, x) = complex.(permutedims(jacobian(central_fdm(5, 1), func, x)[1]))
# Function for approximating the hessian using finite differences.
function approx_hessian(grad, x)
    hess = zeros(ComplexF64, size(grad(x),1), length(x)*(length(x)+1)÷2)
    count = 1
    @views for i ∈ 1:length(x)
        hess[:,count:count+length(x)-i] = approx_gradient(y->grad(y)[:,i], x)[:,i:end]
        count += length(x)-i+1
    end
    return hess
end
function approx_hessian_uni(grad, x)
    full_hess = approx_gradient(grad, x)
    hess = ones(ComplexF64, size(full_hess,1)*(size(full_hess,1)+1)÷2)
    for i ∈ 1:size(full_hess,1), j ∈ 1:i
        hess[WhittleLikelihoodInference.indexLT(i,j,size(full_hess,1))] = full_hess[i,j]
    end
    return hess
end

## run tests
include("typestructure_test.jl")
include("memoryallocation_test.jl")
include("properties/spectraldensity_test.jl")
include("properties/autocovariance_test.jl")
include("properties/expectedperiodogram_test.jl")
include("properties/coherancy_test.jl")

include("whittle/whittledata_test.jl")
include("whittle/generalwhittle_test.jl")
include("whittle/standardwhittle_test.jl")
include("whittle/debiasedwhittle_test.jl")
include("whittle/chainrules_test.jl")
include("fit_test.jl")

include("models/OU_test.jl")
include("models/CorrelatedOU_test.jl")
include("models/OUUnknown_test.jl")
include("models/CorrelatedOUUnknown_test.jl")
include("models/Matern_test.jl")
include("models/MaternAcv_test.jl")

include("simulation_test.jl")
include("nonparametric_test.jl")
include("plotting_test.jl")


text:
const sample_stats_types = Dict(
    "mean_tree_accept" => Float64,
    "energy" => Float64,
    "energy_error" => Float64,
    "max_energy_error" => Float64,
    "step_size" => Float64,
    "step_size_bar" => Float64,
    "tree_size" => Int,
    "depth" => Int,
    "tune" => Bool,
    "diverging" => Bool,
)

@doc forwarddoc(:compare) compare(args...; kwargs...) =
    arviz.compare(args...; kwargs...) |> Pandas.DataFrame

Docs.getdoc(::typeof(compare)) = forwardgetdoc(:compare)

@forwardfun hpd

@doc forwarddoc(:loo) loo(args...; kwargs...) =
    arviz.loo(args...; kwargs...) |> Pandas.Series

Docs.getdoc(::typeof(loo)) = forwardgetdoc(:loo)

@forwardfun loo_pit

@forwardfun psislw

@doc forwarddoc(:r2_score) r2_score(args...; kwargs...) =
    arviz.r2_score(args...; kwargs...) |> Pandas.Series

Docs.getdoc(::typeof(r2_score)) = forwardgetdoc(:r2_score)

@doc forwarddoc(:waic) waic(args...; kwargs...) =
    arviz.waic(args...; kwargs...) |> Pandas.Series

Docs.getdoc(::typeof(waic)) = forwardgetdoc(:waic)

"""
    summarystats(data::Dataset; kwargs...) -> Union{Pandas.DataFrame,Dataset}
    summarystats(data::InferenceData; group = :posterior kwargs...) -> Union{Pandas.DataFrame,Dataset}

Compute summary statistics on `data`.

# Arguments
- `data::Union{Dataset,InferenceData}`: The data on which to compute summary
      statistics. If `data` is an [`InferenceData`](@ref), only the dataset
      corresponding to `group` is used.

# Keywords
- `var_names::Vector{String}=nothing`: Names of variables to include in summary
- `include_circ::Bool=false`: Whether to include circular statistics
- `fmt::String="wide"`: Return format is either `Pandas.DataFrame` ("wide", "long")
      or [`Dataset`](@ref) ("xarray").
- `round_to::Int=nothing`: Number of decimals used to round results.
      Use `nothing` to return raw numbers.
- `stat_funcs::Union{Dict{String,Function},Vector{Function}}=nothing`:
      A vector of functions or a dict of functions with function names as keys
      used to calculate statistics. By default, the mean, standard deviation,
      simulation standard error, and highest posterior density intervals are
      included.
      The functions will be given one argument, the samples for a variable as an
      array, The functions should operate on an array, returning a single
      number. For example, `Statistics.mean`, or `Statistics.var` would both
      work.
- `extend::Bool=true`: If `true`, use the statistics returned by `stat_funcs` in
      addition to, rather than in place of, the default statistics. This is only
      meaningful when `stat_funcs` is not `nothing`.
- `credible_interval::Real=0.94`: Credible interval to plot. This is only
      meaningful when `stat_funcs` is `nothing`.
- `order::String="C"`: If `fmt` is "wide", use either "C" or "F" unpacking order.
- `index_origin::Int=1`: If `fmt` is "wide", select 𝑛-based indexing for
      multivariate parameters.

# Returns
- `Union{Pandas.DataFrame,Dataset}`: Return type dicated by `fmt` argument.
      Return value will contain summary statistics for each variable. Default
      statistics are:
    + `mean`
    + `sd`
    + `hpd_3%`
    + `hpd_97%`
    + `mcse_mean`
    + `mcse_sd`
    + `ess_bulk`
    + `ess_tail`
    + `r_hat` (only computed for traces with 2 or more chains)

# Examples

```@example summarystats
using ArviZ
idata = load_arviz_data("centered_eight")
summarystats(idata; var_names=["mu", "tau"])
```

Other statistics can be calculated by passing a list of functions or a
dictionary with key, function pairs:

```@example summarystats
using StatsBase, Statistics
function median_sd(x)
    med = median(x)
    sd = sqrt(mean((x .- med).^2))
    return sd
end

func_dict = Dict(
    "std" => x -> std(x; corrected = false),
    "median_std" => median_sd,
    "5%" => x -> percentile(x, 5),
    "median" => median,
    "95%" => x -> percentile(x, 95),
)

summarystats(idata; var_names = ["mu", "tau"], stat_funcs = func_dict, extend = false)
```
"""
function StatsBase.summarystats(data::Dataset; index_origin = 1, kwargs...)
    s = arviz.summary(data; index_origin = index_origin, kwargs...)
    s isa Dataset && return s
    return Pandas.DataFrame(s)
end

function StatsBase.summarystats(data::InferenceData; group = :posterior, kwargs...)
    dataset = getproperty(data, Symbol(group))
    return summarystats(dataset; kwargs...)
end

"""
    summary(
        data;
        group = :posterior,
        coords = nothing,
        dims = nothing,
        kwargs...,
    ) -> Union{Pandas.DataFrame,PyObject}

Compute summary statistics on any object that can be passed to
[`convert_to_dataset`](@ref).

# Keywords
- `coords::Dict{String,Vector}=nothing`: Map from named dimension to named
      indices.
- `dims::Dict{String,Vector{String}}=nothing`: Map from variable name to names
      of its dimensions.
- `kwargs`: Keyword arguments passed to [`summarystats`](@ref).
"""
function summary(data; group = :posterior, coords = nothing, dims = nothing, kwargs...)
    idata = convert_to_inference_data(data; group = group, coords = coords, dims = dims)
    return summarystats(idata; group = group, kwargs...)
end


text:
@testset  "select" begin
    values_ = [:a, 1.0, "hello", 2]
    time_ = [0, 2, 3, 5]

    ts = EventSeries(time_, values_)

    ts1 = select(ts, 0, 5)
    @test ts1.timestamps == time_
    @test ts1.values == values_

    ts2 = select(ts, 1, 4)
    @test ts2.timestamps == [1, 2, 3, 4]
    @test ts2.values ==  [:a, 1.0, "hello", "hello"]

    ts3 = select(ts, 0, 4)
    @test ts3.timestamps == [0, 2, 3, 4]
    @test ts3.values ==  [:a, 1.0, "hello", "hello"]

    ts4 = select(ts, 2, 5)
    @test ts4.timestamps == [2, 3, 5]
    @test ts4.values ==  [1.0, "hello", 2]

    @test_throws AssertionError select(ts, -1, 4)
    @test_throws AssertionError select(ts, 1, 10)
    @test_throws AssertionError select(ts, -1, 10)

    values_ = [:a, :b]
    time_ = [0, 10]
    ts = EventSeries(time_, values_)
    ts5 = select(ts, 1,2)
    @test ts5.timestamps == [1,2]
    @test ts5.values ==  [:a, :a]
end

@testset  "align" begin
    t1_ = [0, 2, 4]
    v1_ = [1, 4, 5]

    t2_ = [1, 2, 3]
    v2_ = [4, 3, 6]

    ts1 = EventSeries(t1_, v1_)
    ts2 = EventSeries(t2_, v2_)

    ts1_, = align(ts1)
    @test ts1_.values == ts1.values
    @test ts1_.timestamps == ts1.timestamps

    ts1_, ts2_ = align(ts1, ts2)

    @test ts1_.timestamps == [1,2,3]
    @test ts1_.values == [1,4,4]

    @test ts2_.timestamps == [1,2,3]
    @test ts2_.values == v2_
end

@testset  "cumtime" begin
    t_ = 1:10 |> collect
    v_ = repeat([1,2], 5)
    ts = EventSeries(t_, v_)
    @test cumtime(ts, val->val==1) == 5
    @test cumtime(ts, val->val==2) == 4 # integration stops at last time stamp
end

@testset  "fuse and splice" begin
    t1 = [1,3,5]
    v1 = ['A', 'B', 'C']
    ts1 = EventSeries(t1, v1)

    t2 = [2,4,6]
    v2 = ['a', 'b', 'c']
    ts2 = EventSeries(t2, v2)

    ts3 = fuse(x=ts1, y=ts2)
    ts4 = splice(x=ts1, y=ts2)

    @test ts3.timestamps == sort([t1; t2])
    @test ts4.timestamps == sort([t1; t2])

    @test ts3.values == [(x = 'A', y = nothing)
                         (x = 'A', y = 'a')
                         (x = 'B', y = 'a')
                         (x = 'B', y = 'b')
                         (x = 'C', y = 'b')
                         (x = 'C', y = 'c')]

    @test ts4.values == [:x => 'A',
                         :y => 'a',
                         :x => 'B',
                         :y => 'b',
                         :x => 'C',
                         :y => 'c']
end

@testset  "segments" begin
    t = [1, 3, 5]
    v = ['A', 'B', 'C']
    ts = EventSeries(t, v)
    @test collect(segments(ts)) == [Segment('A', 1,3), Segment('B', 3,5)]
end

@testset  "filter" begin
    t = [1, 3, 5]
    v = ['A', 'B', 'C']
    ts = EventSeries(t, v)
    ts2 = filter(x->x.value!='C', ts)
    @test ts2.values == ['A', 'B']
    @test ts2.timestamps == [1, 3]

end

@testset "fuse duplicate timestamps" begin
    t = [1, 3, 5]
    v1 = ['A', 'B', 'C']
    v2 = ['a', 'b', 'c']
    ts1 = EventSeries(t, v1)
    ts2 = EventSeries(t, v2)
    ts = fuse(a=ts1, b=ts2)
    @test length(ts) == length(t)
end


text:
module Day12

import ..data_dir # from parent module
input = read(joinpath(data_dir, "day12"), String)

export part1, part2

struct Graph
    vertices::Vector{String}
    edges::Vector{Tuple{String,String}}
end

function Graph(input::AbstractString)
    raw_edges = string.(split(input, "\n"))
    edges = Tuple{String,String}.(split.(raw_edges, "-"))
    vertices = String.(vcat(split.(raw_edges, "-")...))
    Graph(vertices, edges)
end

function neighbours(graph::Graph, vertex::AbstractString)
    result = []
    for edge in graph.edges
        x, y = edge
        if x == vertex
            push!(result, y)
        elseif y == vertex
            push!(result, x)
        end
    end
    result
end

cave_size(vertex::String) = uppercase(vertex) == vertex ? "large" : "small"

function has_visited_a_small_cave_twice(path::Vector{String})
    small_caves = filter(vertex -> lowercase(vertex) == vertex, path)
    unique(small_caves) != small_caves
end

function traverse(
    graph::Graph,
    path_so_far = Vector{String}(),
    next_vertex = "start";
    allow_small_cave_revisit = false
)
    updated_path = [path_so_far; next_vertex]
    if next_vertex == "end"
        return Set([updated_path])
    end
    if allow_small_cave_revisit && has_visited_a_small_cave_twice(updated_path)
        return traverse(graph, path_so_far, next_vertex; allow_small_cave_revisit = false)
    end
    valid_neighbours = filter(candidate -> candidate != "start", neighbours(graph, next_vertex))
    if !allow_small_cave_revisit
        valid_neighbours = filter(
            candidate -> cave_size(candidate) == "large" || candidate ∉ updated_path,
            valid_neighbours
        )
        if length(valid_neighbours) == 0 # dead end
            return Set{Vector{String}}()
        end
    end
    branches = [
        traverse(graph, updated_path, neighbour; allow_small_cave_revisit = allow_small_cave_revisit)
        for neighbour in valid_neighbours
    ]
    return ∪(branches...)
end

function part1(input = input)
    traverse(
        Graph(input);
        allow_small_cave_revisit = false
    ) |> length
end

function part2(input = input)
    traverse(
        Graph(input);
        allow_small_cave_revisit = true
    ) |> length
end

end #module


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Tuple{Core.kwftype(typeof(optimize_ticks)),NamedTuple{(:k_min, :k_max),Tuple{Int,Int}},typeof(optimize_ticks),Float64,Float64})
    for C in (RGB, RGBA), T in (Colors.FixedPointNumbers.N0f8, Float32, Float64)
        @assert precompile(cgrad, (ColorScheme{Vector{C{T}},String,String},))
        for V in (Vector{Float64}, typeof(get_range(3)))
            @assert precompile(cgrad, (ColorScheme{Vector{C{T}},String,String}, V))
            @assert precompile(CategoricalColorGradient, (ColorScheme{Vector{C{T}},String,String}, V))
            @assert precompile(ContinuousColorGradient, (ColorScheme{Vector{C{T}},String,String}, V))
        end
    end
    for T in (Float64, Int)
        @assert precompile(zscale, (Vector{T},))
    end
    @assert precompile(optimize_ticks, (Int, Int))
    @assert precompile(optimize_datetime_ticks, (Int, Int))
    @assert precompile(adapted_grid, (Function, Tuple{Int,Int}))
end


text:
# =========================================================================================
"""
    writeskel{T, M <: Model{T}}(
        stream::IOStream,
        model ::M
    )

Creates a SKEL file from a given surface or volume model, representing the model as a
collection of points and polylines.

# Specification
<http://www.geomview.org/docs/html/SKEL.html>

# Return type
`Void`

# Alias

    writeskel{T, M <: Model{T}}(
        fname::String,
        model::M
    )

Creates the SKEL file by name rather than `IOStream` object.
"""
function writeskel(
        stream::IOStream,
        model ::Model{T, Triangle{T}}
    ) where T
    println(stream, "SKEL")
    println(stream, "$(length(model.nodes))\t$(length(model.elements))")
    for node in model.nodes
        println(stream, "$(node[1])\t$(node[2])\t$(node[3])")
    end
    revidx = reverseindex(model.nodes)
    for elem in model.elements
        (v1, v2, v3) = (revidx[objectid(elem.v1)] - 1,
                        revidx[objectid(elem.v2)] - 1,
                        revidx[objectid(elem.v3)] - 1)
        println(stream, "4\t$v1\t$v2\t$v3\t$v1")
    end
    nothing
end

function writeskel(
        stream::IOStream,
        model ::Model{T, Tetrahedron{T}}
    ) where T
    println(stream, "SKEL")
    println(stream, "$(length(model.nodes))\t$(length(model.elements))")
    for node in model.nodes
        println(stream, "$(node[1])\t$(node[2])\t$(node[3])")
    end
    revidx = reverseindex(model.nodes)
    for elem in model.elements
        (v1, v2, v3, v4) = (revidx[objectid(elem.v1)] - 1,
                            revidx[objectid(elem.v2)] - 1,
                            revidx[objectid(elem.v3)] - 1,
                            revidx[objectid(elem.v4)] - 1)
        println(stream, "7\t$v1\t$v2\t$v3\t$v4\t$v2\t$v4\t$v1")
    end
    nothing
end

function writeskel(
        fname::String,
        model::M
    ) where {T, M <: Model{T}}
    open(fh -> writeskel(fh, model), fname, "w")
end


text: negbindata = [0, 1, 4, 0, 2, 2, 5, 0, 1]


text:
LarA = Lar.Arrangement
LAR = LinearAlgebraicRepresentation
#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 1 METHODS
#    frag_edge_channel
#    frag_edge
#    intersect_edges
#    merge_vertices!
#-------------------------------------------------------------------------------

"""                                                                             ## UNTESTED
    frag_edge_channel(in_chan, out_chan, model::Lar.Model, bigPI)

Utility function for parallel edge fragmentation.

This function handles the edge fragmentation in the first part of arrangement's
algorithmic pipeline (see also [`Lar.Arrangement.planar_arrangement_1`](@ref))
if multiprocessing computation has been enabled.
In order to do so, it needs two `Distributed.RemoteChannel`s, one with
the inputs and one for outputs.

See also: [`Lar.Arrangement.planar_arrangement_1`](@ref)
It uses: [`Lar.Arrangement.frag_edge`](@ref)
"""
function frag_edge_channel(
        in_chan::Distributed.RemoteChannel{Channel{Int64}},
        out_chan::Distributed.RemoteChannel{Channel{Tuple}},
        model::Lar.Model,
        bigPI::Array{Array{Int64,1},1}
    )
    run_loop = true
    while run_loop
        edgenum = take!(in_chan)
        if edgenum != -1
            put!(out_chan, (edgenum, frag_edge(model, edgenum, bigPI)))
        else
            run_loop = false
        end
    end
end

"""
    frag_edge(
        model::Lar.Model,
        edge_idx::Int,
        bigPI::Array{Array{Int64,1},1}
    )::Tuple{Lar.Points, Lar.ChainOp}

Splits the `edge_idx`-th edge of the Edge Topology `model.T[1]`.

This method splits the `edge_idx`-th edge between the Edges in the Topology
`EV = model.T[1]` into several parts by confronting it with the others that
intersect its bounding box (see also [`Lar.Arrangement.intersect_edges`](@ref)).
Returns a set made of the detected vertices over the edge (with redundancies)
and the associated cochain (without redundancies).

See also:
 - [`Lar.Arrangement.planar_arrangement_1`](@ref)
 - [`Lar.Arrangement.frag_edge_channel`](@ref)
It uses: [`Lar.Arrangement.intersect_edges`](@ref)
---
# Examples
```jldoctest
julia> model = Lar.Model(hcat([
    [1.0, 0.0], [0.0, 1.0], [0.0, 0.5], [0.5, 1.0], [1.0, 1.0]
]...));
julia> Lar.addModelCells!(model, 1, [[1, 2], [2, 5], [3, 4], [4, 5]])
julia> bigPI = LarA.spaceIndex(model);
julia> LarA.frag_edge(model, 1, bigPI)[1]
4×2 Array{Float64,2}:
 1.0   0.0
 0.0   1.0
 0.25  0.75
 0.0   1.0
julia> Matrix(LarA.frag_edge(model, 1, bigPI)[2])
2×4 Array{Int8,2}:
 1  0  1  0
 0  1  1  0
```
"""
function frag_edge(
        model::Lar.Model,
        edge_idx::Int,
        bigPI::Array{Array{Int64,1},1}
    )::Tuple{Lar.Points, Lar.ChainOp}
    V = convert(Lar.Points, model.G')
    EV = model.T[1]
    alphas = Dict{Float64, Int}()
    edge = EV[edge_idx, :]
    verts = V[edge.nzind, :]
    for i in bigPI[edge_idx]
        if i != edge_idx # && edge_idx in bigPI[i]
            intersection = LarA.intersect_edges(V, edge, EV[i, :])
            for (point, alpha) in intersection
                verts = [verts; point]
                alphas[alpha] = size(verts, 1)
            end
        end
    end
    alphas[0.0], alphas[1.0] = [1, 2]
    alphas_keys = sort(collect(keys(alphas)))
    edge_num = length(alphas_keys)-1
    verts_num = size(verts, 1)
    ev = SparseArrays.spzeros(Int8, edge_num, verts_num)
    for i in 1:edge_num
        ev[i, alphas[alphas_keys[i]]] = 1
        ev[i, alphas[alphas_keys[i+1]]] = 1
    end
    return verts, ev
end

"""
    intersect_edges(
        V::Lar.Points,
        edge1::Lar.Cell,
        edge2::Lar.Cell
    )::Array{Tuple{Lar.Points, Float64}, 1}

Finds the intersection points (if there exist) between the two given edges.

Compute the new points over `edge1` given by the intersection with `edge2`.
For each intersection point it evaluates the scalar corresponding to the
relative distance from the first vertex of the edge.
Note that vertices of `edge1` are therefore never detected (see example 2).

See also: [`Lar.Arrangement.frag_edge`](@ref)
---
# Examples
```jldoctest
# Cross
julia> V = [1.0 0.0; 0.0 1.0; 0.0 0.5; 0.5 1.0];
julia> copEV = LAR.coboundary_0([[1, 2], [3, 4]]);
julia> LarA.intersect_edges(V, copEV[1, :], copEV[2, :])
1-element Array{Tuple{Array{T,2} where T,Float64},1}:
 ([0.25 0.75], 0.75)
```

```jldoctest
# Collinear
julia> V = [1.0 0.0; 0.0 1.0; 0.75 0.25; 0.5 0.5];
julia> copEV = LAR.coboundary_0([[1, 2], [3, 4]]);
julia> LarA.intersect_edges(V, copEV[1, :], copEV[2, :])
2-element Array{Tuple{Array{T,2} where T,Float64},1}:
 ([0.75 0.25], 0.25)
 ([0.5 0.5], 0.5)
julia> LarA.intersect_edges(V, copEV[2, :], copEV[1, :])
0-element Array{Tuple{Array{T,2} where T,Float64},1}
```
"""
function intersect_edges(
        V::Lar.Points,
        edge1::Lar.Cell,
        edge2::Lar.Cell
    )::Array{Tuple{Lar.Points, Float64}, 1}

    err = 10e-8

    x1, y1, x2, y2 = vcat(map(c->V[c, :], edge1.nzind)...)
    x3, y3, x4, y4 = vcat(map(c->V[c, :], edge2.nzind)...)
    ret = Array{Tuple{Lar.Points, Float64}, 1}()

    v1 = [x2-x1, y2-y1]
    v2 = [x4-x3, y4-y3]
    v3 = [x3-x1, y3-y1]
    ang1 = dot(LA.normalize(v1), LA.normalize(v2))
    ang2 = dot(LA.normalize(v1), LA.normalize(v3))
    parallel = 1-err < abs(ang1) < 1+err
    colinear = parallel && (1-err < abs(ang2) < 1+err || -err < norm(v3) < err)
    if colinear
        o = [x1 y1]
        v = [x2 y2] - o
        alpha = 1/dot(v,v')
        ps = [x3 y3; x4 y4]
        for i in 1:2
            a = alpha*dot(v',(reshape(ps[i, :], 1, 2)-o))
            if 0 < a < 1
                push!(ret, (ps[i:i, :], a))
            end
        end
    elseif !parallel
        denom = (v2[2])*(v1[1]) - (v2[1])*(v1[2])
        a = ((v2[1])*(-v3[2]) - (v2[2])*(-v3[1])) / denom
        b = ((v1[1])*(-v3[2]) - (v1[2])*(-v3[1])) / denom

        if -err < a < 1+err && -err <= b <= 1+err
            p = [(x1 + a*(x2-x1))  (y1 + a*(y2-y1))]
            push!(ret, (p, a))
        end
    end
    return ret
end

"""
merge_vertices!(
    model::Lar.model,
    [ edge_map::Array{Array{Int64,1},1} = Array{Array{Int64,1},1}() ],
    [ err::Float64 = 1e-4 ]
)::Nothing

Modify `model` by compacting its points closer than `err` in a single one.

This method check one at time each vertex ``v`` in `model.G` and identifies
each other vertex within `err` with ``v`` itself.
The edges cochain `model.T[1]` is coherently modified (multiple edges between
two vertices are not allowed).
If an `edge_map` is given in input (this could be usefull during the planar
arrangement), then also the map is coherently modified and given back in output.

See also: [`Lar.Arrangement.planar_arrangement_1`](@ref)
---
# Examples
```jldoctest
julia> model = Lar.Model([
    0.5 0.0 0.5 1.0 0.5 1.0
    0.5 0.0 0.5 1.0 0.5 1.0
]);
julia> Lar.addModelCells!(model, 1, [[1, 4], [3, 2], [5, 6], [1, 6], [5, 3]]);
julia> LarA.merge_vertices!(model)

julia> model.G
2×3 Array{Float64,2}:
 0.5  0.0  1.0
 0.5  0.0  1.0

julia> Matrix(model.T[1])
2×3 Array{Int8,2}:
 1  0  1
 1  1  0
```
"""
function merge_vertices!(
        model::Lar.Model,
        edge_map::Array{Array{Int64,1},1} = Array{Array{Int64,1},1}(),
        err::Float64 = 1e-4
    )::Nothing
    V = convert(Lar.Points, model.G')
    EV = model.T[1]
    vertsnum = size(V, 1)
    edgenum = size(EV, 1)
    newverts = zeros(Int, vertsnum)
    # KDTree constructor needs an explicit array of Float64
    V = Array{Float64,2}(V)
    kdtree = KDTree(permutedims(V))

    # merge congruent vertices
    todelete = []
    i = 1
    for vi in 1:vertsnum
        if !(vi in todelete)
            nearvs = LAR.inrange(kdtree, V[vi, :], err)
            newverts[nearvs] .= i
            nearvs = setdiff(nearvs, vi)
            todelete = union(todelete, nearvs)
            i = i + 1
        end
    end
    nV = V[setdiff(collect(1:vertsnum), todelete), :]

    # merge congruent edges
    edges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    oedges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    for ei in 1:edgenum
        v1, v2 = EV[ei, :].nzind
        edges[ei] = Tuple{Int, Int}(sort([newverts[v1], newverts[v2]]))
        oedges[ei] = Tuple{Int, Int}(sort([v1, v2]))
    end
    nedges = union(edges)
    nedges = filter(t->t[1]!=t[2], nedges)
    nedgenum = length(nedges)
    nEV = spzeros(Int8, nedgenum, size(nV, 1))
    # maps pairs of vertex indices to edge index
    etuple2idx = Dict{Tuple{Int, Int}, Int}()
    # builds `edge_map`
    for ei in 1:nedgenum
        nEV[ei, collect(nedges[ei])] .= 1
        etuple2idx[nedges[ei]] = ei
    end
    if !isempty(edge_map)
        for i in 1:length(edge_map)
            row = edge_map[i]
            row = map(x->edges[x], row)
            row = filter(t->t[1]!=t[2], row)
            row = map(x->etuple2idx[x], row)
            edge_map[i] = row
        end
    end
    # return new vertices and new edges
    model.G = convert(Lar.Points, nV')
    model.T[1] = nEV
    model.T[2] = SparseArrays.spzeros(Int8, 0, size(model, 1, 1))
    return
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - CLEAN DECOMPOSITION
#-------------------------------------------------------------------------------

"""
    cleandecomposition(
        model::Lar.Model,
        sigma::Union{Lar.Chain, Array{Int,1}},
        [ edge_map::Array{Array{Int64,1},1} = Array{Array{Int64,1},1}() ]
    )

This function clears the `model` from all edges outside ``σ``.

This function take an arranged-edge `model` and a `Lar.Chain` ``σ`` made of
edges indices of the 1-chains `model.T[1]`.
The method drops from the Topology all the edges that are outside ``σ`` and
gives back the dropped model.
If an `edge_map` is given as an input it is also purged from the edges outside
w.r.t. ``σ`` and it is given as an output as well.

See also: [`Lar.planar_arrangement`](@ref)
---
# Examples
```jldoctest
# Nested Triangles
julia> model = Lar.Model([
    0.0  2.0  4.0  1.0  3.0  2.0  2.0 -2.0  6.0
    0.0  0.0  0.0  1.5  1.5  3.0 -3.0  3.0  3.0
]);
julia> Lar.addModelCells!(model, 1, [
    [1, 3], [1, 6], [3, 6], [2, 4], [2, 5], [4, 5], [7, 8], [7, 9], [8, 9]
]);

julia> σ = SparseArrays.sparse([1; 1; 1; 0; 0; 0; 0; 0; 0]);

julia> model = LarA.cleandecomposition(model, convert(Lar.Chain, σ))

julia> model.G
2×6 Array{Float64,2}:
 0.0  2.0  4.0  1.0  3.0  2.0
 0.0  0.0  0.0  1.5  1.5  3.0

julia> Matrix(model.T[1])
6×6 Array{Int8,2}:
 1  0  1  0  0  0
 1  0  0  0  0  1
 0  0  1  0  0  1
 0  1  0  1  0  0
 0  1  0  0  1  0
 0  0  0  1  1  0
```
"""
function cleandecomposition(
        model::Lar.Model,
        sigma::Union{Lar.Chain, Array{Int,1}},
        edge_map::Array{Array{Int64,1},1} = Array{Array{Int64,1},1}()
    )::Union{Lar.Model, Tuple{Lar.Model, Array{Array{Int64,1},1}}}

    model = deepcopy(model)
    edge_map = deepcopy(edge_map)
    # Model point extraction to use LAR.point_in_face                           ##
    V = convert(Lar.Points, model.G')

    if issparse(sigma)
        sigma = sigma.nzind
    end
    todel = Array{Int,1}()
    sigma_edges = model.T[1][sigma, :]
    for e in 1 : model.T[1].m
        # Since the model is already arranged an edge can only be a sigma-edge,
        #  an intersection-free innner edge w.r.t. sigma face (not to purge),
        #  or an intersection-free outer edge (has to be deleted).
        if !(e in sigma)
            v1, v2 = Lar.getModelCellVertices(model, 1, e)
            centroid = .5*(v1 + v2)

            if ! LAR.point_in_face(centroid, V, sigma_edges)
                push!(todel, e)
            end
        end
    end

    Lar.deleteModelCells!(model, 1, todel)
    Lar.modelPurge!(model)

    if isempty(edge_map)
        return model
    end

    # Edges in edge_map must be updated too according to the same approach
    #  edges are deleted from models.
    for i in reverse(todel)
        for row in edge_map

            filter!(x->x!=i, row)

            for j in 1:length(row)
                if row[j] > i
                    row[j] -= 1
                end
            end
        end
    end
	return model, edge_map
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - BICONNECTED COMPONENTS
#-------------------------------------------------------------------------------

"""
    biconnected_components(EV::Lar.ChainOp)

Compute sets of 2-cells on the same 3-cells biconnected components.

The method evaluates the 2-cells ``σ_i`` wich lies on the same 3-cells
biconnected component ``χ_j`` and gives back an array that contains an array
of ``σ_i`` for each ``j``.
Do note that if the cochain `EV` contains more copies of the same 2-cell then
it will be considered like a 3-cell.

Intermediate part of Planar Arrangement's algorithmic pipeline.

See also:
 - [`Lar.planar_arrangement`](@ref) for the complete pipeline.
 - [`Lar.Arragnement.planar_arrangement_2`](@ref).
---

# Examples
```jldoctest
julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
        [1 1 0 0 0 0] #1 -> 1,2  |
        [1 0 1 0 0 0] #2 -> 1,3  |
        [1 0 0 1 0 0] #3 -> 1,4   |
        [1 0 0 0 1 0] #4 -> 1,5   |
        [1 0 0 0 0 1] #5 -> 1,6
        [0 1 1 0 0 0] #6 -> 2,3  |
        [0 0 0 1 1 0] #7 -> 4,5   |
    ]));
julia> LarA.biconnected_components(copEV)
2-element Array{Array{Int64,1},1}:
 [2, 6, 1]
 [4, 7, 3]
```

```jldoctest
julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
        [1 1 0] #1 -> 1,2  |
        [1 1 0] #2 -> 1,2  |
        [1 0 1] #3 -> 1,2
    ]));
julia> LarA.biconnected_components(copEV)
1-element Array{Array{Int64,1},1}:
 [2, 1]
```
"""
function biconnected_components(EV::Lar.ChainOp)

    ps = Array{Tuple{Int, Int, Int}, 1}()
    es = Array{Tuple{Int, Int}, 1}()
    todel = Array{Int, 1}()
    visited = Array{Int, 1}()
    bicon_comps = Array{Array{Int, 1}, 1}()
    hivtx = 1

    function an_edge(point) # TODO: fix bug
        # error? : BoundsError: attempt to access 0×0 SparseMatrix ...
        edges = setdiff(EV[:, point].nzind, todel)
        if length(edges) == 0
            edges = [false]
        end
        edges[1]
    end

    function get_head(edge, tail)
        setdiff(EV[edge, :].nzind, [tail])[1]
    end

    function v_to_vi(v)
        i = findfirst(t->t[1]==v, ps)
        # seems findfirst changed from 0 to Nothing
        if typeof(i) == Nothing
            return false
        elseif i == 0
            return false
        else
            return ps[i][2]
        end
    end

    push!(ps, (1,1,1))
    push!(visited, 1)
    exit = false
    while !exit
        edge = an_edge(ps[end][1])
        if edge != false
            tail = ps[end][2]
            head = get_head(edge, ps[end][1])
            hi = v_to_vi(head)
            if hi == false
                hivtx += 1
                push!(ps, (head, hivtx, ps[end][2]))
                push!(visited, head)
            else
                if hi < ps[end][3]
                    ps[end] = (ps[end][1], ps[end][2], hi)
                end
            end
            push!(es, (edge, tail))
            push!(todel, edge)
        else
            if length(ps) == 1
                found = false
                pop!(ps)
                for i in 1:size(EV,2)
                    if !(i in visited)
                        hivtx = 1
                        push!(ps, (i, hivtx, 1))
                        push!(visited, i)
                        found = true
                        break
                    end
                end
                if !found
                    exit = true
                end

            else
                if ps[end][3] == ps[end-1][2]
                    edges = Array{Int, 1}()
                    while true
                        edge, tail = pop!(es)
                        push!(edges, edge)
                        if tail == ps[end][3]
                            if length(edges) > 1
                                push!(bicon_comps, edges)
                            end
                            break
                        end
                    end

                else
                    if ps[end-1][3] > ps[end][3]
                        ps[end-1] = (ps[end-1][1], ps[end-1][2], ps[end][3])
                    end
                end
                pop!(ps)
            end
        end
    end
    bicon_comps = sort(bicon_comps, lt=(x,y)->length(x)>length(y))
    return bicon_comps
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - REMOVE MAPPING DANGLING EDGES
#-------------------------------------------------------------------------------

"""
    remove_mapping_dangling_edges(
        edge_map::Array{Array{Int64,1},1},
        bicon_comps::Array{Array{Int64,1},1}
    )::Array{Array{Int64,1},1}

Remove from the edge_map the edges not belonging to a biconnected component.

This utility function deletes from the `edge_map` all the edges that do not
belongs to a biconnected component by ordinatelly recompatting the indices of
the remaning edges.

See also: [`Lar.Arrangement.planar_arrangement`](@ref).
"""
function remove_mapping_dangling_edges(
        edge_map::Array{Array{Int64,1},1},
        bicon_comps::Array{Array{Int64,1},1}
    )::Array{Array{Int64,1},1}

    edge_map = deepcopy(edge_map)

    edges = sort(union(bicon_comps...))
                                  #size(model.T[1, 1]) == max(edge_map...)
    todel = sort(setdiff(collect(1:max(max(edge_map...)...)), edges))

    for i in reverse(todel)
        for row in edge_map

            filter!(x->x!=i, row)

            for j in 1:length(row)
                if row[j] > i
                    row[j] -= 1
                end
            end
        end
    end

    return edge_map
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 2 - COMPONENT GRAPH (TGW) METHODS
#    get_external_cycle
#    pre_containment_test
#    prune_containment_graph
#    transitive_reduction!
#-------------------------------------------------------------------------------

"""
    get_external_cycle(model::Lar.Model)::Union{Int64, Nothing}
    get_external_cycle(V::Lar.Points, EV::Lar.ChainOp, FE::Lar.ChainOp)

Evaluates the index of the external 3D-cell of a Topological Model.

This method looks for and retrieve the external cycle of the `model`.
If the cell does not exist then it returns `nothing`.
The method is also callable by the signature `(V, EV, FE)`.

See also: [`Lar.Arrangement.componentgraph`](@ref).
---
# Examples
```jldoctest
# ``K^4`` graph
julia> model = Lar.Model(hcat([[0.0,0.0], [4.0,0.0], [2.0,3.0], [2.0,1.5]]...));
julia> Lar.addModelCells!(model, 1, [[1,2],[2,3],[3,1],[1,4],[2,4],[3,4]]);
julia> Lar.addModelCells!(model, 2, [[1,4,5],[2,5,6],[3,4,6],[1,2,3]]);

julia> LarA.get_external_cycle(model)
4

julia> Lar.deleteModelCell!(model, 2, 4);
julia> typeof(LarA.get_external_cycle(model))
Nothing
```
"""
function get_external_cycle(model::Lar.Model)::Union{Int64, Nothing}
    LarA.get_external_cycle(
        convert(Lar.Points, model.G'),
        model.T[1],
        model.T[2]
    )
end

function get_external_cycle(
        V::Lar.Points,
        EV::Lar.ChainOp,
        FE::Lar.ChainOp
    )::Union{Int64, Nothing}

    FV = abs.(FE)*EV

    vs = sparsevec(mapslices(sum, abs.(EV), dims=1)').nzind
    minv_x1 = maxv_x1 = minv_x2 = maxv_x2 = pop!(vs)
    for i in vs
        if V[i, 1] > V[maxv_x1, 1]
            maxv_x1 = i
        elseif V[i, 1] < V[minv_x1, 1]
            minv_x1 = i
        end
        if V[i, 2] > V[maxv_x2, 2]
            maxv_x2 = i
        elseif V[i, 2] < V[minv_x2, 2]
            minv_x2 = i
        end
    end
    cells = intersect(
        FV[:, minv_x1].nzind,
        FV[:, maxv_x1].nzind,
        FV[:, minv_x2].nzind,
        FV[:, maxv_x2].nzind
    )
    if length(cells) == 1
        return cells[1]
    else
        for c in cells
            if LAR.face_area(V, EV, FE[c, :]) < 0
                return c
            end
        end
    end
end

"""
    pre_containment_test(
        bboxes::Array{Tuple{Array{Float64,2},Array{Float64,2}},1}
    )::SparseMatrixCSC{Int8,Int64}

Generate the containment graph associated to `bboxes`.

The function evaluate a `SparseArrays{Int8, n, n}` PRE-containmnet graph meaning
it stores `1` when the bbox associated to the `i`-th component contains entirely
the bbox associated to the `j`-th one. (where `n = length(bboxes)`)

See also: [`Lar.Arrangement.componentgraph`](@ref).

---

# Examples
```jldoctest
# Planar Tiles
julia> bboxes = [
		([0.0 0.0], [1.0 1.0])
		([0.0 0.0], [0.4 0.4])
		([0.6 0.0], [1.0 0.4])
		([0.0 0.6], [0.4 1.0])
		([0.6 0.6], [1.0 1.0])
		([2.0 3.0], [2.0 3.0]) # no intersection
	];

julia> LarA.pre_containment_test(bboxes)
6x6 SparseMatrixCSC{Int8,Int64} with 4 stored entries:
  [2, 1]  =  1
  [3, 1]  =  1
  [4, 1]  =  1
  [5, 1]  =  1
```
"""
function pre_containment_test(
        bboxes#::Array{Tuple{Array{Float64,2},Array{Float64,2}},1}              #
    )::SparseMatrixCSC{Int8,Int64}
    n = length(bboxes)
    containment_graph = spzeros(Int8, n, n)

    for i in 1:n
        for j in 1:n
            if i != j && LAR.bbox_contains(bboxes[j], bboxes[i])
                containment_graph[i, j] = 1
            end
        end
    end

    return containment_graph
end

"""
    prune_containment_graph(
            V::Lar.Points,
            EVs::Array{Lar.ChainOp, 1},
            shells::Array{Lar.Chain,1},
            graph::SparseMatrixCSC{Int8, Int64}
        )::SparseMatrixCSC{Int8, Int64}

Prunes the containment `graph` from the non-included faces.

This method prunes the containment graph eliminating the included bouning boxes
that do not corresponds to included faces.

Do note that this method expects to work on biconnectet clusters of faces.

See also: [`Lar.Arrangement.componentgraph`](@ref).
```
"""
function prune_containment_graph(
        V::Lar.Points,
        EVs::Array{Lar.ChainOp, 1},
        shells::Array{Lar.Chain,1},
        graph::SparseMatrixCSC{Int8, Int64}
    )::SparseMatrixCSC{Int8, Int64}

    graph = deepcopy(graph) # copy needed
    n = length(EVs)
    length(shells) == n || throw(ArgumentError("EVs and shells not coherent"))
    size(graph) == (n, n) || throw(ArgumentError("graph dim not coherent"))

    for i in 1 : n
        # Take a point of the i-th shell by taking a vertex of its first edge
        an_edge = shells[i].nzind[1]
        origin_index = EVs[i][an_edge, :].nzind[1]
        origin = V[origin_index, :]

        # Check the containment for every other cell by testing if such a point
        #   is inner to the face defined byt that cell
        for j in 1 : n
            if i != j && graph[i, j] == 1
                shell_edge_indexes = shells[j].nzind
                ev = EVs[j][shell_edge_indexes, :]

                if !LAR.point_in_face(origin, V, ev)
                    graph[i, j] = 0
                end
            end
         end

     end
     return dropzeros(graph)
end

"""
    transitive_reduction!(graph)

Prune the redundancies in containment `graph`.

The ``m``-layer containment `graph` is reduced to the corresponding 1-layer
containment `graph` where only the first level nested components are considered.

See also: [`Lar.Arrangement.componentgraph`](@ref).

---

# Examples
```jldoctest
julia> containment = [
		0 0 0 0 0 0
		0 0 0 0 0 0
		1 0 0 0 0 0
		1 0 0 0 0 0
		1 0 1 0 0 0
		1 0 1 0 1 0
	];

julia> LarA.transitive_reduction!(containment)

julia> containment
6×6 Array{Int64,2}:
 0  0  0  0  0  0
 0  0  0  0  0  0
 1  0  0  0  0  0
 1  0  0  0  0  0
 0  0  1  0  0  0
 0  0  0  0  1  0
```
"""
function transitive_reduction!(graph)
    n = size(graph, 1)
    for j in 1:n
        for i in 1:n
            if graph[i, j] > 0
                for k in 1:n
                    if graph[j, k] > 0
                        graph[i, k] = 0
                    end
                end
            end
        end
    end
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 2 METHODS
#    componentgraph (aka TGW)
#    cell_merging
#-------------------------------------------------------------------------------

"""
    componentgraph(V, copEV, bicon_comps)

Topological Gift Wrapping algorithm on 2D skeletons.

This is the offline part of the TGW algorithm. It takes in input a model and its
biconnected components mapping and evaluates usefull informations:
 1. Number of biconnected components.
 2. Component Graph of the biconnected structure.
 3. The 1-cells structure (UNMODIFIED).                                         # Could be removed?
 4. Association between non-dangling 2-cells and their orientation (∀ component)
 5. Association between 3-cells and 2-cells (with orientation, ∀ component).
 6. Association between 3-cells and their orientation (∀ component).
 7. Shell bounding boxes of the components.

See also: [`Lar.Arrangement.planar_arrangement_2`](@ref).
It uses:
 - [`Lar.Arrangement.get_external_cycle`](@ref)
 - [`Lar.Arrangement.pre_containment_test`](@ref)
 - [`Lar.Arrangement.prune_containment_graph`](@ref)
 - [`Lar.Arrangement.transitive_reduction!`](@ref)
---
# Examples
```jldoctest
# Papillon
julia> V = [0.0 0.0; 0.0 3.0; 2.0 1.5; 4.0 0.0; 4.0 3.0];

julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 0 0 0] #1 -> 1,2
		[0 1 1 0 0] #2 -> 2,3
		[1 0 1 0 0] #3 -> 3,1
		[0 0 1 1 0] #4 -> 3,4
		[0 0 0 1 1] #5 -> 4,5
		[0 0 1 0 1] #6 -> 3,5
	]));

julia> copFE = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 1 0 0 0] #1 -> 1,2,3
		[0 0 0 1 1 1] #2 -> 4,5,6
		[1 1 1 1 1 1] #3 -> 1,2,3,4,5,6 External
	]));

julia> bicon_comps = [[1, 2, 3], [4, 5, 6]];

julia> LarA.componentgraph(V, copEV, bicon_comps)[1]
2

julia> LarA.componentgraph(V, copEV, bicon_comps)[2]
2×2 SparseMatrixCSC{Int8,Int64} with 0 stored entries

julia> LarA.componentgraph(V, copEV, bicon_comps)[4]
2-element Array{SparseMatrixCSC{Int8,Int64},1}:

  [1, 1]  =  -1
  [3, 1]  =  -1
  [1, 2]  =  1
  [2, 2]  =  -1
  [2, 3]  =  1
  [3, 3]  =  1

  [1, 3]  =  -1
  [3, 3]  =  -1
  [1, 4]  =  1
  [2, 4]  =  -1
  [2, 5]  =  1
  [3, 5]  =  1

julia> LarA.componentgraph(V, copEV, bicon_comps)[5]
2-element Array{SparseMatrixCSC{Int8,Int64},1}:

  [1, 1]  =  -1
  [1, 2]  =  -1
  [1, 3]  =  1

  [1, 1]  =  1
  [1, 2]  =  1
  [1, 3]  =  -1

julia> LarA.componentgraph(V, copEV, bicon_comps)[6]
2-element Array{SparseVector{Int8,Int64},1}:
   [1]  =  1
  [2]  =  1
  [3]  =  -1
   [1]  =  -1
  [2]  =  -1
  [3]  =  1

julia> LarA.componentgraph(V, copEV, bicon_comps)[7]
2-element Array{Any,1}:
 ([0.0 0.0], [2.0 3.0])
 ([2.0 0.0], [4.0 3.0])
```
"""
function componentgraph(V, copEV, bicon_comps)
    # arrangement of isolated components
	n = size(bicon_comps, 1)
   	shells = Array{Lar.Chain, 1}(undef, n)
	boundaries = Array{Lar.ChainOp, 1}(undef, n)
	EVs = Array{Lar.ChainOp, 1}(undef, n)
    # for each component
	for p in 1 : n
		ev = copEV[sort(bicon_comps[p]), :]
        # computation of 2-cells
		fe = LarA.minimal_2cycles(V, ev)
        # exterior cycle
		shell_num = LarA.get_external_cycle(V, ev, fe)
        # decompose each fe (co-boundary local to component)
		EVs[p] = ev
		tokeep = setdiff(1:fe.m, shell_num)
		boundaries[p] = fe[tokeep, :]
		shells[p] = fe[shell_num, :]
    end
    # computation of bounding boxes of isolated components
	shell_bboxes = []
	for i in 1 : n
    	vs_indexes = (abs.(EVs[i]')*abs.(shells[i])).nzind
   		push!(shell_bboxes, LAR.bbox(V[vs_indexes, :]))
	end
    # computation and reduction of containment graph
	containment_graph = LarA.pre_containment_test(shell_bboxes)
	containment_graph = LarA.prune_containment_graph(V,EVs,shells,containment_graph)
	LarA.transitive_reduction!(containment_graph)
	return n, containment_graph, V, EVs, boundaries, shells, shell_bboxes
end

"""
    cell_merging(n, containment_graph, V, EVs, boundaries, shells, shell_bboxes)

Cells composing for the Topological Gift Wrapping algorithm.

This is the online part of the TGW algorithm.

See also: [`Lar.Arrangement.planar_arrangement_2`](@ref).
"""
function cell_merging(
        n, containment_graph, V, EVs, boundaries, shells, shell_bboxes
    )
    function bboxes(V::Lar.Points, indexes::Lar.ChainOp)
        boxes = Array{Tuple{Any, Any}}(undef, indexes.n)
        for i in 1:indexes.n
            v_inds = indexes[:, i].nzind
            boxes[i] = LAR.bbox(V[v_inds, :])
        end
        boxes
    end
    # initialization
    sums = Array{Tuple{Int, Int, Int}}(undef, 0)
    # assembling child components with father components
    for father in 1:n
        if sum(containment_graph[:, father]) > 0
            father_bboxes = bboxes(
                V, abs.(EVs[father]')*abs.(boundaries[father]')
            )
            for child in 1:n
                if containment_graph[child, father] > 0
                    child_bbox = shell_bboxes[child]
                    for b in 1:length(father_bboxes)
                        if LAR.bbox_contains(father_bboxes[b], child_bbox)
                            push!(sums, (father, b, child))
                            break
                        end
                    end
                end
            end
        end
    end
    # offset assembly initialization
    EV = vcat(EVs...)
    edgenum = size(EV, 1)
    facenum = sum(map(x->size(x,1), boundaries))
    FE = spzeros(Int8, facenum, edgenum)
    shells2 = spzeros(Int8, length(shells), edgenum)
    r_offsets = [1]
    c_offset = 1
    # submatrices construction
    for i in 1:n
        min_row = r_offsets[end]
        max_row = r_offsets[end] + size(boundaries[i], 1) - 1
        min_col = c_offset
        max_col = c_offset + size(boundaries[i], 2) - 1
        FE[min_row:max_row, min_col:max_col] = boundaries[i]
        shells2[i, min_col:max_col] = shells[i]
        push!(r_offsets, max_row + 1)
        c_offset = max_col + 1
    end
    # offsetting assembly of component submatrices
    for (f, r, c) in sums
        FE[r_offsets[f]+r-1, :] += shells2[c, :]
    end

    return EV, FE
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 1
#-------------------------------------------------------------------------------

"""
	planar_arrangement_1(
        model::Lar.Model,
        [sigma::Lar.Chain = spzeros(Int8, 0)],
        [return_edge_map::Bool = false],
        [multiproc::Bool = false]
    )

First part of arrangement's algorithmic pipeline.

This function computes the pairwise intersection between each edge of a given
2D cellular complex 1-skeleton. The computation is speeded up via the evaluation
of the Spatial Index. See [`Lar.spaceindex`](@ref).

See also: [`Lar.planar_arrangement`](@ref) for the complete pipeline.
It uses:
 - [`Lar.Arrangement.frag_edges_channel`](@ref)
 - [`Lar.Arrangement.frag_edges`](@ref)
 - [`Lar.Arrangement.merge_vertices!`](@ref)
---
# Examples
```jldoctest
julia> model = Lar.Model([
    0.0 0.5 0.0 0.5 0.3 1.0 0.3 1.0;
    0.0 0.0 1.0 1.0 0.5 0.5 1.0 1.0
]);
julia> Lar.addModelCells!(model, 1, [
    [1, 2], [3, 4], [1, 3], [2, 4], [5, 6], [7, 8], [5, 7], [6, 8]
]);
julia> (model, edge_map) =
    LarA.planar_arrangement_1(model, spzeros(Int8, 0), true);
julia> model.G
2×9 Array{Float64,2}:
 0.0  0.5  0.0  0.5  0.3  0.5  0.3  1.0  1.0
 0.0  0.0  1.0  1.0  1.0  0.5  0.5  0.5  1.0

julia> Matrix(model.T[1])
11×9 Array{Int8,2}:
 1  1  0  0  0  0  0  0  0
 0  0  1  1  0  0  0  0  0
 0  0  0  1  1  0  0  0  0
 1  0  1  0  0  0  0  0  0
 0  1  0  0  0  1  0  0  0
 0  0  0  0  1  1  0  0  0
 0  0  0  0  0  1  1  0  0
 0  0  0  0  0  1  0  1  0
 0  0  0  0  1  0  0  0  1
 0  0  0  1  0  0  1  0  0
 0  0  0  0  0  0  0  1  1

julia> edge_map
8-element Array{Array{Int64,1},1}:
 [1]
 [2, 3]
 [4]
 [5, 6]
 [7, 8]
 [3, 9]
 [10]
 [11]

```
"""
function planar_arrangement_1(
        model::Lar.Model,
		sigma::Lar.Chain = spzeros(Int8, 0),
		return_edge_map::Bool = false,
		multiproc::Bool = false
    )::Union{
        Lar.Model,
        Tuple{Lar.Model, Array{Array{Int64,1},1}},
        Tuple{Lar.Model, Array{Int,1}},
        Tuple{Lar.Model, Array{Int,1}, Array{Array{Int64,1},1}}
    }
	# data structures initialization
    # V = convert(Lar.Points, model.G')
    # copEV = model.T[1]
	edgenum = size(model, 1, 1)
	edge_map = Array{Array{Int, 1}, 1}(undef,edgenum)
	rV = Lar.Points(zeros(0, 2))
	rEV = SparseArrays.spzeros(Int8, 0, 0)
	finalcells_num = 0

	# spaceindex computation
	bigPI = LarA.spaceIndex(model, 1)

    # multiprocessing of edge fragmentation
    if (multiproc == true)
        in_chan = Distributed.RemoteChannel(()->Channel{Int64}(0))
        out_chan = Distributed.RemoteChannel(()->Channel{Tuple}(0))
        ordered_dict = DataStructures.SortedDict{Int64,Tuple}()
        @async begin
            for i in 1:edgenum
                put!(in_chan,i)
            end
            for p in Distributed.workers()
                put!(in_chan,-1)
            end
        end
        for p in Distributed.workers()
            @async Base.remote_do(
                frag_edge_channel, p, in_chan, out_chan, V, model.T[1], bigPI
            )
        end
        for i in 1:edgenum
            frag_done_job = take!(out_chan)
            ordered_dict[frag_done_job[1]] = frag_done_job[2]
        end
        for (dkey, dval) in ordered_dict
            i = dkey
            v, ev = dval
            newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))
            edge_map[i] = newedges_nums
            finalcells_num += size(ev, 1)
            rV, rEV = LAR.skel_merge(rV, rEV, v, ev)
        end
    else
        # sequential (iterative) processing of edge fragmentation
        for i in 1:edgenum
            v, ev = LarA.frag_edge(model, i, bigPI)
            newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))
            edge_map[i] = newedges_nums
            finalcells_num += size(ev, 1)
            rV = convert(Lar.Points, rV)
            rV, rEV = LAR.skel_merge(rV, rEV, v, ev)
        end
    end

    # merging of close vertices and edges (2D congruence)
    model = Lar.Model(convert(Lar.Points, rV'))
    Lar.addModelCells!(model, 1, rEV)
    LarA.merge_vertices!(model, edge_map)

    if isempty(sigma)
        if return_edge_map
            return model, edge_map
        else
            return model
        end
    end

    # Sigma update by mapping each of its edges into newly born edges
    new_sigma = Array{Int,1}()
    map(i->new_sigma=union(new_sigma, edge_map[i]), sigma.nzind)
    if return_edge_map
        return model, new_sigma, edge_map
    else
        return model, new_sigma
    end
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE - PART 2
#-------------------------------------------------------------------------------

"""
	planar_arrangement_2(
        model::Lar.Model,
        bicon_comps::{Array{Array,1},1}
    )::Lar.Model

Second part of arrangement's algorithmic pipeline.

This function is the complete Topological Gift Wrapping (TGW) algorithm that
is firstly locally used in order to decompose the 2-cells and then globally
to generate the 3-cells of the arrangement of the ambient space ``E^3``.

During this process each dangling 1-cell is removed.
Do note that the isolated 2-cells are not removed by this procedure.

See also: [`Lar.planar_arrangement`](@ref) for the complete pipeline.
It uses:
 - [`Lar.Arrangement.biconnected_components`](@ref)
 - [`Lar.Arrangement.componentgraph`](@ref)
 - [`Lar.Arrangement.cell_merging`](@ref)
---

# Examples
```jldoctest
# Triforce
julia> model = Lar.Model([
    0.0 0 2 0 0 1 4 4 3 2 4 4
    0.0 2 0 3 4 4 2 0 0 4 4 2
]);

julia> Lar.addModelCells!(model, 1, [
    [ 1,  2], [ 2,  3], [ 3,  1],
    [ 4,  5], [ 5,  6], [ 6,  7], [ 7,  8], [ 8,  9], [ 9, 4],
    [10, 11], [11, 12], [12, 10]
])

julia> bicon_comps = LarA.biconnected_components(model.T[1])
3-element Array{Array{Int64,1},1}:
 [9, 8, 7, 6, 5, 4]
 [3, 2, 1]
 [12, 11, 10]

julia> model = Lar.Arrangement.planar_arrangement_2(model, bicon_comps);

julia> Matrix(model.T[2])
3×12 Array{Int8,2}:
 -1  -1  -1  -1  -1  1   0   0  0   0   0  0
  0   0   0   0   0  0  -1  -1  1   0   0  0
  0   0   0   0   0  0   0   0  0  -1  -1  1
```
"""
function planar_arrangement_2(
        model::Lar.Model,
        bicon_comps::Array{Array{Int64,1},1}
    )::Lar.Model

    #==
    test to remove
    @assert isequal(
        bicon_comps,
        LarA.biconnected_components(model.T[1])
    )
    ==#

	# component graph
	n, containment_graph, V, EVs, boundaries, shells, shell_bboxes =
        LarA.componentgraph(
            convert(Lar.Points, model.G'), model.T[1], bicon_comps
        )

	copEV, FE = LarA.cell_merging(
	   	n, containment_graph, V, EVs, boundaries, shells, shell_bboxes
    )

    model = Lar.Model(convert(Lar.Points, V'))
    Lar.addModelCells!(model, 1, copEV)
    Lar.addModelCells!(model, 2, FE)
    Lar.modelPurge!(model)
	return model
end

#-------------------------------------------------------------------------------
#   PLANAR ARRANGEMENT PIPELINE
#-------------------------------------------------------------------------------

"""
    planar_arrangement(
        model::Lar.Model,
        [sigma::Lar.Chain],
        [return_edge_map::Bool],
        [multiproc::Bool]
    )::Union{Lar.Model, Tuple{Lar.model, Array{Array{Int, 1}, 1}}}

Compute the arrangement on the given cellular complex 1-skeleton in 2D.

A cellular complex is arranged when the intersection of every possible pair
of cell of the complex is empty and the union of all the cells is the
whole Euclidean space.
The methods default input/output is a `Lar.Model` with at least Geometry and
1-cell Topoly.
Mutiprocessing is by default disabled.

It uses:
 - [`Lar.Arrangement.planar_arrangement_1`](@ref)
 - [`Lar.Arrangement.cleandecomposition`](@ref)
 - [`Lar.Arrangement.biconnected_components`](@ref)
 - [`Lar.Arrangement.remove_mapping_dangling_edges`](@ref)
 - [`Lar.Arrangement.planar_arrangement_2`](@ref)

## Additional arguments:
 - `sigma::Chain`: if specified, the method return only cells contained inside
        this cell. _Default_: empty cell.
 - `return_edge_map::Bool`: builds an `edge_map` which maps the input edges
        to the output's. _Defaults_: `false`.
 - `multiproc::Bool`: Runs the method in parallel mode. _Defaults_: `false`.
"""
function planar_arrangement(
        model::Lar.Model,
        sigma::Lar.Chain = spzeros(Int8, 0),                                    ## TEST
        return_edge_map::Bool = false,
        multiproc::Bool = false
    )::Union{Lar.Model, Tuple{Lar.Model, Array{Array{Int, 1}, 1}}}

    # Check dimensional condition
    @assert length(model) == 2

    # Chek multiprocessing
    if multiproc && Distributed.nprocs() <= 1
        multiproc = false
        println("Setting Multiproc to False due to insufficient processors.\n")
    end

    #planar_arrangement_1
    if isempty(sigma)
        model, edge_map =
            LarA.planar_arrangement_1(model, sigma, true, multiproc)
    else
	    model, sigma, edge_map =
            LarA.planar_arrangement_1(model, sigma, true, multiproc)
    end

    print("PLANAR_ARRANGEMENT_1 COMPLETE\n")

    # cleandecomposition
	if !isempty(sigma)
		model, edge_map = LarA.cleandecomposition(model, sigma, edge_map)
        print("DECOMPOSITION CLEANED\n")
	end

    # generation of biconnected components
    bicon_comps = LarA.biconnected_components(model.T[1])

    print("BICONNECTED COMPONENTS FOUND\n")

	if isempty(bicon_comps)
    	println("No biconnected components found.\n")
    	if (return_edge_map)
    	    return (model(), nothing)
    	else
    	    return model()
    	end
	end

    # Removing Dangling edges from edge map
    if return_edge_map
        @assert size(model, 1, 1) == max((edge_map...)...)                           # CTR
        edge_map = LarA.remove_mapping_dangling_edges(edge_map, bicon_comps)
    end

    print("DANGLING EDGES REMOVED\n")

    # Planar_arrangement_2
	model = LarA.planar_arrangement_2(model, bicon_comps)

    print("PLANAR ARRANGEMENT 2 COMPLETE\n")

    if size(model, 0, 2) - size(model, 1, 1) + size(model, 2, 1) != 1
        @show size(model, 0, 2) - size(model, 1, 1) + size(model, 2, 1)
    end

	if return_edge_map
	     return model, edge_map
	else
	     return model
	end
end


text:
struct NodeError{T <: NodeValue} <: NodeValue
    o::T
end
NodeError(value::NapiValue) = NodeError(node_value(value))

function Base.show(io::IO, v::NodeError)
    print(io, "NodeError(")
    print(io, getfield(v, :o))
    print(io, ")")
end
function Base.showerror(io::IO, v::NodeError)
    print(io, "NodeError: ")
    print(io, v.message)
end
function Base.showerror(io::IO, v::NodeError{NodeValueTemp})
    print(io, "NodeError: ")
    print(io, value(v))
end
napi_value(node_error::NodeError) = napi_value(getfield(node_error, :o))

node_throw(err::NapiValue) = @napi_call napi_throw(err::NapiValue)
node_throw(err::Exception) = @with_scope node_throw(napi_value(err))
function node_throw(x)
    err = ErrorException(sprint(showerror, x))
    node_throw(err)
end
napi_value(err::Exception) = @napi_call napi_create_error(
    C_NULL::NapiValue, sprint(showerror, err)::NapiValue
)::NapiValue


text:
# ---
# title: 376. Wiggle Subsequence
# id: problem376
# author: Indigo
# date: 2020-12-12
# difficulty: Medium
# categories: Dynamic Programming, Greedy
# link: <https://leetcode.com/problems/wiggle-subsequence/description/>
# hidden: true
# ---
#
# A sequence of numbers is called a **wiggle sequence** if the differences
# between successive numbers strictly alternate between positive and negative.
# The first difference (if one exists) may be either positive or negative. A
# sequence with fewer than two elements is trivially a wiggle sequence.
#
# For example, `[1,7,4,9,2,5]` is a wiggle sequence because the differences
# `(6,-3,5,-7,3)` are alternately positive and negative. In contrast,
# `[1,4,7,2,5]` and `[1,7,4,5,5]` are not wiggle sequences, the first because
# its first two differences are positive and the second because its last
# difference is zero.
#
# Given a sequence of integers, return the length of the longest subsequence
# that is a wiggle sequence. A subsequence is obtained by deleting some number
# of elements (eventually, also zero) from the original sequence, leaving the
# remaining elements in their original order.
#
# **Example 1:**
#
#
#
#     Input: [1,7,4,9,2,5]
#     Output: 6
#     Explanation: The entire sequence is a wiggle sequence.
#
# **Example 2:**
#
#
#
#     Input: [1,17,5,10,13,15,10,5,16,8]
#     Output: 7
#     Explanation: There are several subsequences that achieve this length. One is [1,17,10,13,10,16,8].
#
# **Example 3:**
#
#
#
#     Input: [1,2,3,4,5,6,7,8,9]
#     Output: 2
#
# **Follow up:**
# Can you do it in O( _n_ ) time?
#
#
## @lc code=start
using LeetCode
function wiggle_max_length(nums::Vector{Int})::Int
    len = length(nums)
    if len < 2
        return len
    end
    pre_diff = nums[2] - nums[1]
    res = (pre_diff == 0 ? 1 : 2)
    for i in 3:len
        diff = nums[i] - nums[i - 1]
        if diff > 0 && pre_diff <= 0 || diff < 0 && pre_diff > 0
            res += 1
            pre_diff = diff
        end
    end
    return res
end
## @lc code=end