text:
@testset "iban                [one]" begin
        result = iban(
                CountryCode="AL",
                BankCode="212",
                BranchCode="1100",
                AccountNumber="0000000235698741",
                NationalCheckDigit="9"
        )
        @test result["CountryCode"] == "AL"
        @test result["BankCode"] == "212"
        @test result["BranchCode"] == "1100"
        @test result["AccountNumber"] == "0000000235698741"
        @test result["NationalCheckDigit"] == "9"
        @test result["value"] == "AL47212110090000000235698741"
end

@testset "iban                [all]" begin
        @test iban(
                CountryCode="AL",
                BankCode="212",
                BranchCode="1100",
                AccountNumber="0000000235698741",
                NationalCheckDigit="9",
        )["value"] == "AL47212110090000000235698741"
        @test iban(
                CountryCode="AD",
                BankCode="0001",
                BranchCode="2030",
                AccountNumber="200359100100",
        )["value"] == "AD1200012030200359100100"
        @test iban(
                CountryCode="AT",
                BankCode="19043",
                AccountNumber="00234573201",
        )["value"] == "AT611904300234573201"
        @test iban(
                CountryCode="AZ",
                BankCode="NABZ",
                AccountNumber="00000000137010001944",
        )["value"] == "AZ21NABZ00000000137010001944"
        @test iban(
                CountryCode="BH",
                BankCode="SCBL",
                AccountNumber="BHD18903608801",
        )["value"] == "BH72SCBLBHD18903608801"
        @test iban(
                CountryCode="BE",
                BankCode="539",
                AccountNumber="0075470",
                NationalCheckDigit="34",
        )["value"] == "BE68539007547034"
        @test iban(
                CountryCode="BA",
                BankCode="129",
                BranchCode="007",
                AccountNumber="94010284",
                NationalCheckDigit="94",
        )["value"] == "BA391290079401028494"
        @test iban(
                CountryCode="BR",
                BankCode="00360305",
                BranchCode="00001",
                AccountNumber="0009795493",
                AccountType="P",
                OwnerAccountType="1",
        )["value"] == "BR9700360305000010009795493P1"
        @test iban(
                CountryCode="BG",
                BankCode="BNBG",
                BranchCode="9661",
                AccountNumber="20345678",
                AccountType="10",
        )["value"] == "BG80BNBG96611020345678"
        @test iban(
                CountryCode="BY",
                BankCode="NBRB",
                BranchCode="3600",
                AccountNumber="900000002Z00AB00",
        )["value"] == "BY13NBRB3600900000002Z00AB00"
        @test iban(
                CountryCode="CR",
                BankCode="0152",
                AccountNumber="02001026284066",
        )["value"] == "CR05015202001026284066"
        @test iban(
                CountryCode="HR",
                BankCode="1001005",
                AccountNumber="1863000160",
        )["value"] == "HR1210010051863000160"
        @test iban(
                CountryCode="CY",
                BankCode="002",
                BranchCode="00128",
                AccountNumber="0000001200527600",
        )["value"] == "CY17002001280000001200527600"
        @test iban(
                CountryCode="CZ",
                BankCode="0800",
                AccountNumber="0000192000145399",
        )["value"] == "CZ6508000000192000145399"
        @test iban(
                CountryCode="DK",
                BankCode="0040",
                AccountNumber="0440116243",
        )["value"] == "DK5000400440116243"
        @test iban(
                CountryCode="DO",
                BankCode="BAGR",
                AccountNumber="00000001212453611324",
        )["value"] == "DO28BAGR00000001212453611324"
        @test iban(
                CountryCode="EE",
                BankCode="22",
                BranchCode="00",
                AccountNumber="22102014568",
                NationalCheckDigit="5",
        )["value"] == "EE382200221020145685"
        @test iban(
                CountryCode="FI",
                BankCode="123456",
                AccountNumber="0000078",
                NationalCheckDigit="5",
        )["value"] == "FI2112345600000785"
        @test iban(
                CountryCode="FR",
                BankCode="20041",
                BranchCode="01005",
                AccountNumber="0500013M026",
                NationalCheckDigit="06",
        )["value"] == "FR1420041010050500013M02606"
        @test iban(
                CountryCode="GE",
                BankCode="NB",
                AccountNumber="0000000101904917",
        )["value"] == "GE29NB0000000101904917"
        @test iban(
                CountryCode="DE",
                BankCode="37040044",
                AccountNumber="0532013000",
        )["value"] == "DE89370400440532013000"
        @test iban(
                CountryCode="GI",
                BankCode="NWBK",
                AccountNumber="000000007099453",
        )["value"] == "GI75NWBK000000007099453"
        @test iban(
                CountryCode="GR",
                BankCode="011",
                BranchCode="0125",
                AccountNumber="0000000012300695",
        )["value"] == "GR1601101250000000012300695"
        @test iban(
                CountryCode="GT",
                BankCode="TRAJ",
                AccountNumber="01020000001210029690",
        )["value"] == "GT82TRAJ01020000001210029690"
        @test iban(
                CountryCode="HU",
                BankCode="117",
                BranchCode="7301",
                AccountNumber="6111110180000000",
                NationalCheckDigit="0",
        )["value"] == "HU42117730161111101800000000"
        @test iban(
                CountryCode="IS",
                BankCode="0159",
                BranchCode="26",
                AccountNumber="007654",
                IdentificationNumber="5510730339",
        )["value"] == "IS140159260076545510730339"
        @test iban(
                CountryCode="IE",
                BankCode="AIBK",
                BranchCode="931152",
                AccountNumber="12345678",
        )["value"] == "IE29AIBK93115212345678"
        @test iban(
                CountryCode="IL",
                BankCode="010",
                BranchCode="800",
                AccountNumber="0000099999999",
        )["value"] == "IL620108000000099999999"
        @test iban(
                CountryCode="IT",
                BankCode="05428",
                BranchCode="11101",
                NationalCheckDigit="X",
                AccountNumber="000000123456",
        )["value"] == "IT60X0542811101000000123456"
        @test iban(
                CountryCode="JO",
                BankCode="CBJO",
                BranchCode="0010",
                AccountNumber="000000000131000302",
        )["value"] == "JO94CBJO0010000000000131000302"
        @test iban(
                CountryCode="KZ",
                BankCode="125",
                AccountNumber="KZT5004100100",
        )["value"] == "KZ86125KZT5004100100"
        @test iban(
                CountryCode="KW",
                BankCode="CBKU",
                AccountNumber="0000000000001234560101",
        )["value"] == "KW81CBKU0000000000001234560101"
        @test iban(
                CountryCode="LC",
                BankCode="HEMM",
                AccountNumber="000100010012001200023015",
        )["value"] == "LC55HEMM000100010012001200023015"
        @test iban(
                CountryCode="LV",
                BankCode="BANK",
                AccountNumber="0000435195001",
        )["value"] == "LV80BANK0000435195001"
        @test iban(
                CountryCode="LB",
                BankCode="0999",
                AccountNumber="00000001001901229114",
        )["value"] == "LB62099900000001001901229114"
        @test iban(
                CountryCode="LI",
                BankCode="08810",
                AccountNumber="0002324013AA",
        )["value"] == "LI21088100002324013AA"
        @test iban(
                CountryCode="LT",
                BankCode="10000",
                AccountNumber="11101001000",
        )["value"] == "LT121000011101001000"
        @test iban(
                CountryCode="LU",
                BankCode="001",
                AccountNumber="9400644750000",
        )["value"] == "LU280019400644750000"
        @test iban(
                CountryCode="MK",
                BankCode="250",
                AccountNumber="1200000589",
                NationalCheckDigit="84",
        )["value"] == "MK07250120000058984"
        @test iban(
                CountryCode="MT",
                BankCode="MALT",
                BranchCode="01100",
                AccountNumber="0012345MTLCAST001S",
        )["value"] == "MT84MALT011000012345MTLCAST001S"
        @test iban(
                CountryCode="MR",
                BankCode="00020",
                BranchCode="00101",
                AccountNumber="00001234567",
                NationalCheckDigit="53",
        )["value"] == "MR1300020001010000123456753"
        @test iban(
                CountryCode="MU",
                BankCode="BOMM01",
                BranchCode="01",
                AccountNumber="101030300200000MUR",
        )["value"] == "MU17BOMM0101101030300200000MUR"
        @test iban(
                CountryCode="MD",
                BankCode="AG",
                AccountNumber="000225100013104168",
        )["value"] == "MD24AG000225100013104168"
        @test iban(
                CountryCode="MC",
                BankCode="11222",
                BranchCode="00001",
                AccountNumber="01234567890",
                NationalCheckDigit="30",
        )["value"] == "MC5811222000010123456789030"
        @test iban(
                CountryCode="ME",
                BankCode="505",
                AccountNumber="0000123456789",
                NationalCheckDigit="51",
        )["value"] == "ME25505000012345678951"
        @test iban(
                CountryCode="NL",
                BankCode="ABNA",
                AccountNumber="0417164300",
        )["value"] == "NL91ABNA0417164300"
        @test iban(
                CountryCode="NO",
                BankCode="8601",
                AccountNumber="111794",
                NationalCheckDigit="7",
        )["value"] == "NO9386011117947"
        @test iban(
                CountryCode="PK",
                BankCode="SCBL",
                AccountNumber="0000001123456702",
        )["value"] == "PK36SCBL0000001123456702"
        @test iban(
                CountryCode="PS",
                BankCode="PALS",
                AccountNumber="000000000400123456702",
        )["value"] == "PS92PALS000000000400123456702"
        @test iban(
                CountryCode="PL",
                BankCode="109",
                BranchCode="0101",
                AccountNumber="0000071219812874",
                NationalCheckDigit="4",
        )["value"] == "PL61109010140000071219812874"
        @test iban(
                CountryCode="PT",
                BankCode="0002",
                BranchCode="0123",
                AccountNumber="12345678901",
                NationalCheckDigit="54",
        )["value"] == "PT50000201231234567890154"
        @test iban(
                CountryCode="RO",
                BankCode="AAAA",
                AccountNumber="1B31007593840000",
        )["value"] == "RO49AAAA1B31007593840000"
        @test iban(
                CountryCode="QA",
                BankCode="DOHB",
                AccountNumber="00001234567890ABCDEFG",
        )["value"] == "QA58DOHB00001234567890ABCDEFG"
        @test iban(
                CountryCode="SC",
                BankCode="SSCB",
                BranchCode="1101",
                AccountNumber="0000000000001497",
                AccountType="USD",
        )["value"] == "SC18SSCB11010000000000001497USD"
        @test iban(
                CountryCode="SM",
                BankCode="03225",
                BranchCode="09800",
                AccountNumber="000000270100",
                NationalCheckDigit="U",
        )["value"] == "SM86U0322509800000000270100"
        @test iban(
                CountryCode="ST",
                BankCode="0001",
                BranchCode="0001",
                AccountNumber="0051845310112",
        )["value"] == "ST68000100010051845310112"
        @test iban(
                CountryCode="SA",
                BankCode="80",
                AccountNumber="000000608010167519",
        )["value"] == "SA0380000000608010167519"
        @test iban(
                CountryCode="RS",
                BankCode="260",
                BranchCode="26",
                AccountNumber="0056010016113",
                NationalCheckDigit="79",
        )["value"] == "RS35260005601001611379"
        @test iban(
                CountryCode="SK",
                BankCode="1200",
                AccountNumber="0000198742637541",
        )["value"] == "SK3112000000198742637541"
        @test iban(
                CountryCode="SV",
                BankCode="CENR",
                AccountNumber="00000000000000700025",
        )["value"] == "SV62CENR00000000000000700025"
        @test iban(
                CountryCode="SI",
                BankCode="26",
                BranchCode="330",
                AccountNumber="00120390",
                NationalCheckDigit="86",
        )["value"] == "SI56263300012039086"
        @test iban(
                CountryCode="ES",
                BankCode="2100",
                BranchCode="0418",
                AccountNumber="0200051332",
                NationalCheckDigit="45",
        )["value"] == "ES9121000418450200051332"
        @test iban(
                CountryCode="SE",
                BankCode="500",
                AccountNumber="00000058398257466",
        )["value"] == "SE4550000000058398257466"
        @test iban(
                CountryCode="CH",
                BankCode="00762",
                AccountNumber="011623852957",
        )["value"] == "CH9300762011623852957"
        @test iban(
                CountryCode="TN",
                BankCode="10",
                BranchCode="006",
                AccountNumber="035183598478831",
        )["value"] == "TN5910006035183598478831"
        @test iban(
                CountryCode="TR",
                BankCode="00061",
                AccountNumber="0519786457841326",
                NationalCheckDigit="0",
        )["value"] == "TR330006100519786457841326"
        @test iban(
                CountryCode="AE",
                BankCode="033",
                AccountNumber="1234567890123456",
        )["value"] == "AE070331234567890123456"
        @test iban(
                CountryCode="GB",
                BankCode="NWBK",
                BranchCode="601613",
                AccountNumber="31926819",
        )["value"] == "GB29NWBK60161331926819"
        @test iban(
                CountryCode="VG",
                BankCode="VPVG",
                AccountNumber="0000012345678901",
        )["value"] == "VG96VPVG0000012345678901"
        @test iban(
                CountryCode="TL",
                BankCode="008",
                AccountNumber="00123456789101",
                NationalCheckDigit="57",
        )["value"] == "TL380080012345678910157"
        @test iban(
                CountryCode="XK",
                BankCode="10",
                BranchCode="00",
                AccountNumber="0000000000",
                NationalCheckDigit="53",
        )["value"] == "XK051000000000000053"
        @test iban(
                CountryCode="IR",
                BankCode="017",
                AccountNumber="0000000000123456789",
        )["value"] == "IR200170000000000123456789"
        @test iban(
                CountryCode="FO",
                BankCode="5432",
                AccountNumber="038889994",
                NationalCheckDigit="4",
        )["value"] == "FO9754320388899944"
        @test iban(
                CountryCode="GL",
                BankCode="6471",
                AccountNumber="0001000206",
        )["value"] == "GL8964710001000206"
        @test iban(
                CountryCode="UA",
                BankCode="354347",
                AccountNumber="0006762462054925026",
        )["value"] == "UA573543470006762462054925026"

end


text:
import SatSolver: is_clause_satisfied, is_clause_unit, PartialAssignment,
  determine_variable_count, is_clause_falsified, find_unit_literal

@testset "Formula" begin
  @testset "Determine variable count from a vector of clauses" begin
    function test_case(clauses::Vector{Clause}, ref)
      @test determine_variable_count(clauses) == ref
    end

    test_case(Clause[Clause(-1)], 1)
    test_case(Clause[Clause(-2)], 2)
    test_case(Clause[Clause(2)], 2)
    test_case(Clause[Clause(2), Clause(40)], 40)
  end

  @testset "Is clause satisfied" begin
    cl = Clause(1, 4, -6)
    function test_case(pa::PartialAssignment, ref::Bool)
      @test is_clause_satisfied(cl, pa) == ref
    end

    pa = PartialAssignment(6)
    test_case(pa, false)
    pa[1] = 1
    test_case(pa, true)
    pa[1] = 0
    pa[4] = -1
    test_case(pa, false)
    pa[1] = -1
    test_case(pa, false)
    pa[6] = 1
    test_case(pa, false)
    pa[6] = -1
    test_case(pa, true)

    cl = Clause(-1, 2)
    test_case(PartialAssignment([-1, 1]), true)
    test_case(PartialAssignment([-1, -1]), true)
    test_case(PartialAssignment([-1, 0]), true)
    test_case(PartialAssignment([0, -1]), false)
    test_case(PartialAssignment([0, 0]), false)
    test_case(PartialAssignment([0, 1]), true)
    test_case(PartialAssignment([1, -1]), false)
    test_case(PartialAssignment([1, 0]), false)
    test_case(PartialAssignment([1, 1]), true)
  end

  @testset "Is clause falsified" begin
    cl = Clause(1, -2)
    function test_case(pa, ref)
      @test ref == is_clause_falsified(cl, PartialAssignment(pa))
    end

    test_case([-1, 1], true)
    test_case([-1, -1], false)
    test_case([-1, 0], false)
    test_case([0, -1], false)
    test_case([0, 0], false)
    test_case([0, 1], false)
    test_case([1, -1], false)
    test_case([1, 0], false)
    test_case([1, 1], false)
  end

  @testset "Is clause unit" begin
    function test_case(cl::Clause, pa::PartialAssignment, ref::Bool)
      @test ref == is_clause_unit(cl, pa)
    end

    cl = Clause(1, 2, -3)
    pa = PartialAssignment(3)
    test_case(cl, pa, false)
    pa[1] = 1
    test_case(cl, pa, false)
    pa[1] = -1
    pa[3] = 1
    test_case(cl, pa, true)
  end

  @testset "Find unit literal in a clause" begin
    cl = Clause(-1, 2)
    function test_case(pa, ref)
      @test ref == find_unit_literal(cl, PartialAssignment(pa))
    end

    test_case([-1, 1], 0)
    test_case([-1, -1], 0)
    test_case([-1, 0], 0)
    test_case([0, -1], -1)
    test_case([0, 0], 0)
    test_case([0, 1], 0)
    test_case([1, -1], 0)
    test_case([1, 0], 2)
    test_case([1, 1], 0)
  end
end

ext: jl


text:
# span.jl

@doc """
### findSpan

Determines the knot span index, i

**Inputs**

*  `u` : coordinate value (u,v,w is the coordinate space in 3D)
*  `U` : Knot vector
*  `k` : B-spline order
*  `nctl` : Number of control points

**Outputs**

*  `span` : Knot span index

SOURCE: The NURBS book 2nd Edition, Algorithm A2.1

"""->

function findSpan(u::Tffd, U::AbstractArray{Tffd,1}, k::Int, nctl::Int) where Tffd

  if u >= U[nctl+1]
    return nctl  # Special case when u = last term of knot vector
  elseif u < U[k]
    return k  # When u lies at the starting point of the curve
  end

  low = k-1
  high = nctl+1
  span = div(low + high, 2)
  # Do a binary search
  while u < U[span] || u >= U[span+1]
    if u < U[span]
      high = span
    else
      low = span
    end  # End if
    span = div(low+high, 2)
  end  # End while

  return span
end  # End function findSpan


text:
####################
# LU Factorization #
####################
immutable LU{T,S<:AbstractMatrix{T}} <: Factorization{T}
    factors::S
    ipiv::Vector{BlasInt}
    info::BlasInt
end

# StridedMatrix
function lufact!{T<:BlasFloat}(A::StridedMatrix{T}; pivot = true)
    !pivot && return generic_lufact!(A, pivot=pivot)
    lpt = LAPACK.getrf!(A)
    return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
end
lufact!(A::StridedMatrix; pivot = true) = generic_lufact!(A, pivot=pivot)
function generic_lufact!{T}(A::StridedMatrix{T}; pivot = true)
    m, n = size(A)
    minmn = min(m,n)
    info = 0
    ipiv = Array(BlasInt, minmn)
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if pivot
                amax = real(zero(T))
                for i = k:m
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            ipiv[k] = kp
            if A[kp,k] != 0
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                for i = k+1:m
                    A[i,k] *= Akkinv
                end
            elseif info == 0
                info = k
            end
            # Update the rest
            for j = k+1:n
                for i = k+1:m
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    LU{T,typeof(A)}(A, ipiv, convert(BlasInt, info))
end
lufact{T<:BlasFloat}(A::AbstractMatrix{T}; pivot = true) = lufact!(copy(A), pivot=pivot)
lufact{T}(A::AbstractMatrix{T}; pivot = true) = (S = typeof(one(T)/one(T)); S != T ? lufact!(convert(AbstractMatrix{S}, A), pivot=pivot) : lufact!(copy(A), pivot=pivot))
lufact(x::Number) = LU(fill(x, 1, 1), BlasInt[1], x == 0 ? one(BlasInt) : zero(BlasInt))
lufact(F::LU) = F

lu(x::Number) = (one(x), x, 1)
function lu(A::AbstractMatrix; pivot = true)
    F = lufact(A, pivot = pivot)
    F[:L], F[:U], F[:p]
end

function convert{T}(::Type{LU{T}}, F::LU)
    M = convert(AbstractMatrix{T}, F.factors)
    LU{T,typeof(M)}(M, F.ipiv, F.info)
end
convert{T,S}(::Type{LU{T,S}}, F::LU) = LU{T,S}(convert(S, F.factors), F.ipiv, F.info)
convert{T}(::Type{Factorization{T}}, F::LU) = convert(LU{T}, F)

size(A::LU) = size(A.factors)
size(A::LU,n) = size(A.factors,n)

function ipiv2perm{T}(v::AbstractVector{T}, maxi::Integer)
    p = T[1:maxi]
    @inbounds for i in 1:length(v)
        p[i], p[v[i]] = p[v[i]], p[i]
    end
    return p
end

function getindex{T,S<:StridedMatrix}(A::LU{T,S}, d::Symbol)
    m, n = size(A)
    if d == :L
        L = tril!(A.factors[1:m, 1:min(m,n)])
        for i = 1:min(m,n); L[i,i] = one(T); end
        return L
    end
    d == :U && return triu!(A.factors[1:min(m,n), 1:n])
    d == :p && return ipiv2perm(A.ipiv, m)
    if d == :P
        p = A[:p]
        P = zeros(T, m, m)
        for i in 1:m
            P[i,p[i]] = one(T)
        end
        return P
    end
    throw(KeyError(d))
end

A_ldiv_B!{T<:BlasFloat, S<:StridedMatrix}(A::LU{T, S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('N', A.factors, A.ipiv, B) A.info
A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, b::StridedVector) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), b[ipiv2perm(A.ipiv, length(b))]))
A_ldiv_B!{T,S<:StridedMatrix}(A::LU{T,S}, B::StridedMatrix) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), B[ipiv2perm(A.ipiv, size(B, 1)),:]))
At_ldiv_B{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, copy(B)) A.info
Ac_ldiv_B{T<:BlasComplex,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, copy(B)) A.info
At_ldiv_Bt{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, transpose(B)) A.info
Ac_ldiv_Bc{T<:BlasComplex,S<:StridedMatrix}(A::LU{T,S}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, ctranspose(B)) A.info

function det{T,S}(A::LU{T,S})
    n = chksquare(A)
    A.info > 0 && return zero(typeof(A.factors[1]))
    return prod(diag(A.factors)) * (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
end

function logdet2{T<:Real,S}(A::LU{T,S})  # return log(abs(det)) and sign(det)
    n = chksquare(A)
    dg = diag(A.factors)
    s = (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T)) * prod(sign(dg))
    sum(log(abs(dg))), s
end

function logdet{T<:Real,S}(A::LU{T,S})
    d,s = logdet2(A)
    s>=0 || error("DomainError: determinant is negative")
    d
end

function logdet{T<:Complex,S}(A::LU{T,S})
    n = chksquare(A)
    s = sum(log(diag(A.factors))) + (bool(sum(A.ipiv .!= 1:n) % 2) ? complex(0,pi) : 0)
    r, a = reim(s)
    a = pi-mod(pi-a,2pi) #Take principal branch with argument (-pi,pi]
    complex(r,a)
end

inv{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}) = @assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info

cond{T<:BlasFloat,S<:StridedMatrix}(A::LU{T,S}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm((A[:L]*A[:U])[A[:p],:], p)))
cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p)

# Tridiagonal

# See dgttrf.f
function lufact!{T}(A::Tridiagonal{T}; pivot = true)
    n = size(A, 1)
    info = 0
    ipiv = Array(BlasInt, n)
    dl = A.dl
    d = A.d
    du = A.du
    du2 = A.du2

    @inbounds begin
        for i = 1:n
            ipiv[i] = i
        end
        for i = 1:n-2
            # pivot or not?
            if !pivot || abs(d[i]) >= abs(dl[i])
                # No interchange
                if d[i] != 0
                    fact = dl[i]/d[i]
                    dl[i] = fact
                    d[i+1] -= fact*du[i]
                    du2[i] = 0
                end
            else
                # Interchange
                fact = d[i]/dl[i]
                d[i] = dl[i]
                dl[i] = fact
                tmp = du[i]
                du[i] = d[i+1]
                d[i+1] = tmp - fact*d[i+1]
                du2[i] = du[i+1]
                du[i+1] = -fact*du[i+1]
                ipiv[i] = i+1
            end
        end
        if n > 1
            i = n-1
            if !pivot || abs(d[i]) >= abs(dl[i])
                if d[i] != 0
                    fact = dl[i]/d[i]
                    dl[i] = fact
                    d[i+1] -= fact*du[i]
                end
            else
                fact = d[i]/dl[i]
                d[i] = dl[i]
                dl[i] = fact
                tmp = du[i]
                du[i] = d[i+1]
                d[i+1] = tmp - fact*d[i+1]
                ipiv[i] = i+1
            end
        end
        # check for a zero on the diagonal of U
        for i = 1:n
            if d[i] == 0
                info = i
                break
            end
        end
    end
    LU{T,Tridiagonal{T}}(A, ipiv, convert(BlasInt, info))
end
factorize(A::Tridiagonal) = lufact(A)

# See dgtts2.f
function A_ldiv_B!{T}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat)
    n = size(A,1)
    n == size(B,1) || throw(DimensionMismatch(""))
    nrhs = size(B,2)
    dl = A.factors.dl
    d = A.factors.d
    du = A.factors.du
    du2 = A.factors.du2
    ipiv = A.ipiv
    @inbounds begin
        for j = 1:nrhs
            for i = 1:n-1
                ip = ipiv[i]
                tmp = B[i+1-ip+i,j] - dl[i]*B[ip,j]
                B[i,j] = B[ip,j]
                B[i+1,j] = tmp
            end
            B[n,j] /= d[n]
            if n > 1
                B[n-1,j] = (B[n-1,j] - du[n-1]*B[n,j])/d[n-1]
            end
            for i = n-2:-1:1
                B[i,j] = (B[i,j] - du[i]*B[i+1,j] - du2[i]*B[i+2,j])/d[i]
            end
        end
    end
    return B
end

function At_ldiv_B!{T}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat)
    n = size(A,1)
    n == size(B,1) || throw(DimentionsMismatch(""))
    nrhs = size(B,2)
    dl = A.factors.dl
    d = A.factors.d
    du = A.factors.du
    du2 = A.factors.du2
    ipiv = A.ipiv
    @inbounds begin
        for j = 1:nrhs
            B[1,j] /= d[1]
            if n > 1
                B[2,j] = (B[2,j] - du[1]*B[1,j])/d[2]
            end
            for i = 3:n
                B[i,j] = (B[i,j] - du[i-1]*B[i-1,j] - du2[i-2]*B[i-2,j])/d[i]
            end
            for i = n-1:-1:1
                if ipiv[i] == i
                    B[i,j] = B[i,j] - dl[i]*B[i+1,j]
                else
                    tmp = B[i+1,j]
                    B[i+1,j] = B[i,j] - dl[i]*tmp
                    B[i,j] = tmp
                end
            end
        end
    end
    return B
end

# Ac_ldiv_B!{T<:Real}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) = At_ldiv_B!(A,B)
function Ac_ldiv_B!{T}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat)
    n = size(A,1)
    n == size(B,1) || throw(DimentionsMismatch(""))
    nrhs = size(B,2)
    dl = A.factors.dl
    d = A.factors.d
    du = A.factors.du
    du2 = A.factors.du2
    ipiv = A.ipiv
    @inbounds begin
        for j = 1:nrhs
            B[1,j] /= conj(d[1])
            if n > 1
                B[2,j] = (B[2,j] - conj(du[1])*B[1,j])/conj(d[2])
            end
            for i = 3:n
                B[i,j] = (B[i,j] - conj(du[i-1])*B[i-1,j] - conj(du2[i-2])*B[i-2,j])/conj(d[i])
            end
            for i = n-1:-1:1
                if ipiv[i] == i
                    B[i,j] = B[i,j] - conj(dl[i])*B[i+1,j]
                else
                    tmp = B[i+1,j]
                    B[i+1,j] = B[i,j] - conj(dl[i])*tmp
                    B[i,j] = tmp
                end
            end
        end
    end
    return B
end

/(B::AbstractMatrix,A::LU) = At_ldiv_Bt(A,B).'

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# include("functions_align.jl")
# include("functions_collapse.jl")

### --------- functions for additive addition of new seqs

# struct redundancy_set
#     tip::String
#     redundants::Array{String, 1}
# end

function m_n_reduced_is_same(nuc_bit_array_A, nuc_bit_array_B, ref_array)

    width = size(nuc_bit_array_A, 2)
    height = size(nuc_bit_array_B, 2)

    same_array = trues(height, width)

    # println("number of threads for sequence comparisons: ", Threads.nthreads())

    difference_list_A = get_list_of_differences(nuc_bit_array_A, ref_array)
    difference_list_B = get_list_of_differences(nuc_bit_array_B, ref_array)

    for a in 1:size(nuc_bit_array_A, 2)

        for b in 1:size(nuc_bit_array_B, 2)

            start_again = true

            for r in difference_list_A[a]
                @inbounds x = nuc_bit_array_A[r,a]
                @inbounds y = nuc_bit_array_B[r,b]
                different = (x & y) < 16

                # store the different True/False result in an array:
                if different
                    same_array[b,a] = !different
                    start_again = false
                    break
                end
            end

            # need to check the second seq's sites too, if we get this far:
            if start_again
                for r in difference_list_B[b]
                    if r in difference_list_A[a]
                        continue
                    end

                    @inbounds x = nuc_bit_array_A[r,a]
                    @inbounds y = nuc_bit_array_B[r,b]
                    different = (x & y) < 16

                    # store the different True/False result in an array:
                    if different
                        same_array[b,a] = !different
                        break
                    end
                end
            end
        end
    end

    return same_array

end

function read_relationships(filepath)
    # filepath should point to a tip_2_redundants.csv format file

    all_seqs = Array{String,1}()
    sets = Dict{String, Array{String, 1}}()

    first = true

    open(filepath, "r") do io
        while !eof(io)

            if first
                readline(io)
                first = false
                continue
            end

            line = readline(io)

            fields = split(line, ",")
            tip = fields[1]

            push!(all_seqs, tip)

            if length(fields) != 2
                e = error("badly formed additions file")
                throw(e)
            elseif fields[2] == ""
                redundants = Array{String,1}()
            else
                redundants = split(fields[2], "|")
                for redundant in redundants
                    push!(all_seqs, redundant)
                end
            end

            sets[tip] = redundants
        end
    end

    all_seqs = Set(all_seqs)

    return sets, all_seqs
end

function update_relationships_based_on_appearance(sets, all_seqs, fasta_IDs)
    #=
    sets = a dict of string(tip):array of strings(redundants)
    all_seqs = a set of strings (all the sequences from the relationships file)
    fasta_IDs = an array of strings which are the IDs for all the seqs in the fasta file

    A function to update the relationships that have been read in from the
    previous "tip_to_redundants.csv" file, to deal with the case where sequences
    that were in the previous dataset are missing in the new dataset.

    The following things should happen:

    1) for sets with tips that are not in the new dataset, a new tip should be chosen
       from within that set. If this is the case this redundant also needs to be removed
       from any other sets where it is a redundant (as would happen in the full routine).
       If the set has no other sequences, it can be thrown out.
    2) for sets with redundants that are not in the new dataset, those redundants
       should be stripped out of the sets.
    =#

    # make a set of the fasta IDs for faster lookup
    fasta_ID_set = Set(fasta_IDs)

    # sets whose tips is missing in these data. will check that all their
    # redundants are not redundants in sets_that_dont_need_new_tips already
    # (unchanged from the input data), because if so, can throw out this
    # whole set
    sets_that_need_new_tips = Dict{String, Array{String, 1}}()
    sets_that_dont_need_new_tips = Dict{String, Array{String, 1}}()

    # make a set of redundants that have become tips and therefore need to be taken
    # out of any other sets that they are in.

    # step 1: sort out sets that need new tips (or could potentially be thrown away)
    for (key, value) in sets
        # if this key doesn't exist any more
        if !(key in fasta_ID_set)
            # there are any redundants, flag for followup (else ignore)
            if length(value) >= 1
                sets_that_need_new_tips[key] = value
            end
        else
            sets_that_dont_need_new_tips[key] = value
        end
    end

    # step 2:
    # a) get the set of all the redundant seqs in the unchanged data:
    # make a set of the current redundants for faster lookup
    unchanged_set_seqs = Set{String}()
    for (key, value) in sets_that_dont_need_new_tips
        for v in value
            if v in fasta_ID_set
                push!(unchanged_set_seqs, v)
            end
        end
    end

    # b) find new tips for sets whose tip is missing, or throw them out if
    # they're already represented
    redundants_that_are_now_tips = Set{String}()
    sets_with_new_tips = Dict{String, Array{String, 1}}()
    for (key, value) in sets_that_need_new_tips
        # key doesn't exist in the new dataset
        for v in value
            # v might not:
            if !(v in fasta_ID_set)
                continue
            # or it might already be represented elsewhere:
            elseif v in unchanged_set_seqs
                continue
            # otherwise just
            else
                sets_with_new_tips[v] = setdiff(value, v)
                push!(redundants_that_are_now_tips, v)
                break
            end
        end
    end

    # step 3: merge the unchanged and the new sets
    new_sets_temp = Dict{String, Array{String, 1}}()

    for (key, value) in sets_with_new_tips
        new_sets_temp[key] = value
    end

    for (key, value) in sets_that_dont_need_new_tips
        new_sets_temp[key] = value
    end

    # step 4: remove redundants that are missing or have become tips
    # and populate the new set of old seqs at the same time. ???

    new_sets = Dict{String, Array{String, 1}}()
    new_all_seqs = Set{String}()

    for (key, value) in new_sets_temp
        push!(new_all_seqs, key)

        newredundants = Array{String, 1}()
        for redundant in value
            if !(redundant in redundants_that_are_now_tips)
                if redundant in fasta_ID_set
                    push!(newredundants, redundant)
                    push!(new_all_seqs, redundant)
                end
            end
        end

        new_sets[key] = newredundants
    end

    return new_sets, new_all_seqs
end

function update_relationships_based_on_sequences(sets, all_seqs, nucleotide_array, fasta_IDs, ref_array)
    #=
    sets = a dict of string(tip):array of strings(redundants)
    all_seqs = a set of strings (all the sequences from the relationships file)
    nucleotide_array = the alignment
    fasta_IDs = an array of strings which are the IDs for all the seqs in the fasta file

    A function to update the relationships that have been read in from the
    previous "tip_to_redundants.csv" file, to deal with the case where sequences
    that were in the previous dataset are different in the new dataset (the only
    way you can tell is by testing the previous sets for still being good sets).

    check that the sets are still good sets. TODO im not sure how to deal with the situation
    where they are not. bad apples (sequences) can be pushed into the new set, which means that
    they will get reassigned additively, but im not sure how to identify the bad apples in
    the first place. I could throw the whole set out. - this might be easiest?
    =#

    # need to subset the nucleotide array, so need to know what columns to use:
    ID_to_column_dict = get_id_to_column_dict(fasta_IDs) # returns a dict fasta IDs to their column in the alignment

    new_sets = Dict{String, Array{String, 1}}()
    new_all_seqs = Set{String}()

    for (key, value) in sets
        # an array containing all the seq names for this set:
        IDs = Array{String, 1}()
        push!(IDs, key)
        for v in value
            push!(IDs, v)
        end

        # an array containing all the column numbers for this set:
        colmns = Array{Int32, 1}()
        for ID in IDs
            push!(colmns, ID_to_column_dict[ID])
        end

        # now check the set:
        same_array = reduced_is_same(view(nucleotide_array, :, colmns), ref_array) # returns the same_array

        set_is_good = check_the_view_is_good(same_array) # returns true/false the same_array is ok

        if set_is_good
            # if the set is still fine, we can keep it
            new_sets[key] = value

            # and we add its members to the set of all sequences
            push!(new_all_seqs, key)
            for v in value
                push!(new_all_seqs, v)
            end
        else
            # Otherwise, we don't keep it. But we don't need to strip its members from other sets,
            # because they will be tested anyway...
            continue
        end

    end

    return new_sets, new_all_seqs
end

#
function split_data(nuc_bit_array, fasta_ID_array, old_IDs)

    old_indices = Array{Int, 1}()
    sizehint!(old_indices, length(old_IDs))
    new_indices = Array{Int, 1}()
    sizehint!(new_indices, length(old_IDs))

    for (i, ID) in enumerate(fasta_ID_array)
        if ID in old_IDs
            push!(old_indices, i)
        else
            push!(new_indices, i)
        end
    end

    old_array = nuc_bit_array[:,old_indices]
    new_array = nuc_bit_array[:,new_indices]

    old_IDs = fasta_ID_array[old_indices]
    new_IDs = fasta_ID_array[new_indices]

    return old_array, new_array, old_IDs, new_IDs
end

# function get_set_dict(old_sets)
#     set_dict = Dict{String, Array{String, 1}}()
#     for set in old_sets
#         set_dict[set.tip] = set.redundants
#     end
#
#     return set_dict
# end

# # comparison between two vectors (at sites that are diff from ref in either)
# function check_identity(nuc_bit_V_A, nuc_bit_V_B, diff_list_A, diff_list_B)
#
#     start_again = true
#     all_same = true
#
#     for r in diff_list_A
#         @inbounds x = nuc_bit_V_A[r]
#         @inbounds y = nuc_bit_V_B[r]
#         different = (x & y) < 16
#
#         # store the different True/False result in an array:
#         if different
#             all_same = false
#             start_again = false
#             break
#         end
#     end
#
#     # need to check the second seq's sites too, if we get this far:
#     if start_again
#         for r in diff_list_B
#             if r in diff_list_A
#                 continue
#             end
#
#             @inbounds x = nuc_bit_V_A[r]
#             @inbounds y = nuc_bit_V_B[r]
#             different = (x & y) < 16
#
#             # store the different True/False result in an array:
#             if different
#                 all_same = false
#                 break
#             end
#         end
#     end
#
#     return all_same
# end

function assign_new_seqs(m_n_pairwise_identity, unplaced_pairwise_identity, old_IDs, new_IDs, old_sets)

    old_ID_set = Set(old_IDs)

    # old_set_dict = get_set_dict(old_sets)
    old_ID_to_old_cl_dict = get_id_to_column_dict(old_IDs)
    new_ID_to_new_cl_dict = get_id_to_column_dict(new_IDs)

    placed_new_tips = Array{String, 1}()

    for (i, ID) in enumerate(new_IDs)

        for tip in keys(old_sets)

            tip_column = old_ID_to_old_cl_dict[tip]

            tip_same = m_n_pairwise_identity[tip_column, i]

            if tip_same
                 # check all the other members of the set...

                 set_same = true

                 other_members = old_sets[tip]

                 for member in other_members
                     # possible that this member is a new one (that has already been added to this set),
                     # so need to check what to do here:
                     if member in old_ID_set
                         redundant_column = old_ID_to_old_cl_dict[member]
                         set_same = m_n_pairwise_identity[redundant_column, i]
                         if ! set_same
                             break
                         end
                     else
                         # CHANGE BELOW - TO LOOK UP IN NEW NUC ARRAY
                         redundant_column = new_ID_to_new_cl_dict[member]
                         set_same = unplaced_pairwise_identity[redundant_column, i]
                         if ! set_same
                             break
                         end
                     end
                 end

                 if set_same
                     # update the sets
                     old_sets[tip] = vcat(old_sets[tip], ID)
                     #
                     push!(placed_new_tips, ID)
                     #
                 end
            end

        end
    end

    placed_new_tips = Set(placed_new_tips)

    unplaced_new_tip_indices = Array{Int, 1}()
    for (i, tip) in enumerate(new_IDs)
        if ! (tip in placed_new_tips)
            push!(unplaced_new_tip_indices, i)
        end
    end

    updated_sets = old_sets

    return updated_sets, unplaced_new_tip_indices
end

function TEST_get_dict_levels(relationships_dict)
    s = Set{String}()
    for (key, value) in relationships_dict
        push!(s, key)
        for v in value
            push!(s, v)
        end
    end
    return s
end

function TEST_get_sets_from_relationships(fasta_IDs, relationships_dict)

    col_dict = get_id_to_column_dict(fasta_IDs)

    A = Array{BitSet,1}()

    for (key, value) in relationships_dict
        temp_A = BitSet()
        push!(temp_A, col_dict[key])
        for v in value
            push!(temp_A, col_dict[v])
        end

        push!(A, temp_A)
    end

    return A
end

#
# function assign_new_seqs(old_nuc_bit_array, new_nuc_bit_array, old_IDs, new_IDs, old_sets, ref_array)
#
#     # old_set_dict = get_set_dict(old_sets)
#     old_ID_to_old_cl_dict = get_id_to_column_dict(old_IDs)
#     new_ID_to_new_cl_dict = get_id_to_column_dict(new_IDs)
#
#     difference_list_new = get_list_of_differences(new_nuc_bit_array, ref_array)
#     difference_list_old = get_list_of_differences(old_nuc_bit_array, ref_array)
#
#     # updated_sets = Array{redundancy_set, 1}()
#
#     placed_new_tips = Array{String, 1}()
#
#     for (i, ID) in enumerate(new_IDs)
#
#         for tip in keys(old_sets)
#
#             tip_column = old_ID_to_old_cl_dict[tip]
#
#             tip_same = check_identity(old_nuc_bit_array[:,tip_column], new_nuc_bit_array[:,i], difference_list_old[tip_column], difference_list_new[i])
#
#             if tip_same
#                  # check all the other members of the set...
#
#                  set_same = true
#
#                  other_members = old_sets[tip]
#
#                  for member in other_members
#                      # possible that this member is a new one (that has already been added to this set),
#                      # so need to check what to do here:
#                      if member in old_IDs
#                          redundant_column = old_ID_to_old_cl_dict[member]
#                          set_same = check_identity(old_nuc_bit_array[:,redundant_column], new_nuc_bit_array[:,i], difference_list_old[redundant_column], difference_list_new[i])
#                          if ! set_same
#                              break
#                          end
#                      else
#                          # CHANGE BELOW - TO LOOK UP IN NEW NUC ARRAY
#                          redundant_column = new_ID_to_new_cl_dict[member]
#                          set_same = check_identity(new_nuc_bit_array[:,redundant_column], new_nuc_bit_array[:,i], difference_list_new[redundant_column], difference_list_new[i])
#                          if ! set_same
#                              break
#                          end
#                      end
#                  end
#
#                  if set_same
#                      # update the sets
#                      old_sets[tip] = vcat(old_sets[tip], ID)
#                      #
#                      push!(placed_new_tips, ID)
#                      #
#                  end
#             end
#
#         end
#     end
#
#     placed_new_tips = Set(placed_new_tips)
#
#     unplaced_new_tip_indices = Array{Int, 1}()
#     for (i,tip) in enumerate(new_IDs)
#         if ! (tip in placed_new_tips)
#             push!(unplaced_new_tip_indices, i)
#         end
#     end
#
#     updated_sets = old_sets
#
#     unplaced_nuc_bit_array = new_nuc_bit_array[:,unplaced_new_tip_indices]
#     unplaced_IDs = new_IDs[unplaced_new_tip_indices]
#
#     return updated_sets, unplaced_nuc_bit_array, unplaced_IDs
# end

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
using Distributed
#addprocs(2)

using TaskMaster

@everywhere f(x) = x^2

@info "Testing evaluate"

master = WorkMaster(f)
learner = IgnorantLearner(1:10)
loop = Loop(master,learner)
evaluate!(loop)

@info "Testing evaluate with source"

master = WorkMaster(f)
learner = IgnorantLearner(1:10)
loop = Loop(master,learner)
output1 = evaluate!(loop,1:6) ### In a way looks like a transducer

@info "Testing the debugger for Learner"

master = HistoryMaster(output1,length(master.slaves))
learner = IgnorantLearner(1:10)
loop = Loop(master,learner,loop->println("Learner state $(loop.learner.state)"))
output2 = evaluate!(loop,1:6)

@info "Testing evaluate with stopping condition"

master = WorkMaster(f)
learner = IgnorantLearner(1:10)
loop = Loop(master,learner)
evaluate!(loop,l->l.state==4)

@info "Testing capturing and releasing of the slave"

@everywhere f(x) = x^2
learner = IgnorantLearner(1:10)
master = WorkMaster(f,WorkerPool())

captureslave!(2,f,master)
captureslave!(3,f,master)
@show releaseslave!(master)
@show releaseslave!(master)

@info "Testing a more coincise API"

@everywhere f(x) = x^2

learner = IgnorantLearner(1:10)
results = play!(learner,f,x->x.state==4)

@info "Testing replay of the learner"

learner = IgnorantLearner(1:10)
replay!(learner, results)
play!(learner,f,x->x.state==8)

@info "Success!!!"


text:
using LinearAlgebra: norm

"""
    DiatomicHarmonic(r₀=1.0)

Harmonic interaction between two particles.
"""
Parameters.@with_kw struct DiatomicHarmonic <: AdiabaticModel
    r₀::Float64 = 1.0
end

NQCModels.ndofs(::DiatomicHarmonic) = 3

function NQCModels.potential(model::DiatomicHarmonic, R::AbstractMatrix)
    return (norm(R[:,1] .- R[:,2]) - model.r₀)^2 / 2
end

function NQCModels.derivative!(model::DiatomicHarmonic, D::AbstractMatrix, R::AbstractMatrix)
    diff = R[:,1] .- R[:,2]
    leng = norm(diff)
    D .= (leng - model.r₀) / leng .* diff
    D[:,2] .*= -1
end


text:
#bellcurve2.jl
#Algorithm 2.1 of Random Eigenvalues by Alan Edelman

#Experiment:    Generate random samples from the normal distribution
#Plot:          Histogram of random samples
#Theory:        The normal distribution curve

## Experiment
t = 10000     # Trials
dx = 0.25       # binsize
v = randn(t,2)  # samples
x = -4:dx:4   # range

## Plot
using Winston
module Plot3D
    using Winston
    # Warning: the 3D interface is very new at the time
    # of writing this and is likely to not be available
    # and/or have changed
    if Winston.output_surface == :gtk
        include(Pkg.dir("Winston","src","canvas3d_gtk.jl"))
    else
        include(Pkg.dir("Winston","src","canvas3d.jl"))
    end
end
edg1, edg2, count = hist2(v, x-dx/2)
#cntrs = ((first(x)+dx/2), (last(x)-dx/2), length(x)-1)
#x,y = meshgrid(cntrs, cntrs)
cntrs = (first(x)+dx/2):dx:(last(x)-dx/2)
x = Float64[x for x = cntrs, y = cntrs]
y = Float64[y for x = cntrs, y = cntrs]
Plot3D.plot3d(Plot3D.surf(x,count./(t*dx^2) *50,y))
Plot3D.plot3d(Plot3D.surf(x,exp(-(x.^2+y.^2)/2)/(2*pi) *50,y))


text:
using DataFrames
using Statistics
using Rmath
using MultipleTesting
using Optim

################################################################################
# Functions that normalise data like M3Drop expects
################################################################################

#Find no. detected or undetected genes/cell. This is a faff in Julia :(
function FindNumDetectedGenesPerCell(countsMatrix::DataFrame, Detected::Bool)

    #This allows us to use the same function for detected + undetected genes
    if Detected == true
        GreaterThanZero = colwise(x -> x .> 0 , countsMatrix[:,2:ncol(countsMatrix)])
    else
        GreaterThanZero = colwise(x -> x .== 0 , countsMatrix[:,2:ncol(countsMatrix)])
    end

    tmp_arr = Float64[]

    for arr in GreaterThanZero
        push!(tmp_arr, sum(arr))
    end

    outputDf = DataFrame(Cells = names(countsMatrix)[2:ncol(countsMatrix)],
    DetectedGenes = tmp_arr)

    return outputDf

end

#Find mean
function FindMean(DetectedGenesDf::DataFrame)
    return mean(DetectedGenesDf.DetectedGenes)
end

#Find standard deviation
function FindStd(DetectedGenesDf::DataFrame)
    return std(DetectedGenesDf.DetectedGenes)
end

#Remove low quality cells based on no. undetected genes
function RemoveLowQuality(DetectedGenesDf, mean, sd)

    pnorms = (DetectedGenesDf.DetectedGenes .- mean) ./ sd
    pnorms = pnorm.(pnorms, false)
    padjust = adjust(pnorms, BenjaminiHochberg())

    DetectedGenesDf.pvals = padjust
    filtered_df = DetectedGenesDf[DetectedGenesDf.pvals .> 0.05, :]

    return(filtered_df[:, 1:2])
end

#Only keep filtered cells
function KeepFiltered(DetectedGenesDf::DataFrame, countsMatrix::DataFrame)

    arr = [Symbol("Isoforms")]
    append!(arr, DetectedGenesDf.Cells)

    return countsMatrix[arr]

end

#Helper function to make an empty counts matrix
function MakeEmptyDf(countsMatrix)

    #Make an empty df to store results in
    namelist = names(countsMatrix)
    df = DataFrame()

    for (i, name) in enumerate(namelist)
        if i == 1
            df[name] = String[]
        else
            df[name] = Float64[]
        end
    end

    return df
end

#Remove undetected genes (separate function)
function RemoveUndetectedGenes(countsMatrix::DataFrame)

    #Make an empty df to store results in
    df = MakeEmptyDf(countsMatrix)

    #Iterate over rows, keep rows where isoforms is expressed in >3 cells
    for row in eachrow(countsMatrix)

        if sum(row[2:ncol(countsMatrix)] .> 0 ) > 3
            push!(df, row)
        end
    end

    return df
end

#Convert to cpm
function ConvertToCpm(countsMatrix::DataFrame)

    #Find cpms
    cpms = colwise(x -> (x / sum(x)) * 1000000, countsMatrix[2:ncol(countsMatrix)])

    #Make back into dataframe...
    namelist = names(countsMatrix[2:ncol(countsMatrix)])
    df = DataFrame(Isoforms = countsMatrix.Isoforms)

    for (i, name) in enumerate(namelist)
        df[name] =  cpms[i]
    end

    return df

end

#Remove lowly expressed genes (function)
function RemoveLowlyExpressedGenes(countsMatrix::DataFrame)

    #Make an empty df to store results in
    df = MakeEmptyDf(countsMatrix)

    #Iterate over rows, keep rows where isoforms isn't lowly expressed
    for row in eachrow(countsMatrix)

        if sum(row[2:ncol(countsMatrix)]) > 0.00001
            push!(df, row)
        end
    end

    return df
end

#Make a function to join everything up
function CleanDataForMM(countsMatrix::DataFrame)

    #Find no. undetected genes/cell
    NumZero = FindNumDetectedGenesPerCell(countsMatrix, false)

    #Find mean and sd no. undetected genes per cell
    Mu = FindMean(NumZero)
    Sigma = FindStd(NumZero)

    #Remove low quality cells based on no. undetected genes per cell
    highQualityCells = RemoveLowQuality(NumZero, Mu, Sigma)

    #Only keep high quality cells in countsMatrix
    outputDf = KeepFiltered(highQualityCells, countsMatrix)

    #Remove undetected genes
    outputDf = RemoveLowlyExpressedGenes(outputDf)

    #Convert counts to cpm
    outputDf = ConvertToCpm(outputDf)

    #Remove lowly expressed genes
    outputDf = RemoveLowlyExpressedGenes(outputDf)

    return outputDf

end

################################################################################
# Functions that calculate Km
################################################################################

#Calculate p (M-M parameter)
function Findp(countsMatrix::DataFrame)

    df = DataFrame()
    df[:Isoforms] = countsMatrix[:Isoforms]
    arr = Float64[]

    for row in eachrow(countsMatrix)
        p = 1 - (sum(row[2:ncol(countsMatrix)] .> 0) / (ncol(countsMatrix) - 1))
        push!(arr, p)
    end

    df[:p] = arr
    return df

end

#Calculate s (M-M parameter)
function Finds(df::DataFrame, countsMatrix::DataFrame)

    arr = Float64[]

    for row in eachrow(countsMatrix)
        s = mean(row[2:ncol(countsMatrix)])
        push!(arr, s)
    end

    df[:s] = arr
    return df
end

#Function that makes function to calculate log likelihood
function LogLikelihood(p, s)
    f(params) = -sum(dnorm.(p .- (1.0 .- (s ./ (params[1] .+ s))), 0.0, params[2], true))
    f
end

#MLE
function MaxLikelihoodEst(p, s)

    #Check p and s make sense
    if length(p) != length(s)
        throw(ArgumentError("Length of p and s don't match, something went wrong\n"))
    end

    #NB. MLE starts get inconsistent below 5,000 genes
    if length(p) < 5000
        println("Warning: Calculating Km for less than 5000 isoforms... MLE
        estimate of Km may not be accurate.")
    end

    #This line basically internalises p and s - useful for optimize function below
    LogLikelihoodFn = LogLikelihood(p,s)

    #Perform MLE and return parameters
    res = optimize(LogLikelihoodFn, [3.0, 0.25], method = NelderMead())
    parameters = Optim.minimizer(res)
    return parameters
end

#Function to join up all the other functions
function FindKm(countsMatrix::DataFrame)

    df = Findp(countsMatrix)
    df = Finds(df, countsMatrix)
    parameters = MaxLikelihoodEst(df.p, df.s)
    return parameters
end

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
using Unitful, UnitfulAssets

module ProductionUnits
    using Unitful
    using Unitful: @unit
    @unit tshirt "tshirt" TShirt 1 false
    @unit worker "worker" Worker 1 false
end

Unitful.register(ProductionUnits);

fabric = 15u"USD"/8u"yd"/44u"inch"
dyes = 20u"USD/lb"
fixer = 8u"USD"/5u"lb"
thread = 19u"USD"/1000u"yd"

# raw_material_cost = 1.6u"m^2" * fabric + 2u"oz" * dyes + 1u"oz" * fixer + 48u"yd" * thread

raw_price_table = (fabric, dyes, fixer, thread)

"""
    raw_material(n::Unitful.Quantity)

Return the amount of each raw material needed to manufacture `n` T-shirts.

The argument `n` must be given in `tshirt` units.

Returns a tuple with the following quantities, respectively:
* The necessary amount of cotton fabric.
* The necessary amount of dye.
* The necessary amount of fixer.
* The necessary amount of thread.

"""
raw_material(n::Unitful.Quantity) = (1.6u"m^2" * n / u"tshirt", 2u"oz" * n / u"tshirt", 1u"oz" * n / u"tshirt", 48u"yd" * n / u"tshirt")

raw_material_price = sum(raw_material(1u"tshirt") .*  raw_price_table)

eletricity_price = 0.13u"USD/kW/hr"
labor_price = 10.50u"USD/worker/hr"
fixed_cost = 24000u"USD/yr"

"""
    manufacturing_cost(n::Unitful.Quantity, t::Unitful.Quantity, tlim::Unitful.Quantity=40u"hr/worker/wk")

Return the cost of manufacturing `n` T-shirts during a time period `t`.

The argument `n` must be given in `tshirt` units, and `t`, in time units.
The optional argument `tlim` is the time limit of work per worker, which
defaults to `40u"hr/worker/wk"`.

Return a tuple with the following quantities, respectively:
* The cost of the production, in US Dollars.
* The cost per T-shirt.
* The number of labor hours required to produce `n` t-shirts.
* The minimum number of workers considering the limit given by `tlim`.
* The eletricity required for the whole manufacturing process.
"""
function production_cost(n::Unitful.Quantity, t::Unitful.Quantity, tlim::Unitful.Quantity=40u"hr/worker/wk")
    labor_hours = 2u"hr/tshirt" * n
    eletricity_spent = 2u"kW * hr/tshirt" * n
    total_cost = n * raw_material_price + labor_hours * labor_price + eletricity_spent * eletricity_price + fixed_cost * t
    cost_per_tshirt = total_cost / n
    min_num_workers = Int(ceil(labor_hours/tlim/t)) * u"worker"
    return total_cost, cost_per_tshirt, labor_hours, min_num_workers, eletricity_spent
end

production_cost(50u"tshirt", 1u"wk")
raw_material(50u"tshirt")

production_cost(2000u"tshirt", 30u"d", 44u"hr/worker/wk")
raw_material(2000u"tshirt")


text:
module WhiskerTracking

#Standard Library
using Statistics,Random,Distributed,SharedArrays,DelimitedFiles,LinearAlgebra, Libdl, Dates

#Deep Learning Libraries
#using CuArrays, CuArrays.CUFFT, CUDAnative, Knet
using CUDA, Knet

using Gtk.ShortNames, Cairo, Images, ImageFiltering, MAT, JLD2, Interpolations, DSP,
FFTW, IterTools, FFMPEG, StackedHourglass, Polynomials, StatsBase

if VERSION > v"0.7-"
    const Void = Nothing
    const setproperty! = set_gtk_property!
    const getproperty = get_gtk_property
    const unshift! = pushfirst!
    const is_windows() = Sys.iswindows()
    const is_unix() = Sys.isunix()
    const find = findall
    const indmax = argmax
end

include("config.jl")

function __init__()
    ccall((:Load_Params_File,libwhisk_path),Int32,(Cstring,),jt_parameters)
end

include("types.jl")

include("deep_learning/helper.jl")
include("deep_learning/save_load.jl")
include("deep_learning/prediction.jl")
include("deep_learning/config.jl")
include("deep_learning/training.jl")

include("gui.jl")
include("janelia_tracker.jl")
include("processing.jl")
include("image_processing.jl")
include("analysis.jl")
include("save_load.jl")
include("discrete.jl")

include("gui/manual.jl")
include("gui/analog.jl")
include("gui/whisker_pad.jl")
include("gui/discrete.jl")
include("gui/mask.jl")
include("gui/pole.jl")
include("gui/view.jl")
include("gui/tracing.jl")
include("gui/image.jl")
include("gui/janelia.jl")
include("gui/export.jl")
include("gui/contact.jl")
include("gui/deeplearning.jl")

include("drawing_tools/draw.jl")
include("drawing_tools/erase.jl")
include("drawing_tools/shapes.jl")

include("analysis/follicle.jl")
include("analysis/kinematics.jl")
include("analysis/forces.jl")
include("analysis/touch.jl")
include("analysis/math.jl")

include("precompile.jl")
_precompile()

end


text:
module LSDescent

export bfgs, bfgs_Stop, bfgs_StopLS
export L_bfgs, L_bfgs_Stop, L_bfgs_StopLS
export M_bfgs, M_bfgs_Stop, M_bfgs_StopLS
export Newton_StopLS, Newton_Stop, Newton_Spectral
export CG_generic, formula_FR, formula_PR, formula_HS, formula_HZ
export CG_FR, CG_PR, CG_HS, CG_HZ

using LinearAlgebra

using NLPModels
using OneDmin

using LinearOperators
using Logging
using SolverCore   # Pour avoir les utilitaires d'affichage log_header et log_row

using Stopping

include("Newton/NewtonStopLS.jl")
include("Newton/NewtonStop.jl")
include("Newton/NewtonSolver.jl")

export Newton_G_StopLS, Newton_G_Stop
include("Newton/Newton_G_StopLS.jl")
include("Newton/Newton_G_Stop.jl")

include("Newton/HessianDense.jl")
include("Newton/HessianOp.jl")
include("Newton/HessianSparse.jl")

include("Newton/NewtonSpectralAbs.jl")

include("Newton/ldlt_symm.jl")
export ldlt_symm
include("Newton/NewtonLDLtAbs.jl")
include("Newton/NewtonCG.jl")

include("BFGS/bfgsSolver.jl")
include("BFGS/bfgsStop.jl")
include("BFGS/bfgsStopLS.jl")
include("BFGS/FormuleN2.jl")

include("BFGS/Type.jl")
include("BFGS/FormuleN2Op.jl")

include("ConjugateGradient/CG_generic.jl")
include("ConjugateGradient/formulae.jl")

include("ConjugateGradient/CG_FR.jl")
include("ConjugateGradient/CG_PR.jl")
include("ConjugateGradient/CG_HS.jl")
include("ConjugateGradient/CG_HZ.jl")

end # module

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
using Cthulhu
using Test

function process(@nospecialize(f), @nospecialize(TT); optimize=true)
    mi = Cthulhu.first_method_instance(f, TT)
    (ci, rt, slottypes) = Cthulhu.do_typeinf_slottypes(mi, optimize, Cthulhu.current_params())
    Cthulhu.preprocess_ci!(ci, mi, optimize)
    ci, mi, rt, slottypes
end

function find_callsites_by_ftt(@nospecialize(f), @nospecialize(TT); optimize=true)
    ci, mi, _, slottypes = process(f, TT; optimize=optimize)
    callsites = Cthulhu.find_callsites(ci, mi, slottypes)
end

# Testing that we don't have spurious calls from `Type`
let callsites = find_callsites_by_ftt(Base.throw_boundserror, Tuple{UnitRange{Int64},Int64})
    @test length(callsites) == 1
end

function test()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

let callsites = find_callsites_by_ftt(test, Tuple{})
    @test length(callsites) == 4
end

let callsites = find_callsites_by_ftt(test, Tuple{}; optimize=false)
    @test length(callsites) == 4
end

# Check that we see callsites that are the rhs of assignments
@noinline bar_callsite_assign() = nothing
function foo_callsite_assign()
    x = bar_callsite_assign()
    x
end
let callsites = find_callsites_by_ftt(foo_callsite_assign, Tuple{}; optimize=false)
    @test length(callsites) == 1
end

@eval function call_rt()
    S = $(Core.Compiler.return_type)(+, Tuple{Int, Int})
end
let callsites = find_callsites_by_ftt(call_rt, Tuple{}; optimize=false)
    @test length(callsites) == 1
end

function call_by_apply(args...)
    identity(args...)
end
let callsites = find_callsites_by_ftt(call_by_apply, Tuple{Tuple{Int}}; optimize=false)
    @test length(callsites) == 1
end

if VERSION >= v"1.1.0-DEV.215" && Base.JLOptions().check_bounds == 0
Base.@propagate_inbounds function f(x)
    @boundscheck error()
end
g(x) = @inbounds f(x)
h(x) = f(x)

let (CI, _, _, _) = process(g, Tuple{Vector{Float64}})
    @test length(CI.code) == 3
end

let (CI, _, _, _) = process(h, Tuple{Vector{Float64}})
    @test length(CI.code) == 2
end
end

f(a, b) = a + b
let callsites = find_callsites_by_ftt(f, Tuple{Any, Any})
    @test length(callsites) == 1
    callinfo = callsites[1].info
    @test callinfo isa Cthulhu.MultiCallInfo
end

# Failed return_type
only_ints(::Integer) = 1
return_type_failure(::T) where T = Base._return_type(only_ints, Tuple{T})
let callsites = find_callsites_by_ftt(return_type_failure, Tuple{Float64}, optimize=false)
    @test length(callsites) == 1
    callinfo = callsites[1].info
    @test callinfo isa Cthulhu.ReturnTypeCallInfo
    callinfo = callinfo.called_mi
    @test callinfo isa Cthulhu.FailedCallInfo
end

# tasks
ftask() = @sync @async show(io, "Hello")
let callsites = find_callsites_by_ftt(ftask, Tuple{})
    @test !isempty(filter(c->c.info isa Cthulhu.TaskCallInfo, callsites))
end

callf(f::F, x) where F = f(x)
let callsites = find_callsites_by_ftt(callf, Tuple{Union{typeof(sin), typeof(cos)}, Float64})
    @test !isempty(callsites)
    @test first(callsites).info isa Cthulhu.MultiCallInfo
end


text:
module ObjectTile

using ..Ahorn, Maple

displayName = "Object Tiles"
group = "Placements"

drawingLayers = Ahorn.Layer[]

toolsLayer = nothing
targetLayer = nothing
entityLayer = nothing

material = -1

position = nothing

function idToDisplayName(id::Int)
    return get(Ahorn.ObjectTileNames.names, id, string(id))
end

function displayNameToId(name::String)
    for (k, v) in Ahorn.ObjectTileNames.names
        if v == name
            return k
        end
    end

    try
        return parse(Int, name)

    catch
        return -1
    end
end

function drawBrushes(layer::Ahorn.Layer, room::Maple.Room)
    if position !== nothing
        ctx = Ahorn.getSurfaceContext(layer.surface)

        Ahorn.drawObjectTile(ctx, position..., material)
    end
end

function cleanup()
    global position = nothing

    Ahorn.redrawLayer!(toolsLayer)
end

function setMaterials!(layer::Ahorn.Layer)
    materials = String[Ahorn.ObjectTileNames.names[-1]]
    objSprite = Ahorn.getSprite("tilesets/scenery", "Gameplay")

    width = objSprite.realWidth
    height = objSprite.realHeight

    cellWidth = floor(Int, width / 8)
    cellHeight = floor(Int, height / 8)

    surface = Ahorn.CairoImageSurface(zeros(UInt32, width, height), Ahorn.Cairo.FORMAT_ARGB32)
    ctx = Ahorn.CairoContext(surface)

    Ahorn.drawImage(ctx, "tilesets/scenery", 0, 0, atlas="Gameplay")

    data = surface.data

    for x in 0:cellWidth - 1
        for y in 0:cellHeight - 1
            if !all(data[x * 8 + 1:x * 8 + 8, y * 8 + 1:y * 8 + 8] .== 0)
                push!(materials, idToDisplayName(y * cellWidth + x))
            end
        end
    end

    Ahorn.setMaterialList!(sort(materials), row -> row[1] == string(material))
end

function mouseMotion(x::Number, y::Number)
    global position = (x, y)

    Ahorn.redrawLayer!(toolsLayer)
end

function leftClick(x::Number, y::Number)
    Ahorn.History.addSnapshot!(Ahorn.History.RoomSnapshot("Obj Tile Test", Ahorn.loadedState.room))

    Maple.updateTileSize!(Ahorn.loadedState.room, '0', '0')

    data = Ahorn.loadedState.room.objTiles.data
    height, width = size(data)

    dr = Ahorn.getDrawableRoom(Ahorn.loadedState.map, Ahorn.loadedState.room)
    entityLayer = Ahorn.getLayerByName(dr.layers, "entities")

    if 1 <= x <= width && 1 <= y <= height
        data[y, x] = material

        Ahorn.redrawLayer!(targetLayer)
        Ahorn.redrawLayer!(entityLayer)
    end
end

function middleClick(x::Number, y::Number)
    global material = get(Ahorn.loadedState.room.objTiles.data, (y, x), -1)

    Ahorn.persistence["objtiles_material"] = material
    Ahorn.selectMaterialList!(string(material))
end

function toolSelected(subTools::Ahorn.ListContainer, layers::Ahorn.ListContainer, materials::Ahorn.ListContainer)
    global material = get(Ahorn.persistence, "objtiles_material", -1)

    Ahorn.updateLayerList!(["entities"], 1)

    Ahorn.redrawingFuncs["tools"] = drawBrushes
    Ahorn.redrawLayer!(toolsLayer)
end

function layerSelected(list::Ahorn.ListContainer, materials::Ahorn.ListContainer, selected::String)
    global targetLayer = Ahorn.getLayerByName(drawingLayers, selected)

    setMaterials!(targetLayer)
end

function materialSelected(list::Ahorn.ListContainer, selected::String)
    selectedId = displayNameToId(selected)

    Ahorn.persistence["objtiles_material"] = selectedId

    global material = selectedId
end

function materialFiltered(list::Ahorn.ListContainer)
    targetName = idToDisplayName(material)

    Ahorn.selectRow!(list, row -> row[1] == targetName)
end

function layersChanged(layers::Array{Ahorn.Layer, 1})
    global drawingLayers = layers

    global toolsLayer = Ahorn.getLayerByName(layers, "tools")
    global entityLayer = Ahorn.getLayerByName(layers, "entities")
end

end


text:
const Path = Vector{Synset}
const SimilarityGroups = Vector{Set{Synset}}

struct BasicSynset
    words::Vector{String}
end
Base.hash(s::BasicSynset, h::UInt) = hash(s.words, h)
Base.isequal(s1::BasicSynset, s2::BasicSynset) = s1.words === s2.words

WordNet.words(b::BasicSynset) = b.words

const BasicPath = Vector{BasicSynset}
const BasicSimilarityGroups = Vector{Set{BasicSynset}}

struct Cache
    # db::WordNet.DB
    synsets::Dict{String, Vector{BasicSynset}}
    paths::Dict{String, Vector{BasicPath}}
    similarity_groups::Dict{String, Vector{BasicSimilarityGroups}}
end

function push(v::AbstractVector, x)
    result = copy(v)
    push!(result, x)
    result
end

function paths_to_synset(db::DB, synset::Synset)
    complete_paths = Vector{Path}()
    active_set = [[synset]]
    while !isempty(active_set)
        path = pop!(active_set)
        parents = WordNet.relation(db, path[end], WordNet.HYPERNYM)
        if isempty(parents)
            push!(path, WordNet.∅)
            reverse!(path)
            push!(complete_paths, path)
        else
            for parent in parents
                if parent ∉ path
                    if length(parents) == 1
                        push!(path, parent)
                        new_path = path
                    else
                        new_path = push(path, parent)
                    end
                    push!(active_set, new_path)
                end
            end
        end
    end
    complete_paths
end

similar_to(db::DB, synset::Synset) = WordNet.relation(db, synset, WordNet.SIMILAR_TO)

function similarity_groups(db::DB, synset::Synset, max_depth=10)
    explored = Set([synset])
    groups = [Set([synset])]
    for i in 1:max_depth
        active_set = groups[i]
        frontier = Set{Synset}()
        for synset in active_set
            for related in similar_to(db, synset)
                if related in explored
                    continue
                end
                push!(explored, related)
                push!(frontier, related)
            end
        end
        if isempty(frontier)
            break
        end
        push!(groups, frontier)
    end
    groups
end

function basic_stem(word)
    if endswith(word, 's')
        word[1:end-1]
    else
        word
    end
end

strip_wordnet_pointers(s::Synset, c::AbstractDict{Synset, BasicSynset}) =
    get!(() -> BasicSynset(collect(words(s))), c, s)
strip_wordnet_pointers(v::Set{Synset}, c::AbstractDict{Synset, BasicSynset}) =
    Set([strip_wordnet_pointers(x, c) for x in v])
function strip_wordnet_pointers(v::Union{
            Vector{Synset},
            Vector{Set{Synset}},
            Vector{Vector{Synset}},
            Vector{Vector{Set{Synset}}}},
        c::AbstractDict{Synset, BasicSynset})
    [strip_wordnet_pointers(x, c) for x in v]
end

function Cache()
    db = DB()
    synsets = Dict{String, Vector{BasicSynset}}()
    paths = Dict{String, Vector{BasicPath}}()
    groups = Dict{String, Vector{BasicSimilarityGroups}}()
    synset_cache = Dict{Synset, BasicSynset}()
    @showprogress "Caching similarities" for (pos, part_of_speech_synsets) in db.synsets
        for synset in values(part_of_speech_synsets)
            for word in words(synset)
                push!(get!(Vector{BasicSynset},
                           synsets,
                           basic_stem(normalize(word))),
                      strip_wordnet_pointers(synset, synset_cache))
            end
            if synset.pos ∈ ('a', 's')
                g = strip_wordnet_pointers(similarity_groups(db, synset), synset_cache)
                for word in words(synset)
                    push!(get!(Vector{BasicSimilarityGroups}, groups, basic_stem(normalize(word))),
                          g)
                end
            else
                p = strip_wordnet_pointers(paths_to_synset(db, synset), synset_cache)
                for word in words(synset)
                    append!(get!(Vector{BasicPath}, paths, basic_stem(normalize(word))), p)
                end
            end
        end
    end
    Cache(synsets, paths, groups)
end


text:
using OceanDistributions
using Test

@testset "OceanDistributions.jl" begin
    Tcensus=readoceandistribution("../examples/Tcensus.txt")
    tmp=sum(Tcensus[!,:V])
    @test isapprox(tmp,1.3349978769392906e18; atol=1e16)
end


text:
using CellListMap
using StaticArrays
import Random

#
# In this test we compute the minimum distance between two independent sets of particles,
# without periodic conditions. We need to define a special reducing function to the
# parallel computation.
#
function nearest_neighbor_nopbc(;N1=1_500,N2=1_500_000,parallel=true,x=nothing,y=nothing)

    # Particle positions
    Random.seed!(321)
    if x === nothing
        x = [ rand(SVector{3,Float64}) for i in 1:N1 ]
    else
        N1 = length(x)
    end
    if y === nothing
        y = [ rand(SVector{3,Float64}) for i in 1:N2 ]
    else
        N2 = length(y)
    end

    # Obtain one upper bound for dmin by computing one distance for each element
    # of the smallest vector
    cutoff = +Inf
    for v in x
        iy = rand(1:N2)
        cutoff = min(CellListMap.norm(v - y[iy]), cutoff)
    end

    # Define box, since no PBC are used, define sizes with limits(x,y)
    box = Box(limits(x,y), cutoff)

    # Initialize auxiliary linked lists
    cl = CellList(x, y, box, parallel=parallel)

    # Function that keeps the minimum distance
    f(i, j, d2, mind) = d2 < mind[3] ? (i, j, d2) : mind

    # We have to define our own reduce function here (for the parallel version)
    function reduce_mind(output, output_threaded)
        mind = output_threaded[1]
        for i in 2:length(output_threaded)
            if output_threaded[i][3] < mind[3]
                mind = output_threaded[i]
            end
        end
        return mind
    end

    # Initialize
    mind = (0, 0, +Inf)

    # Run pairwise computation
    mind = map_pairwise(
        (x, y, i, j, d2, mind) -> f(i, j, d2, mind),
        mind,box,cl;reduce=reduce_mind,parallel=parallel
    )

    # Take the square root of the minimum distance to return
    return (mind[1], mind[2], sqrt(mind[3]))

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## CHOLMOD
const TRUE  = Int32(1)
const FALSE = Int32(0)

## itype defines the types of integer used:
const INT     = Int32(0)  # all integer arrays are int
const INTLONG = Int32(1)  # most are int, some are SuiteSparse_long
const LONG    = Int32(2)  # all integer arrays are SuiteSparse_long

## dtype defines what the numerical type is (double or float):
const DOUBLE = Int32(0)        # all numerical values are double
const SINGLE = Int32(1)        # all numerical values are float
dtyp(::Type{Float32}) = SINGLE
dtyp(::Type{Float64}) = DOUBLE
dtyp(::Type{ComplexF32}) = SINGLE
dtyp(::Type{ComplexF64}) = DOUBLE

## xtype defines the kind of numerical values used:
const PATTERN = Int32(0)       # pattern only, no numerical values
const REAL    = Int32(1)       # a real matrix
const COMPLEX = Int32(2)       # a complex matrix (ANSI C99 compatible)
const ZOMPLEX = Int32(3)       # a complex matrix (MATLAB compatible)
xtyp(::Type{Float32})    = REAL
xtyp(::Type{Float64})    = REAL
xtyp(::Type{ComplexF32}) = COMPLEX
xtyp(::Type{ComplexF64}) = COMPLEX

## Scaling modes, selected by the scale input parameter:
const SCALAR = Int32(0)        # A = s*A
const ROW    = Int32(1)        # A = diag(s)*A
const COL    = Int32(2)        # A = A*diag(s)
const SYM    = Int32(3)        # A = diag(s)*A*diag(s)

## Types of systems to solve
const CHOLMOD_A    = Int32(0)          # solve Ax=b
const CHOLMOD_LDLt = Int32(1)          # solve LDL'x=b
const CHOLMOD_LD   = Int32(2)          # solve LDx=b
const CHOLMOD_DLt  = Int32(3)          # solve DL'x=b
const CHOLMOD_L    = Int32(4)          # solve Lx=b
const CHOLMOD_Lt   = Int32(5)          # solve L'x=b
const CHOLMOD_D    = Int32(6)          # solve Dx=b
const CHOLMOD_P    = Int32(7)          # permute x=Px
const CHOLMOD_Pt   = Int32(8)          # permute x=P'x

# Symmetry types
const EMPTY                 =-1
const MM_RECTANGULAR        = 1
const MM_UNSYMMETRIC        = 2
const MM_SYMMETRIC          = 3
const MM_HERMITIAN          = 4
const MM_SKEW_SYMMETRIC     = 5
const MM_SYMMETRIC_POSDIAG  = 6
const MM_HERMITIAN_POSDIAG  = 7

# check the size of SuiteSparse_long
if Int(ccall((:jl_cholmod_sizeof_long, :libsuitesparse_wrapper),Csize_t,())) == 4
    const SuiteSparse_long = Int32
    const IndexTypes = (:Int32,)
    const ITypes = Union{Int32}
else
    const SuiteSparse_long = Int64
    const IndexTypes = (:Int32, :Int64)
    const ITypes = Union{Int32, Int64}
end
ityp(::Type{SuiteSparse_long}) = LONG

const VTypes = Union{ComplexF64, Float64}
const VRealTypes = Union{Float64}

struct CHOLMODException <: Exception
    msg::String
end

macro isok(A)
    :($(esc(A)) == TRUE || throw(CHOLMODException("")))
end

max_stars_repo_name: sg-s/julia

max_issues_repo_name: sg-s/julia

max_forks_repo_name: sg-s/julia


text:
using Basins
using DynamicalSystems
using Test
using DifferentialEquations

@testset "All tests" begin

@testset "Test basin_discrete_map" begin
    v_bnd=[1.2727272727272727 1.1111111111111112;
1.6767676767676767 1.7575757575757576;
-0.2222222222222222 0.30303030303030304;
-0.6666666666666666 -1.595959595959596;
-0.26262626262626265 0.42424242424242425;
1.1515151515151516 1.1919191919191918;
0.5454545454545454 -1.797979797979798;
-1.4747474747474747 0.8282828282828283;
1.2727272727272727 -0.2222222222222222;
0.9494949494949495 1.5151515151515151;
-1.2323232323232323 -0.020202020202020204;
1.1111111111111112 -0.5454545454545454;
0.5858585858585859 -1.878787878787879;
0.7878787878787878 -1.3131313131313131;
-1.4747474747474747 0.9494949494949495;
1.7171717171717171 1.797979797979798;
0.5858585858585859 0.9090909090909091;
1.7171717171717171 1.7575757575757576;
-1.1111111111111112 2.0;
-0.7070707070707071 -1.5555555555555556;
1.0707070707070707 -0.6666666666666666;
1.0707070707070707 0.3838383838383838;
-1.2727272727272727 0.020202020202020204;
1.1515151515151516 1.878787878787879;
-1.393939393939394 0.7474747474747475;
-0.9090909090909091 1.4747474747474747;
0.8282828282828283 -1.3131313131313131;
1.3535353535353536 1.3131313131313131;
-0.8686868686868687 1.3131313131313131;
1.5151515151515151 1.0303030303030303;
0.9494949494949495 -0.9494949494949495;
1.5555555555555556 1.0303030303030303;
0.3434343434343434 0.46464646464646464;
-0.5454545454545454 0.6666666666666666;
1.6767676767676767 1.7575757575757576;
1.3131313131313131 1.6363636363636365;
1.4343434343434343 1.4343434343434343;
-1.2323232323232323 0.020202020202020204;
0.9090909090909091 -1.2323232323232323;
0.7474747474747475 1.1515151515151516;
-1.6363636363636365 1.5555555555555556;
1.4747474747474747 1.5555555555555556;
1.2727272727272727 0.9090909090909091;
-0.3434343434343434 0.42424242424242425;
-0.7070707070707071 -1.4343434343434343;
-1.2727272727272727 0.26262626262626265;
1.0707070707070707 1.7575757575757576;
0.46464646464646464 -2.0;
1.4343434343434343 1.5555555555555556;
0.30303030303030304 0.42424242424242425;
-0.8282828282828283 -1.2727272727272727;
-1.5151515151515151 1.0303030303030303;
1.3535353535353536 0.3838383838383838;
-0.8686868686868687 1.3535353535353536;
1.4343434343434343 0.7070707070707071;
-0.98989898989899 1.5555555555555556;
1.4747474747474747 0.6262626262626263;
1.4343434343434343 1.9595959595959596;
-0.6666666666666666 0.7878787878787878;
1.1111111111111112 -0.7070707070707071;
-1.2323232323232323 0.10101010101010101;
-1.393939393939394 0.5050505050505051;
1.3535353535353536 1.6767676767676767;
0.8686868686868687 1.3535353535353536;
0.42424242424242425 0.5050505050505051;
0.8686868686868687 1.393939393939394;
-0.5050505050505051 0.6262626262626263;
1.0707070707070707 1.6363636363636365;
1.3131313131313131 1.0303030303030303;
1.0303030303030303 -0.9494949494949495;
0.5454545454545454 -1.9191919191919191;
-1.4747474747474747 0.98989898989899;
1.1919191919191918 0.6262626262626263;
1.0707070707070707 0.7474747474747475;
0.6262626262626263 -1.7171717171717171;
0.42424242424242425 0.5454545454545454;
-1.2727272727272727 0.18181818181818182;
-0.9494949494949495 1.595959595959596;
-0.7878787878787878 -1.2323232323232323;
1.1515151515151516 0.5050505050505051;
0.9494949494949495 1.5151515151515151;
1.7575757575757576 1.7575757575757576;
0.46464646464646464 0.7070707070707071;
1.0303030303030303 2.0;
-0.5858585858585859 0.7070707070707071;
1.1111111111111112 1.0303030303030303;
1.1919191919191918 -0.5050505050505051;
1.5151515151515151 0.7878787878787878;
-0.5050505050505051 -1.8383838383838385;
1.1919191919191918 -0.2222222222222222;
-0.9090909090909091 1.3535353535353536;
1.1111111111111112 0.9090909090909091;
-1.595959595959596 1.595959595959596;
0.7070707070707071 -1.5151515151515151;
1.1111111111111112 1.878787878787879;
-1.0707070707070707 -0.5050505050505051;
-1.2727272727272727 0.020202020202020204;
1.0303030303030303 1.7575757575757576;
-0.3838383838383838 -1.9191919191919191;
0.7474747474747475 -1.5151515151515151;
1.1919191919191918 -0.3838383838383838;
1.3131313131313131 0.9494949494949495;
-1.1919191919191918 -0.1414141414141414;
-0.9494949494949495 1.5555555555555556;
1.1111111111111112 0.7474747474747475;
-0.10101010101010101 0.3434343434343434;
1.0707070707070707 1.9595959595959596;
-1.0707070707070707 1.6767676767676767;
0.42424242424242425 0.5050505050505051;
1.1919191919191918 -0.3838383838383838;
0.8686868686868687 1.5555555555555556;
1.4343434343434343 0.5050505050505051;
-0.6666666666666666 0.9090909090909091;
1.3131313131313131 1.6767676767676767;
1.0303030303030303 1.9595959595959596;
1.3535353535353536 1.9191919191919191;
1.4343434343434343 2.0;
1.1111111111111112 1.0303030303030303;
0.6262626262626263 -1.6767676767676767;
1.0707070707070707 1.5151515151515151;
-0.9494949494949495 -0.9494949494949495;
0.18181818181818182 0.3434343434343434;
-0.5858585858585859 0.7878787878787878;
1.0707070707070707 0.5858585858585859;
1.1919191919191918 0.6262626262626263;
0.6666666666666666 -1.7575757575757576;
1.3535353535353536 1.9595959595959596;
-0.42424242424242425 0.5858585858585859;
1.1919191919191918 -0.3434343434343434;
-1.6363636363636365 1.878787878787879;
-1.1919191919191918 -0.1414141414141414;
-0.9494949494949495 -0.9494949494949495;
-1.1515151515151516 -0.3838383838383838;
1.1919191919191918 1.0707070707070707;
0.98989898989899 1.595959595959596;
1.0707070707070707 0.42424242424242425;
-0.7070707070707071 -1.4747474747474747;
-1.1919191919191918 0.020202020202020204;
0.5454545454545454 0.6666666666666666;
1.797979797979798 1.9191919191919191;
0.3434343434343434 0.5050505050505051;
1.7575757575757576 1.9595959595959596;
1.7575757575757576 1.878787878787879;
-0.8686868686868687 1.2727272727272727;
1.3131313131313131 -0.06060606060606061;
1.1111111111111112 0.5050505050505051;
1.595959595959596 1.1919191919191918;
-0.7878787878787878 1.1111111111111112;
1.1515151515151516 1.6767676767676767;
-1.0707070707070707 1.9595959595959596;
-0.9090909090909091 -0.98989898989899;
-1.1919191919191918 -0.06060606060606061;
0.5050505050505051 0.6666666666666666;
-0.46464646464646464 -1.878787878787879;
-0.5050505050505051 0.6262626262626263;
-0.5858585858585859 -1.6767676767676767;
0.9494949494949495 1.5555555555555556;
1.3535353535353536 1.7575757575757576;
-0.7474747474747475 -1.4343434343434343;
1.1515151515151516 1.9191919191919191;
-0.6262626262626263 -1.5151515151515151;
0.020202020202020204 0.30303030303030304;
0.98989898989899 -0.8686868686868687;
1.0303030303030303 0.30303030303030304;
1.393939393939394 1.2323232323232323;
-1.6767676767676767 1.7171717171717171;
1.3131313131313131 1.2323232323232323;
0.98989898989899 -0.98989898989899;
-1.6363636363636365 1.5151515151515151;
-0.7070707070707071 -1.393939393939394;
1.3535353535353536 0.2222222222222222;
0.7474747474747475 1.0303030303030303;
1.1919191919191918 1.1919191919191918;
0.98989898989899 -1.0707070707070707;
-0.6666666666666666 -1.5555555555555556;
0.5050505050505051 0.7070707070707071;
-0.9090909090909091 -1.0707070707070707;
0.7070707070707071 -1.6767676767676767;
0.9494949494949495 0.30303030303030304;
1.0303030303030303 0.7474747474747475;
1.2727272727272727 0.10101010101010101;
-0.9494949494949495 -1.0303030303030303;
0.8686868686868687 1.5151515151515151;
1.393939393939394 0.3434343434343434;
1.393939393939394 1.9595959595959596;
1.5555555555555556 1.9595959595959596;
1.4747474747474747 1.5151515151515151;
-0.7878787878787878 -1.3535353535353536;
0.6666666666666666 -1.595959595959596;
-1.4343434343434343 0.8282828282828283;
1.0303030303030303 1.7171717171717171;
1.0303030303030303 1.8383838383838385;
0.1414141414141414 0.30303030303030304;
-1.0303030303030303 -0.5454545454545454;
-0.42424242424242425 -1.9595959595959596;
0.9494949494949495 0.3838383838383838;
0.98989898989899 0.46464646464646464;
1.1919191919191918 1.3131313131313131;
-1.0707070707070707 -0.6262626262626263;
1.6363636363636365 1.1919191919191918];
    ds = Systems.henon(zeros(2); a = 1.4, b = 0.3)
    integ_df  = integrator(ds)
    xg = range(-2.,2.,length=100)
    yg = range(-2.,2.,length=100)
    bsn = Basins.basins_map2D(xg, yg, integ_df)

    # Compute boundary
    ind  = findall(iseven.(bsn.basin) .== true)
    basin_test = bsn.basin
    for k in ind; basin_test[k] =basin_test[k]+1; end
    bnd = Basins.get_boundary_filt(basin_test)
    I1 = findall(bnd .== 1);
    v = hcat([[xg[ind[1]]; yg[ind[2]]] for ind in I1 ]...)

    hd = Basins.haussdorff_dist(v,v_bnd');
    @test (hd < 0.3)

    #@test abs(count(bsn_nfo.basin .== 3) - 4127) < 10
    #@test abs(count(bsn_nfo.basin .== -1) - 5730) < 10
end

@testset "Test basin_entropy" begin
    ω=1.; F = 0.2
    ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1)
    integ_df  = integrator(ds; alg=Tsit5(),  reltol=1e-8, save_everystep=false)
    xg = range(-2.2,2.2,length=150)
    yg = range(-2.2,2.2,length=150)
    bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω)
    Sb,Sbb = basin_entropy(bsn.basin; eps_x=20, eps_y=20)
    @test (abs(trunc(Sb;digits=3) - 0.663) < 0.001)
    @test (abs(trunc(Sbb;digits=3) - 0.663) < 0.001)
end

@testset "Test Wada detection" begin
    # Equations of motion:
    function forced_pendulum!(du, u, p, t)
        d = p[1]; F = p[2]; omega = p[3]
        du[1] = u[2]
        du[2] = -d*u[2] - sin(u[1])+ F*cos(omega*t)
    end
    # We have to define a callback to wrap the phase in [-π,π]
    function affect!(integrator)
        if integrator.u[1] < 0
            integrator.u[1] += 2*π
        else
            integrator.u[1] -= 2*π
        end
    end
    condition(u,t,integrator) = (integrator.u[1] < -π  || integrator.u[1] > π)
    cb = DiscreteCallback(condition,affect!)

    F = 1.66;   ω = 1.;    d=0.2
    df = ODEProblem(forced_pendulum!,rand(2),(0.0,20.0), [d, F, ω])
    integ  = init(df, alg=AutoTsit5(Rosenbrock23()); reltol=1e-9, save_everystep=false, callback=cb)
    xg = range(-pi,pi,length=100); yg = range(-2.,4.,length=100)
    bsn = Basins.basins_map2D(xg, yg, integ; T=2*pi/ω)

    # Wada merge Haussdorff distances
    # First remove attractors
    max_dist,min_dist = detect_wada_merge_method(xg,yg,bsn)
    epsilon = xg[2]-xg[1]
    dmax = max_dist/epsilon
    dmin = min_dist/epsilon
    @test (round(dmax) == 6.)
    @test (round(dmin) == 2.)

    W = detect_wada_grid_method(integ, bsn; max_iter=5)
    @test abs(trunc(W[3];digits=2)) ≥ 0.9

end

@testset "Test basin_stability" begin
    ω=1.; F = 0.2
    ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1)
    integ_df  = integrator(ds; alg=Tsit5(),  reltol=1e-8, save_everystep=false)
    xg = range(-2.2,2.2,length=150)
    yg = range(-2.2,2.2,length=150)
    bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω)
    bs = basin_stability(bsn.basin)
    @test abs((trunc(bs[1];digits=3) - 0.492)) < 0.005
    @test abs((trunc(bs[2];digits=3) - 0.507)) < 0.005
end

@testset "Test box_counting_dimension" begin
    ω=1.; F = 0.2
    ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1)
    integ_df  = integrator(ds; alg=Tsit5(),  reltol=1e-8, save_everystep=false)
    xg = range(-2.2,2.2,length=150)
    yg = range(-2.2,2.2,length=150)
    bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω)
    bd = box_counting_dim(xg, yg, bsn)
    @test abs((trunc(bd;digits=2) - 1.9))<0.1
end

@testset "Test compute saddle straddle" begin
    v = [-0.7024272244361562 0.7358943669783878;
    0.6884108286512226 0.2883122496731336;
    0.12248777086662276 0.17266855751576238;
    0.21416001598855067 0.1918853810113003;
    -0.13345893435577877 -0.10181462673868892;
    -0.24612685358870395 0.40609732169904256;
    0.5443525330629784 0.3637084477432412;
    -0.8543699001914774 -0.5504191054878991;
    -0.81424974452034 0.3876573128398122;
    -0.3759850591885575 0.3029466252309591;
    1.2695931242637941 0.5649763379043485;
    0.12319548967885363 0.1988038193210813;
    0.08993268394324054 0.08763876246458675;
    0.5998194677317075 0.42267589676742806;
    -0.7471048204006432 0.4173581434012174;
    -0.3371249654066878 0.41666798624746265;
    0.7345532189737 0.3144499987176077;
    0.24865664622871408 0.28838520561505615;
    -0.7047154958073086 -0.5213541656840109;
    -0.9056118879414881 0.33965597519541424;
    -0.4547496618520956 0.21744361580215943;
    -0.2921283747093624 0.4214569617780614;
    0.6291128784123545 0.3716326220739068;
    -0.6010738837992531 -0.43440268456863923;
    -0.9652688621193279 0.24069456477396253;
    -0.5501113728424021 0.15418201057091677;
    -0.4754568961525241 0.10081275951313767;
    -0.46112021536727277 0.09967906030053533;
    -0.4391298368793252 0.11974265732897998;
    -0.3662749654329448 0.20493692475101746;
    -0.4603329215220823 0.6299851440826713;
    0.41313450130864204 0.28437379172110344;
    -0.7299720705099347 -0.5197447840442747;
    -0.890354901563265 0.3447067895447532;
    -0.44687447282979653 0.22139280314316853;
    -0.2966913365698902 0.45089092769960826;
    0.5501212667337514 0.3597852144857481;
    -0.8200621499512287 -0.5412846282080438;
    -0.8349858140921245 0.3758250755769921;
    -0.39618264893265825 0.2735219498521972;
    -0.9437331293087026 0.7574461204326174;
    0.9293611309012172 -0.15957232746111116;
    1.14599339117127 0.4412985174267645;
    1.2396517289691749 0.5188989178649657;
    0.8959919160009866 0.7090718967680737;
    -1.3461532547680735 0.008447716620129747;
    -0.6398916478847173 0.6224442163197923;
    0.8276869977268546 0.06600919221071469;
    1.0607570462732214 0.6160653104393968;
    -0.20258432736829904 -0.10525350252623214;
    -0.3349774815559638 0.19748108608150314;
    -0.663004994554249 0.724099207141831;
    0.6512663688101122 0.31465908699556383;
    -0.19176348994343262 -0.1154502893224026;
    -0.33183913666489595 0.1997476599045356;
    -0.7339332858876773 0.7450732722276884;
    0.7126484478736383 0.2607963013619859;
    0.3565945556966931 0.376716105526132;
    -1.1424583292686683 -0.5181306982978956;
    -0.6691798553091624 0.48524940958270907;
    -0.7750640627426556 0.7447635013295555;
    0.7617537229889098 0.18145588387180991;
    0.7785640313451194 0.6364791644564665;
    -1.3274874780888364 -0.13525661761177976;
    -0.6135443463613546 0.580073323763926;
    0.8897837152389886 -0.06286426899098416;
    1.1440178842431123 0.4993844802023786;
    1.071146877796164 0.6852925880897688;
    -1.0130674106682127 -0.5202505944779693;
    -0.7194964747274598 0.40984327478863253;
    -0.3337755688232476 0.43117917669428807;
    0.6995015529974392 0.3448638592935367;
    -0.09956822741216353 -0.02900414846930688;
    -0.8926894003760993 0.7683209764135135;
    0.8340505900639127 0.011479529745785512;
    1.0830174110506972 0.5892351694272199;
    0.21204311031642775 0.27421769347596214;
    -0.5644430315888005 -0.4442700978076312;
    -0.9833207542715674 0.22818523879153269;
    -0.5545537845226697 0.1593390081784191;
    -0.472231790321022 0.10490889580215614;
    -0.44780564466696554 0.11274653830575736;
    -0.39238113473288266 0.17168250535822036;
    -0.2884018330542508 0.44862621970732874;
    0.5312147137503425 0.35474079716608226;
    -0.8467395756009193 -0.5497951927484416;
    -0.8188633312213228 0.38589323804591835;
    -0.37959560576794266 0.29752006205944304;
    0.4696497379499575 0.5119045223664362;
    -1.3699491532902666 0.0061234201083442685;
    -0.7483772772691258 0.6806762370674923;
    0.8713539043026773 -0.04321198642007509;
    1.1294232129600181 0.527761159044075;
    0.8977550302213696 0.7036254552586797;
    -1.3389562103162267 -0.04653014856135522;
    -0.6203974969982126 0.6005519054452689;
    0.8447656866536019 0.030093293866878423;
    1.096993725067534 0.5832862071649934;
    0.3206215383741723 0.3693331736195564;
    -1.1173678907969715 -0.5383392012987194;
    -0.6806544022816443 0.47871124803080173;
    -0.626554882949195 0.6976130759887653;
    0.6516324418974114 0.31993125496392566;
    -0.2171788716252682 -0.13785785370964124;
    -0.43298341824677933 0.08012326524258628;
    -0.42802570899573394 0.12623349965496122;
    -0.33932641866474456 0.2438970723609983;
    -1.2365333502089975 -0.2700884633323442];

    v = v';
    ω=1.; F = 0.2
    ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1)
    integ_df  = integrator(ds; alg=AutoTsit5(Rosenbrock23()), reltol=1e-8, abstol=1e-8, save_everystep=false)
    xg = range(-2.2,2.2,length=200)
    yg = range(-2.2,2.2,length=200)
    @time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω)
    sa,sb = compute_saddle(integ_df, bsn, [1], [2]; N=100)
    s = hcat(sa...)
    hd = Basins.haussdorff_dist(s,v)
    @test (hd < 0.4)

end

end


text:
"""
    read_pdos_bin(f::FortranFile)

Read the `pdos_bin` file.
"""
function read_pdos_bin(f::FortranFile)

    fversion = read(f, Float64)
    fheader = strip(convert(String, read(f, FString{80})))

    nk = read(f, IntF)
    ns = read(f, IntF)
    num_popn_orb = read(f, IntF)
    max_eignenv = read(f, IntF)

    # Now start the read the actual data
    species = read(f, (IntF, num_popn_orb))
    ion = read(f, (IntF, num_popn_orb))
    am_channel = read(f, (IntF, num_popn_orb))

    # Allocate storage space
    pdos_weights = zeros(Float64, (num_popn_orb, max_eignenv, nk, ns))
    kpos = zeros(Float64, (3, nk))
    num_eignvalues = zeros(Int, ns)

    for ik in 1:nk
        _, kpos[:, ik] = read(f, IntF, (Float64, 3))
        for is in 1:ns
            read(f, IntF)
            num_eignvalues[is] = read(f, IntF)
            for ib in 1:num_eignvalues[is]
                pdos_weights[:, ib, ik, is] = read(f, (Float64, num_popn_orb))
            end
        end
    end
    return (
        fversion=fversion,
        fheader=fheader,
        nk=nk,
        ns=ns,
        norb=num_popn_orb,
        max_eignenv=max_eignenv,
        species=species,   # Indexing array for species
        ion=ion,   # Indexing array for the ion numbers
        am_channel=am_channel,  # Indexing array for the angular momentum channels (l)
        pdos_weights=pdos_weights,
        kpoints=kpos,
        num_eignvalues=num_eignvalues,
    )
end

"""
    read_pdos_bin(fname::AbstractString)

Read the `pdos_bin` file.
"""
function read_pdos_bin(fname::AbstractString)
    open(fname) do fh
        f = CastepFortranFile(fh)
        read_pdos_bin(f)
    end
end

"""
Reorder the PDOS information

Return a dictionary which can be used as

```julia
output = pdos_index_by_site(pdos_data)
idx = output[1][Orbitals.s]   # Index of the s orbital in the pdos_data.weights
weights = pdos_data.pdos_weights[idx, :, :, :]  # Slice the weight array to obtain the weigths for this site/am combination
```
"""
function pdos_index_by_site(parsed_items)
    mapping = Orbitals.castep_orbital_order
    unique_species = sort(unique(parsed_items.species))
    site_index = 1
    output_data = Dict{Int, Dict{Orbital, Vector{Int}}}()
    for ns in unique_species
        specie_mask = parsed_items.species .== ns
        total_ions = maximum(parsed_items.ion[specie_mask])  # Total number of ions for this species
        # Treat each ion for this speice
        for nion in 1:total_ions
            ion_mask = (parsed_items.ion .== nion) .& specie_mask
            max_am = maximum(parsed_items.am_channel[ion_mask])
            orbs = Dict{Orbital, Vector{Int}}()
            # Angular momentum channel for each site
            for am in 0:max_am
                ion_am_mask = (parsed_items.am_channel .== am) .& ion_mask
                ion_am_idx = findall(ion_am_mask)
                for (iam, iloc) in enumerate(ion_am_idx)
                    #iloc is the index of the orbital
                    #iam is the m channel
                    this_orb = mapping[am+1][((iam - 1) % (2 * am + 1)) + 1]  # m number of this orbital
                    if this_orb in keys(orbs)
                        push!(orbs[this_orb], iloc)
                    else
                        orbs[this_orb] = [iloc]
                end
            end
        end
        output_data[site_index] = orbs
        site_index += 1
        end
    end
    output_data
end

"""
Get an masking for the m channels inferred from the repeated appearance of l channels

CASTEP does not write the m channel labels explicity, but they can be inferred from repeated l channels.
This assumes that l channels of the same atom is written consecutively, which seem to be the case.
"""
function m_channel_mask(pdos_data)
    mapping = Orbitals.castep_orbital_order
    norb = pdos_data.norb
    last_am = pdos_data.am_channel[1]
    current_m = 0
    m_channels = Array{Orbital, 1}(undef, norb)
    m_channels[1] = mapping[last_am + 1][1]
    for i in 2:norb
        current_am = pdos_data.am_channel[i]
        if current_am == last_am
            current_m += 1
        else
            current_m = 0
        end
        m_channels[i] = mapping[current_am+1][(current_m % (2 * current_am + 1)) + 1]
        last_am = current_am
    end
    m_channels
end

precompile(read_pdos_bin, (String,))

export m_channel_mask, pdos_index_by_site