text:
include("basics.jl")
# include("IntSet32.jl")
import DataStructures
import Base: copy, union, in, intersect, parse
import Combinatorics: Combinations
# import Base: complement

##############
#Game Engine
##############

#TODO 1: UX
    #- User enters game state, or chooses a pre-defined one (with a number), or a random pre-defined one ('random')
      #hands are evaluated using eval(parse("..."))
    #- GameState.show emits hands, discard piles, table, contract, etc.
    #- User can input
        #a) card to play (e.g. zt)
        #b) enter for ai to move, n2 for 2 ai moves (so he can control one player),
            #run for full ai mode
        #c) move to trick (goto 0 or restart, goto 2 second trick, goto -1 or 'p' previous trick)
    #At the game end the ai scores the round
    #+1: Az ai kiirja a varhato legrosszabb minimax erteket

#TODO 2: terites - az ai terit, ha minden agon ugyanaz az eredmeny (lejatszastol fuggetlenul)

#IDEA 3: for pruning we can pick key cards,
    #distribute them based on hypergeometric and/or based on existing information
    #and generate hands. Also merge branches with irrelevant differences (drop z9 vs. drop zl)

#a set of cards
#implemented as a memoised set of sets to avoid duplication
#Note: could be implemented as a full-blown class if needed

#TODO: Performance: everything Int and native IntSet32

# const CardSet = IntSet32
# CardSet(cards::Array{Card, 1}) = IntSet32([Int(card) for card in cards])
# const CardSet = DataStructures.IntSet
# CardSet(cards::Array{Card, 1}) = DataStructures.IntSet([Int(card) for card in cards])

const CardSet = CardSet32

hanyAsz(cs::CardSet) = length(intersect(cs, A))
hanyAsz(ca::Vector{Card}) = hanyAsz(CardSet(ca))
hanyTizes(cs::CardSet) = length(intersect(cs, T))
hanyTizes(ca::Vector{Card}) = hanyTizes(CardSet(ca))
hanyAT(cs::CardSet) = length(intersect(cs, AT))
hanyAT(ca::Vector{Card}) = hanyAT(CardSet(ca))
sameSuit(cs::CardSet, suit::Suit) = intersect(cs, suit)
sameSuit(cs1::CardSet, cs2::CardSet, suit::Suit) = union(cs1, cs2) in suit
function sameSuitLarger(cs::CardSet, suit::Suit, card::Card)
#TODO - supporting structure?
end
aduk(cs::CardSet, adu::Suit) = intersect(cs, adu)

#PlayerState
struct PlayerState
    player::Player
    hand::CardSet
    discard::CardSet
    negyven::UInt8 #0 vagy 1
    husz::UInt8 #0-3
end

function PlayerState(player::Player, hand::CardSet)
    PlayerState(player, hand, CardSet(), UInt8(0), UInt8(0))
end

function newHand(ps::PlayerState, hand::CardSet)
    PlayerState(ps.player, hand, ps.discard, ps.negyven, ps.husz)
end

function newDiscard(ps::PlayerState, discard::CardSet)
    PlayerState(ps.player, ps.hand, discard, ps.negyven, ps.husz)
end

function newNegyvenHusz(ps::PlayerState, negyven::UInt8, husz::UInt8)
    PlayerState(ps.player, ps.hand, ps.discard, negyven, husz)
end

#State for a whole gaming session including scores, etc.
mutable struct sessionState
    scores::Vector{Tuple{Player, Int}}
end

#Game State for contract phase
mutable struct LicitState
    contract::Contract  #mi a bemondas

    pakli::Vector{CardSet}       # leosztatlan pakli (egyesevel a sorrend miatt)

    playerStates::Vector{PlayerState}

    currentPlayer::Player
    firstPlayer::Player
    contractHistory::Vector{Contract}
end

function LicitState(numberOfPlayers, firstplayer)
    pakli = toArray(fulldeck)
    playerStates = Vector{PlayerState}()
    for i in 1:numberOfPlayers
        push!(playerStates, PlayerState(Player(i), nocard))
    end
    # Array{Tuple{Player, Contract}}()
    return LicitState(nocontract, pakli, playerStates, firstplayer, firstplayer, Vector{Contract}())
end

#TODO test
function dealInitialCards!(licitState::LicitState)
    @assert (length(licitState.pakli) == 32)
    newPlayerStates = Vector{PlayerState}()
    for pl in licitState.playerStates
        hand = dealCards!(5, licitState.pakli)
        push!(newPlayerStates, newHand(pl, hand))
    end
    licitState.playerStates = newPlayerStates

    return licitState
end

#TODO test
function startGamePhase(licitState::LicitState)
    if length(licitState.contractHistory) < 1
        return nothing
    end

    numberOfPlayers = length(licitState.playerStates)

    #figure out who is playing: felvevo, cotractHistory, first, second
    playersPlaying = reverse([contract.felvevo for contract in licitState.contractHistory])
    playersPlaying = vcat(playersPlaying,
        licitState.firstPlayer,
        nextPlayer(licitState.firstPlayer, numberOfPlayers),
        nextPlayerOffset(licitState.firstPlayer, 2, numberOfPlayers))
    playersPlaying = unique(playersPlaying) #if someone has multiple contracts
    playersPlaying = playersPlaying[1:3] #first three
    #make sure the players are in order with the felvevo being the first
    playersPlaying = sort(playersPlaying, by= x-> (x + numberOfPlayers - licitState.contract.felvevo) % numberOfPlayers)

    #get the nonplayer cards back into pakli and reshuffle
    for playerState in licitState.playerStates
        if !(playerState.player in playersPlaying)
            licitState.pakli = vcat(licitState.pakli, toArray(playerState.hand))
        end
    end

    #2. kor osztas
    newPlayerStates = Vector{PlayerState}()
    for player in playersPlaying
        #7 to felvevo 5 to others playing
        if player == licitState.contract.felvevo
            newCards = dealCards!(7, licitState.pakli)
        else
            newCards = dealCards!(5, licitState.pakli)
        end
        currentCards = licitState.playerStates[player].hand #TODO, BUG - really player? same next line
        push!(newPlayerStates, newHand(licitState.playerStates[player], currentCards + newCards))
    end

    #create GameState
    @assert (length(newPlayerStates) == 3)
    @assert (isempty(licitState.pakli))

    contract = licitState.contract
    playerStates = (newPlayerStates[1], newPlayerStates[2], newPlayerStates[3])
    gs = GameState(contract, playerStates, CardSet(licitState.pakli), nocard, nocard)

    return gs
end

#someone is making a contract
function licit!(licitState::LicitState, lct::Contract)
    if lct != nocontract && lct.totalvalue > licitState.contract.totalvalue
        push!(licitState.contractHistory, lct)
        licitState.contract = lct
    end
    #rotate to next player
    licitState.currentPlayer = nextPlayer(licitState.currentPlayer, length(licitState.playerStates))
    return licitState
end

function dealCards!(n::Integer, pakli::Vector{CardSet})
    cardsDealt = Vector{CardSet}()
    shuffle!(rng, pakli)
    for i in 1:n
        push!(cardsDealt, pop!(pakli))
    end
    return CardSet(cardsDealt)
end

#The game state
#IMPORTANT: as it is copied with shallow copy
#   any field (such as contract) that changes across moves must not be referenced
#   or should be explicitly deep copied by the copy constructor

#TODO: remove pakli? add felvevo (to make ordering more flexible)?
struct GameState
    contract::Contract  #mi a bemondas
    playerStates::Tuple{PlayerState,PlayerState,PlayerState}
    pakli::CardSet       # leosztatlan pakli
    asztal::CardSet  # asztal kozepe - here order matters (pl. 🌰K-ra adu 🎃9 nem ugyanaz mint adu 🎃9-re 🌰K dobas: az egyikre adut kell tenni, a masikra makkot ha van) thus it is not a set
    talon::CardSet      # talon
    currentPlayer::Player
    currentSuit::Suit
    whoseTrick::Player  #o viszi az utest eddig
    lastTrick::Player       #utolso utes
    lastTrick7::Player      #utolso utes adu hetessel
    butLastTrick8::Player   #utolso elotti utes adu 8-assal
    adu7kiment::Player
    adu8kiment::Player
    lastTrickFogottUlti::Player #elbukta a (csendes) ultit
    tricks::UInt8             #number of tricks
    felvevoTricks::UInt8      #number of tricks
    ellenvonalTricks::UInt8   #number of tricks
    felvevoTizesek::UInt8     #tizesek + aszok + utolso utes
    ellenvonalTizesek::UInt8  #tizesek + aszok + utolso utes
    felvevoOsszes::UInt8      #tizesek + aszok + utolso utes + 20 + 40
    ellenvonalOsszes::UInt8   #tizesek + aszok + utolso utes + 20 + 40

    #default constructor
    function GameState(contract::Contract, playerStates::Tuple{PlayerState,PlayerState,PlayerState},
        pakli::CardSet, asztal::CardSet, talon::CardSet,
        currentPlayer::Player=p1, currentSuit::Suit=notrump, whoseTrick::Player=pX,
        lastTrick::Player = pX , lastTrick7::Player = pX,
        butLastTrick8::Player = pX , adu7kiment::Player=pX, adu8kiment::Player=pX,
        lastTrickFogottUlti::Player=pX, tricks::UInt8 = 0x0, felvevoTricks = 0x0, ellenvonalTricks = 0x0,
        felvevoTizesek::UInt8=0x0, ellenvonalTizesek::UInt8=0x0, felvevoOsszes::UInt8=0x0, ellenvonalOsszes::UInt8=0x0)

      new(contract, playerStates, pakli, asztal, talon, currentPlayer, currentSuit, whoseTrick, lastTrick, lastTrick7, butLastTrick8, adu7kiment, adu8kiment, lastTrickFogottUlti, tricks, felvevoTricks, ellenvonalTricks, felvevoTizesek, ellenvonalTizesek, felvevoOsszes, ellenvonalOsszes)
    end

end

#parser
function GameState(bemondas::String, lapok::String)
  #example: ('piros nsz ulti', 't7 t8 p9 pa pf ps pk', 'z8 z9', ... <<<all 8 sets>>>)
  #idea: allow for generic cards and suits, such as zx (any zold) or xx (any card) and deal them randomly
  #checks: are all cards accounted for?
  #nobody has the right number of cards

end

#create from serialised value
function GameState(s::String)
    parseGameState(s)
end

#efficiend storage for game state:
#list of 32 cardsets, 3 bits each -> 96 bits + bemondas a #tablazatbol
# function compact()
# function decompact()

#TODO
# #create signature
# function serialize(stream, g::GameState)
# end

#copy constructor
# function copy(g::GameState)
#     ps2 = (copy(g.playerStates[1]), copy(g.playerStates[2]), copy(g.playerStates[3]))

#     GameState(copy(g.contract), copy(g.pakli), copy(g.asztal),
#         copy(g.talon), ps2, g.currentPlayer,
#         g.currentSuit, g.lastTrick, g.lastTrick7, g.butLastTrick8,
#         g.adu7kiment, g.adu8kiment, g.lastTrickFogottUlti, g.tricks,
#         g.felvevoTricks, g.ellenvonalTricks, g.felvevoTizesek,
#         g.ellenvonalTizesek, g.felvevoOsszes, g.ellenvonalOsszes)
# end

#parse a list of cards
#Format: known cards (e.g. p9) are marked suit (non-capital t,z,m,p) and face (capital 7,8,9,U,F,K,A)
#Special suits: x: any suit, tz: tok vagy zold, pmz: piros vagy makk vagy zold
#Special faces: X:any face Z:'nem fontos' [9UFK] .:Y [789UFK], explicit: pl. [AT] az asz v tizes
#Examples: "pATYY mAKX z.. tK[U987][U987]":
#Tricky case: a kicsi t tokot, a nagy T Tizest jelol a felreertesek elkerulese vegett
#Alternativ jelolesek:  also: UVJ, tizes: T vagy 10, tok: t vagy o
#Players are separated by '|' or '/'
#Contract and cards are separated by ';'
function parseCards(s::AbstractString)

    #Ideas: suit then cards x, a(du)bcd, [80%]bFK [20%E] eros(AT), xN10 (fill until 10 cards)
    #first match abcd then single cards then probabilities greedily.
    #Q: try again if fails? how many times? maybe give warning? Try ALL?
    #bemondas: regexp for p[EHR].x[Ult|Rep|40|4A...]K*  (vagy RK, SK MK)

    regexp = r"(?<suit>(t|z|m|p|n|x|a|b|c|d)+)(?<face>(A|T|K|F|U|9|8|7|X|Y|X|N[0-9]+)+)\s*"

    # matches = [mtch for mtch in eachmatch(regexp, s)]
    cs = Vector{CardSet}()
    for mtch in eachmatch(regexp, s)
        suits = CardSet32()
        for suit in mtch[:suit]
            suits += suitStrings[suit]
        end

        faces = CardSet32()
        for face in mtch[:face]
            faces += faceStrings[face]
        end

        push!(cs, suits * faces)
    end

    #TODO: N cards, [80%], X,Y,Z, a,b,c,d,etc.
    asSet = union(CardSet(), CardSet(), cs...) #extra empty sets if cs is empty or has one element
    return (cs, asSet)
end

function parse(CardSet, s)
    parseCards(s)[2]
end

function parseManyCards(s::AbstractString)
        return [parseCards(cards)[2] for cards in split(s, "|")]
end

function parsePlayerStates(s)
    paramStrings = split(s, "|")
    #from typeof(additionalParams)
    types = [UInt8, CardSet32, CardSet32, UInt8, UInt8]
    pStates = []
    for playerString in Iterators.partition(paramStrings, length(types))
        params = []
        for (paramString, ptype) in zip(playerString, types)
            #exceptions are propagated
            param = parse(ptype, paramString)
            push!(params, param)
        end
        push!(pStates, PlayerState(params...))
    end

    #return as tuple
    return (pStates[1], pStates[2], pStates[3])
end

function parseParams(s)
    paramStrings = split(s, "|")
    #from typeof(additionalParams)
    params = []
    types = [CardSet32, CardSet32, CardSet32, UInt8,CardSet32,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8,UInt8]
    for (paramString, ptype) in zip(paramStrings, types)
        #exceptions are propagated
        param = parse(ptype, paramString)
        push!(params, param)
    end
    #give them names
    # (asztal, talon, currentPlayer, currentSuit, whoseTrick, lastTrick, lastTrick7, butLastTrick8, adu7kiment, adu8kiment, lastTrickFogottUlti, tricks, felvevoTricks, ellenvonalTricks, felvevoTizesek, ellenvonalTizesek, felvevoOsszes, ellenvonalOsszes) = params
    return params
end

function parseGameState(s::String)
    elements = split(s, ';')
    if length(elements) > 2
        #full canonical state
        contract = parseContract(elements[1])
        playerStates = parsePlayerStates(elements[2])
        additionalParams = parseParams(elements[3])

        #create state
        gs = GameState(contract, playerStates, additionalParams...)
    else #TODO
        #heuristic parsing - hands, talon defined, rest is pakli
        if(length(elements)) == 1
            contract = nocontract
        else
            contract = parseContract(elements[1])
        end

        cardSets = parseManyCards(elements[end]) #last
        cards = matchCards(cardSets, fulldeck)
        playerStates = (
            PlayerState(Player(1), cards[1]),
            PlayerState(Player(2), cards[2]),
            PlayerState(Player(3), cards[3]),
            )
        if length(cardSets) == 3
            gs = GameState(contract, playerStates, CardSet(), CardSet(), CardSet())
        elseif length(cardSets) == 4 #hands + talon
            gs = GameState(contract, playerStates, CardSet(), CardSet(), cards[4])
        else #TODO 4-5-6 players?
            throw("invalid state")
        end
    end

    return gs
end

#match cards to a pattern
#Example: "pAKxx zT7x mKF tx" , full deck
#Eliminate by rows (every spot must have one card)
#and columns (one card can only be maximum at one place at the end)
#if no obvious step, match a card randomly
#TODO test
# old commment: deal fix cards, then coloured cards (px), then multicouloured cards (pz.) then the rest (xx, xA)
function matchCards(cardSlots::Vector{CardSet}, cards::CardSet)
    #build 2D binary (logical) matrix
    #see Sodoku solver at: http://www.math.cornell.edu/~mec/Summer2009/meerkamp/Site/Solving_any_Sudoku_II.html

    isDone = [false for cs in cardSlots] #mark elements committed

    (slotsByWeight, cards) = commitAllUnique(cardSlots, isDone, cards) #do twice at the beginning to make sure all simple ones are resolved before the complicated searches start

    while !all(isDone)

        #for (cs, dn) in zip(cardSlots, isDone); (print(dn), print(cs), println()) end; (print(cards), println()); for (k,v) in slotsByWeight print(" $(k) $(v)") end #DEBUG

        (slotsByWeight, cards) = commitAllUnique(cardSlots, isDone, cards)

        if all(isDone) break end

        #if all remaining slots are the same, distribute -> done
        if length(slotsByWeight) == 1
            for (n, indices) in slotsByWeight
                @assert (length(indices) > 0)
                firstElem = cardSlots[first(indices)]
                if all(x -> cardSlots[x] == firstElem, indices) #all equal
                    if(length(firstElem * cards) < length(indices))
                        throw("not enough cards to distribute indices: $(indices) cards: $(cards)")
                    end

                    #distribute randomly
                    cardsToCommit = shuffle(rng, toArray(firstElem * cards))
                    cards = commitMany!(cardSlots, isDone, indices, cardsToCommit, cards)
                end
            end
            @assert (all(isDone))
        end
        #for (cs, dn) in zip(cardSlots, isDone); (print(dn), print(cs), println()) end; (print(cards), println()); for (k,v) in slotsByWeight print(" $(k) $(v)") end #DEBUG

        if all(isDone) break end

        #look for preemptive sets (for example 3 slots with pATK, pATK, pATK that we can commit)
        foundPreemptiveSet = false
        for (n, indices) in slotsByWeight
            if length(indices) < n
                continue
            end #for length x we need at least x in the set

            (firstPreemptiveSet, cardsInSet) = findPreemptiveSet(n, indices, cardSlots) #if we find one we commit immediately

            if !isempty(firstPreemptiveSet)
                @assert (length(firstPreemptiveSet) == n && length(cardsInSet) == n) #exactly n cards needed
                foundPreemptiveSet = true
                #
                #TODO: we just randomly do that
                #we could generate all permutations if needed using Combinatorics.Permutations(toArray(cards), n)
                cardsToCommit = shuffle(rng, toArray(cards))
                cards = commitMany!(cardSlots, isDone, firstPreemptiveSet, cardsToCommit, cards)
                break
            end
        end

        #for (cs, dn) in zip(cardSlots, isDone); (print(dn), print(cs), println()) end; (print(cards), println()); for (k,v) in slotsByWeight print(" $(k) $(v)") end #DEBUG

        #all hard criteria fullfilled - randomly fill the rest
        if !foundPreemptiveSet
            #TODO
        end
    end

    return (cardSlots, cards)
end

# deal the rest of the cards based on current cards already dealt
function sampleCards(cards::CardSet)
    cards = CardSet([pF,p9,p8, mT,mF])
    restOfCards = fulldeck - cards
    deck_to_deal = collect([x for x in restOfCards])
    deck_to_deal = shuffle(deck_to_deal)
    my_cards = cards + CardSet32(deck_to_deal[1:7])
end

#find single cards and commit them
@inline function commitAllUnique(cardSlots, isDone, cards)
    slotsByWeight = Dict{Int, Vector{Int}}()
    for i in 1:length(cardSlots)
        if isDone[i] #already done
            continue
        end

        cardSlots[i] *= cards #we cannot wish for cards not in the deck anymore
        cs = cardSlots[i]
        len = length(cs)
        if 0 == len #commit
            throw("no (matching) card - cannot fulfill criteria. index: $(i), card: $(cs), deck: $(cards)")
        elseif 1 == len #commit
            cards = commit!(isDone, i, cs, cards)
        else
             slotsByWeight[len] = push!(get(slotsByWeight, len, Vector{Int}()), i)
        end
    end

    return (slotsByWeight, cards)
end

#find the first preemptive set
function findPreemptiveSet(n, indices, cardSlots)
    for comb in combinations(indices, n) #try all combinations
        allCardsInComb = CardSet(cardSlots[comb]) #array constructor -> union
        if length(allCardsInComb) == n #n cards in n slots
            return (comb, allCardsInComb)
        end
    end

    return ([], CardSet()) # no preemptive set
end

function commit!(isDone, i, card, cards)
    if !(card in cards) throw("card not in cards - cannot fulfill criteria") end
    cards = cards - card
    isDone[i] = true

    return cards
end

function commitMany!(cardSlots, isDone, indices, cardsToCommit, cards)
    for (i, card) in zip(indices, cardsToCommit)
        @assert (card in cardSlots[i])
        cardSlots[i] *= card

        cards = commit!(isDone, i, card, cards)
    end

    return cards
end

function ps(g::GameState, player::Player)
    return g.playerStates[player]
end

#Evaluate the score based on the game state
#if the game is over, it returns the final score
#otherwise the sum of contracts that have been decided
#notrump contracts yield 0
#TODO instead of Int, return (10,-5,-5) for proper handling of csendesUlti and csendesDuri
function score(g::GameState, ce::ContractElement)
    #TODO: maybe check after full tricks only?
    if ce.bem == parti
        if  g.felvevoOsszes >=100 && g.felvevoTizesek > 0 #ha szaz van nyertunk de egy tizes mindig kell
            return 2 * ce.val #szazzal van meg
        end
        if g.ellenvonalOsszes >=100 && g.ellenvonalTizesek > 0
            return -2 * ce.val #szazzal bukott
        end
        if g.tricks == 10 #vege a jateknak
            @assert (g.felvevoTizesek + g.ellenvonalTizesek == 90)
            if g.felvevoOsszes > g.ellenvonalOsszes
                return ce.val
            else
                return -ce.val
            end
        end
        kiadoUtesek = 90 - (g.felvevoTizesek + g.ellenvonalTizesek)
        if g.felvevoOsszes > kiadoUtesek + g.ellenvonalOsszes #Uint(0) - Uint(20) = 236!
            return ce.val #mar barhogyan megnyertuk, akkor is ha a kiado utesek az ellenhez kerulnek
        end
        if g.ellenvonalOsszes > kiadoUtesek + g.felvevoOsszes
            return -ce.val #mar barhogyan megnyertek, akkor is ha a kiado utesek hozzam kerulnek
        end
        return 0

    elseif ce.bem == ulti
        if g.lastTrick7 == p1 return ce.val end #megvan
        if g.tricks == 10 || (g.tricks < 9 && g.adu7kiment!=pX) || (g.tricks == 9 && g.adu7kiment!=pX && !in(Card(g.contract.suit, _7), g.asztal))#TODO fix
            if ce.modosito == elolrol return -ce.val end #elolrol szimplan bukik
            if isempty(ce.kon) return -2 * ce.val end #hatulrol duplan bukik #TODO kiveve ha nem volt ott a hetes
            return -3 * div(ce.val, 2) #hatulrol kotraval triplan bukik, rekontraval 6x, stb.
        end
        return 0    #meg nem tudjuk

    elseif ce.bem == repulo
        if g.butLastTrick8 == p1 return ce.val end #megvan
        if g.tricks >= 9 || (g.tricks < 8 && g.adu8kiment!=pX) return -ce.val end #bukott
        return 0    #meg nem tudjuk

    elseif ce.bem == negyvenSzaz
        if g.felvevoTizesek >= 60 && g.playerStates[1].negyven > 0 #megvan a szaz
            return ce.val
        end
        if g.ellenvonalTizesek >= 40 || g.playerStates[1].negyven == 0 #bukott vagy sose volt
            return -ce.val
        end
        return 0

    elseif ce.bem == huszSzaz
        if g.felvevoTizesek >= 80 && g.playerStates[1].husz > 0 #megvan a szaz
            return ce.val
        end
        if g.ellenvonalTizesek >= 20 || g.playerStates[1].husz == 0 #bukott vagy sose volt
            return -ce.val
        end
        return 0

    elseif ce.bem == negyAsz
        if hanyAsz(g.playerStates[1].discard) == 4 #negy aszt vitt
            return ce.val
        end
        ellenvonalUtesek = union(g.playerStates[2].discard, g.playerStates[3].discard, g.talon)
        if hanyAsz(ellenvonalUtesek) > 0 #legalabb egy asz az ellenvonal utesei kozott
            return -ce.val
        end
        return 0

    elseif ce.bem == negyTizes
        if hanyTizes(g.playerStates[1].discard) == 4 #negy tizest vitt
            return ce.val
        end
        ellenvonalUtesek = union(g.playerStates[2].discard, g.playerStates[3].discard, g.talon)
        if hanyTizes(ellenvonalUtesek) > 0 #legalabb egy tizes az ellenvonal utesei kozott
            return -ce.val
        end
        return 0

    elseif ce.bem == durchmars || ce.bem == redurchmars
        if length(g.playerStates[2].discard) + length(g.playerStates[3].discard) > 0
            return -ce.val #ellenvonal utes -> bukott
        end
        if g.tricks == 10 #megvan
            @assert (length(g.playerStates[1].discard) == 30) #TODO: FELVEVO maybe or count both tricks
            return ce.val #megvan
        end
        return 0    #meg nem tudjuk

    elseif ce.bem == betli || ce.bem == rebetli
        if length(g.playerStates[1].discard) > 0
            return -ce.val #felvevo utes -> bukott
        end
        if g.tricks == 10 #megvan
            @assert (length(g.playerStates[2].discard) + length(g.playerStates[3].discard) == 30)
            return ce.val #megvan
        end
        return 0    #meg nem tudjuk

    elseif ce.bem == csendesUlti
        if g.lastTrick7 == p1 return ce.val end #megvan
        if g.lastTrick7 == p2 return -div(ce.val, 2) end #csak 1 jatekosnak fizetem, ezert a fele
        #TODO: WRONG! fix for p2, p3 to optimize their play
        if g.lastTrick7 == p3 return -div(ce.val, 2) end

        if g.lastTrickFogottUlti == g.playerStates[1].player #bukott #TODO - new state for that
            return -2 * ce.val #hatulrol duplan bukik
        end
        return 0    #meg nem tudjuk vagy nincs

    elseif ce.bem == csendesDuri
        if g.tricks == 10
            if length(g.playerStates[1].discard) == 30 #megvan
                return ce.val
            end #megvan
            if length(g.playerStates[1].discard) == 27 && g.lastTrick != p1
                return -ce.val #bukott, ha csak az utolso utest bukja #TODO fogott ellenvonal duri
            end
        end
        return 0    #meg nem tudjuk vagy nincs
    else
        #TODO: raise exception here
    end
end

#Evaluate the score based on the game state
function score(g::GameState)
    totalScore = 0
    for contractElement in g.contract.contracts
        totalScore += score(g, contractElement)
    end
    return totalScore
end
##############
#UX/UI
##############

function flushScreen()
    println("\n"^100)
end

function printPlayerState(io::IO, pstate, shortFormat::Bool=true, licitPhase::Bool=false)
    if shortFormat
        for pl in pstate

        end
        for pl in pstate
            print(io, pl.player); print(io, "|")
            print(io, pl.hand, true); print(io, "|")
            print(io, pl.discard, true); print(io, "|")
            print(io, pl.negyven); print(io, "|")
            print(io, pl.husz); print(io, "|")
        end
    else
        for pl in pstate
            print(io, playerNames[pl.player]); print(io, ": ")
            print(io, pl.hand, false)
            println(io); println(io)
        end
        if !licitPhase
            println(io)
            for pl in pstate
                print(io, playerNames[pl.player]); print(io, " ütései: ")
                print(io, pl.discard, false)
                println(io)
            end
            for pl in pstate
                if pl.negyven + pl.husz > 0
                    print(io, playerNames[pl.player])
                    if pl.negyven > 0
                        print(io, " negyven:"); print(io, pl.negyven)
                    end
                    if pl.husz > 0
                        print(io, " husz:"); print(io, pl.husz)
                    end
                    println(io)
                end
            end
        end
    end
end

function print(io::IO, pstate::Tuple{PlayerState, PlayerState, PlayerState}, shortFormat::Bool=true)
    printPlayerState(io, pstate, shortFormat, false)
end
function print(io::IO, pstate::Vector{PlayerState}, shortFormat::Bool=true)
    printPlayerState(io, pstate, shortFormat, true)
end

display(pstate::Tuple{PlayerState, PlayerState, PlayerState}) = print(stdout, pstate, false)

function print(io::IO, g::GameState, shortFormat::Bool=true)
    if shortFormat
        print(io, "G*")
        print(io, g.contract, true); print(io, ";")
        print(io, g.playerStates, true); print(io, ";")
        print(io, g.pakli, true); print(io, "|")
        print(io, g.asztal, true); print(io, "|")
        print(io, g.talon, true); print(io, "|")
        additionalParams = [g.currentPlayer, g.currentSuit, g.whoseTrick, g.lastTrick, g.lastTrick7, g.butLastTrick8, g.adu7kiment, g.adu8kiment, g.lastTrickFogottUlti, g.tricks, g.felvevoTricks, g.ellenvonalTricks, g.felvevoTizesek, g.ellenvonalTizesek, g.felvevoOsszes, g.ellenvonalOsszes]
        for param in additionalParams
            # print(Int(param))
            print(io, param)
            print(io, "|")
        end

        #TODO: additional parameters separated by "|"
    else
        println(io); print("Bemondás: "); print(io, g.contract, false); println(io)
        print(io, "\nAsztal: "); print(io, g.asztal, false); println(io); println(io); println(io)
        print(io, g.playerStates, false)
        print("Talon: "); print(io, g.talon, false); println(io)

        print(io, "\n$(playerNames[g.currentPlayer]) jön. Lehetséges hívasok: "); print(io, validMoves(g), false)
        print(io, "\n\nÜtések száma: $(g.tricks)")
        print(io, "\nNyeremény: $(score(g))")
        # print(io, "\nadu: $(g.contract.suit)")
        print(io, "\nFelvevő    Tízesek:$(g.felvevoTizesek)    Összes:$(g.felvevoOsszes)    Ütések:$(g.felvevoTricks)")
        print(io, "\nEllenvonal Tízesek:$(g.ellenvonalTizesek)    Összes:$(g.ellenvonalOsszes)    Ütések:$(g.ellenvonalTricks)")
        print(io, "\nUtolsó ütés:$(playerNames[g.lastTrick])    adu 7 utolsó:$(playerNames[g.lastTrick7])    adu 8 utolsó előtti:$(playerNames[g.butLastTrick8])    csendes ulti bukott:$(playerNames[g.lastTrickFogottUlti]) \nadu 7 kiment:$(playerNames[g.adu7kiment]) adu 8 kiment:$(playerNames[g.adu8kiment])")
        print(io, "\n\nID: \""); print(io, g, true); print(io, "\"")
    end
end

display(g::GameState) = print(stdout, g, false)

#TODO
function print(io::IO, ls::LicitState, shortFormat::Bool=true)
    if shortFormat
        print(io, "L*")
        print(io, ls.contract, true); print(io, ";")
        print(io, ls.playerStates, true); print(io, ";")
        print(io, CardSet(ls.pakli), true); print(io, "|")
        print(io, ls.firstPlayer); print(io, "|")
        print(io, ls.currentPlayer); print(io, "|")
        for bem in ls.contractHistory
            print(io, bem, true); print(io, "|")
        end
    else
        print(io, ls.playerStates, false)
        print(io, "$(playerNames[ls.currentPlayer]) jön.");  println(io)
        println(); print("Bemondás: "); print(io, ls.contract, false)
        print(io, "\nKezdő játékos: $(playerNames[ls.firstPlayer])")
        print(io, "\nEddigi bemondások:")
        for bem in ls.contractHistory
            println();print(io, bem, false);
        end
        println(io)
        print(io, "\nID:\""); print(io, ls, true); print(io, "\"")
        println(io)
    end
end

display(ls::LicitState) = print(stdout, ls, false)

#The valid card to play in this trick (if there is already a card on the table)
function validMoves(g)
    # valid = copy(ps(g, g.currentPlayer).hand)
    valid = ps(g, g.currentPlayer).hand #copy not necessary if on stack

    if length(valid) <= 1 return valid end

    #sorrend: felulutni szinbol, szin, felulutni adu, adu, egyeb
    if length(g.asztal) > 0 #mar van egy szin
        @assert (g.currentSuit != notrump)
        enSuit = sameSuit(valid, g.currentSuit)
        if length(enSuit) > 0
            if g.currentSuit != g.contract.suit && !isempty(aduk(g.asztal, g.contract.suit)) #ha mar aduval utottek nem kell felul utni
                valid = enSuit
            else
                largestEnSuite = largestCard(g.asztal, g.currentSuit) #best in suite - ezt kell felulutni
                enSuitLarger = whichTrumps(enSuit, largestEnSuite, g.contract.suit)
                valid = length(enSuitLarger) > 0 ? enSuitLarger : enSuit
            end
        elseif g.contract.suit <= p #szines jatek -> adu?
            myTrumps = intersect(g.contract.suit, valid) #ha nincs szin adu
            if length(myTrumps) > 0 #nincs szin de van adu
                largest = largestCard(g.asztal, g.contract.suit)
                valid = myTrumps
                if SuitFace(largest)[1] == g.contract.suit #mar van adu felul kell utni
                    myTrumpsLarger = whichTrumps(myTrumps, largest, g.contract.suit) #ha van mar adu azt is felul kell utni
                    valid = length(myTrumpsLarger) > 0 ? myTrumpsLarger : myTrumps
                end
            end
        end
    end

    #adu 7-8 ulti rep nem lehet ido elott

    #ulti-repulonel valaszthatunk
    if length(valid) == 2 && g.contract.suit <= p &&
        in(Card(g.contract.suit, _7), valid) && in(Card(g.contract.suit, _8), valid) &&
        isUlti(g.contract) && isRepulo(g.contract)
        #NOP
    elseif length(valid) > 1 && isUlti(g.contract)
        valid = setdiff(valid, Card(g.contract.suit, _7))
    end
    if length(valid) > 1 && isRepulo(g.contract)
        adu8 = Card(g.contract.suit, _8)
        if g.tricks < 9
            valid -= adu8 #nem szabad kijatszani
        elseif g.tricks == 9 && adu8 in valid
            valid = adu8 #kotelezo kijatszani
        end
    end

    return valid
end

function negyvenHusz(ps::PlayerState, trump::Suit)
    negyven = husz = 0
    for suit in [t,z,m,p]
        if issubset(union(Card(suit, F), Card(suit, K)), ps.hand)
            if suit == trump negyven += 1
            else husz +=1
            end
        end
    end
end

function newState(gOld::GameState, contract::Contract)
    throw("implement me!")
end

function newState(gOld::GameState, negyven::UInt8, husz::UInt8)
    if g.tricks > 1 throw(ArgumentError("negyven, husz csak a jatek elejen")) end
    if g.contract.suit > p throw(ArgumentError("szines jatek kell")) end
    if negyven > 1 || husz > 3 || husz < 0 || negyven < 0 throw(ArgumentError()) end
    (negy, hu) = negyvenHusz(ps(g, g.currentPlayer), g.contract.suit)
    if negy < negyven || hu < husz
        throw(ArgumentError("nincs negyven/husz"))
    end

    g = GameState(gOld)
    ps(g, g.currentPlayer).negyven += negyven
    ps(g, g.currentPlayer).husz += husz

    return g
end

function newState(g::GameState, card::Card)
    @assert (length(card) == 1)

    #copy them one by one, change and reconstruct - poor man's update for immutables
    contract=g.contract; pakli=g.pakli; asztal=g.asztal; talon=g.talon;
    currentPlayer=g.currentPlayer; currentSuit=g.currentSuit; whoseTrick=g.whoseTrick;
    lastTrick=g.lastTrick ; lastTrick7=g.lastTrick7;
    butLastTrick8=g.butLastTrick8 ; adu7kiment=g.adu7kiment; adu8kiment=g.adu8kiment;
    lastTrickFogottUlti=g.lastTrickFogottUlti; tricks=g.tricks;
    felvevoTricks=g.felvevoTricks; ellenvonalTricks=g.ellenvonalTricks;
    felvevoTizesek=g.felvevoTizesek; ellenvonalTizesek=g.ellenvonalTizesek;
    felvevoOsszes=g.felvevoOsszes; ellenvonalOsszes=g.ellenvonalOsszes

    playerStates = [g.playerStates[1], g.playerStates[2], g.playerStates[3]]
    # g = GameState(gOld) #TODO make sure consistency - deep copy for immutability, or custom copy constructor for speed?

    if length(talon) < 2
        talon = add(talon, card)
        hand = remove(playerStates[currentPlayer].hand, card)
        playerStates[currentPlayer] = newHand(playerStates[currentPlayer], hand)
        if hanyAT(card) > 0 # talon tizese az ellenfele
            ellenvonalTizesek += UInt8(10)
            ellenvonalOsszes+= UInt8(10)
        end
    else
        asztal = add(asztal, card)
        hand = remove(playerStates[currentPlayer].hand, card)
        playerStates[currentPlayer] = newHand(playerStates[currentPlayer], hand)

        if card == Card(contract.suit, _7) adu7kiment = currentPlayer end
        if card == Card(contract.suit, _8) adu8kiment = currentPlayer end

        if length(asztal) == 1 #uj szin
            currentSuit = SuitFace(card)[1]
            whoseTrick = currentPlayer
        else
            if trumpsAll(card, asztal, contract.suit)
                whoseTrick = currentPlayer
            end
        end

        if length(asztal) == 3 #kor vege
            tizasz = hanyAT(asztal)

            #update variables
            tricks += UInt8(1)
            if tricks == UInt8(10)
                lastTrick = whoseTrick

                if lastTrick == p1 #utolso utes = 10
                    felvevoTizesek += UInt8(10)
                    felvevoOsszes += UInt8(10)
                else
                    ellenvonalTizesek += UInt8(10)
                    ellenvonalOsszes+= UInt8(10)
                end

                if contract.suit <= p
                    trump7 = Card(contract.suit, _7)
                    if trump7 == largestCard(asztal, contract.suit)
                        lastTrick7 = whoseTrick
                    else
                        lastTrickFogottUlti = adu7kiment #(csendes) ulti bukott
                    end
                end
            end

            #repulo
            if tricks == 9 && contract.suit <= p
                trump8 = Card(contract.suit, _8)
                if trump8 == largestCard(asztal, contract.suit)
                    butLastTrick8 = whoseTrick
                end
            end

            if whoseTrick == p1
                felvevoTricks += UInt8(1)
                felvevoTizesek += UInt8(tizasz * 10)
                felvevoOsszes+= UInt8(tizasz * 10)
            else
                ellenvonalTricks += UInt8(1)
                ellenvonalTizesek += UInt8(tizasz * 10)
                ellenvonalOsszes+= UInt8(tizasz * 10)
            end

            discard = add(playerStates[whoseTrick].discard, asztal)
            playerStates[whoseTrick] = newDiscard(playerStates[whoseTrick], discard)
            asztal = CardSet()

            currentSuit = notrump
            currentPlayer = whoseTrick
            whoseTrick = pX

        else #no trick yet
            currentPlayer = nextPlayer(g.currentPlayer)
        end
    end

    return GameState(contract, (playerStates[1],playerStates[2],playerStates[3]), pakli, asztal, talon, currentPlayer, currentSuit, whoseTrick, lastTrick, lastTrick7, butLastTrick8, adu7kiment, adu8kiment, lastTrickFogottUlti, tricks, felvevoTricks, ellenvonalTricks, felvevoTizesek, ellenvonalTizesek, felvevoOsszes, ellenvonalOsszes)

end

##############
#Tests
##############

function engineTest()
    matchCards([tA, tK, tF+tA+tU+t8, tU, tF+tA+tU+t8, tF+tA+tU+t8+t7], t)
end

function parseTest()
    parseManyCards("pApTpKpUp7mKmUm8z9tT")
    parseManyCards("pATKU7 mKU8 z9 tT")
    # parseManyCards("pApTpKpUp7mKmUm8z9tT;pFp9p8mAmFm7zKzFz7tF;mTm9zAzTzUz8tAtKt8t7")
    parseManyCards("pApTpKpUp7mKmUm8z9tT|pFp9p8mAmFm7zKzFz7tF|mTm9zAzTzUz8tAtKt8t7")
    parseManyCards("pATKU7 mKU8 z9 tT | pF98 mAF7 zKF7 tF | mT9 zATU8 tAK87")
    # parseManyCards("❤ A ❤ T ❤ K ❤ U ❤ 7 🌰 K 🌰 U 🌰 8 🍃 9 🎃 T | ❤ F ❤ 9 ❤ 8 🌰 A 🌰 F 🌰 7 🍃 K 🍃 F 🍃 7 🎃 F |  T 🌰 9 🍃 A 🍃 T 🍃 U 🍃 8 🎃 A 🎃 K 🎃 8 🎃 7", )
    parseManyCards("pAKXX mTXX tX zA7N")
    parseManyCards("xXXXXXXXXXX")
    parseManyCards("xXXXXXXXXXXXX|xXXXXXXXXXX|xXXXXXXXXXX")
    parseManyCards("aAXX7 bATXX cXX")
    parseManyCards("aAXX7 bATXX cXX dXX")
    parseManyCards("[80%]aXXXX [10%]aXXXXX [5%]aXXX [5%]aXXXXXX xN10")
    #hibas

    parseContract("pEUltKE")
    parseContract("pRep")
    parseContract("Ep") #piros parti
    parseContract("EpEK") #elolrol kontra parti
    parseContract("EpHK") #elolrol kontra parti
    parseContract("EpEKHRK") #elolrol kontra hatulrol rekontra
    parseContract("pUltRep4A") #parti
    parseContract("p EP RRep RCsu") #parti
    parseContract("4T") #szintelen 4T
    parseContract("H Bet") #szintelen betli
    parseContract("H pTbet")
end

max_stars_repo_name: albapa/ulti

max_issues_repo_name: albapa/ulti

max_forks_repo_name: albapa/ulti


text:
"""
`DataValueArray{T, N}` is an efficient alternative to `Array{DataValue{T}, N}`.
"""
struct DataValueArray{T,N} <: AbstractArray{DataValue{T},N}
    values::Array{T,N}
    isna::Array{Bool,N}

    function DataValueArray{T,N}(d::AbstractArray{T, N}, m::AbstractArray{Bool, N}) where {T,N}
        if size(d) != size(m)
            msg = "values and missingness arrays must be the same size"
            throw(ArgumentError(msg))
        end
        new(d, m)
    end
end

const DataValueVector{T} = DataValueArray{T, 1}
const DataValueMatrix{T} = DataValueArray{T, 2}

ext: jl


text:
"""
```
measurement(m::KrusellSmithCT{T}, TTT::Matrix{T},
                              RRR::Matrix{T}, CCC::Vector{T})
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::KrusellSmithCT{T},
                     TTT::Matrix{T},
                     RRR::Matrix{T},
                     CCC::Vector{T},
                     inverse_basis::SparseMatrixCSC{T,Int64}) where {T<:AbstractFloat}
    return measurement(m, TTT, RRR, CCC, full(inverse_basis))
end

function measurement(m::KrusellSmithCT{T},
                     TTT::Matrix{T},
                     RRR::Matrix{T},
                     CCC::Vector{T},
                     inverse_basis::Matrix{T}) where {T<:AbstractFloat}

    endo = m.endogenous_states

    track_lag = get_setting(m, :track_lag)
    freq = get_setting(m, :state_simulation_freq)

    QQ = eye(freq) * 1/freq # freq = no. of partitions of quarter
    EE = zeros(1,1)
    DD = [0.]
    if track_lag
        ZZ = zeros(1, n_states(m))
        ZZ[endo[:output][1]] = 1
        augmented_inverse_basis = zeros(n_states(m), (freq+1)*size(TTT, 1))
        augmented_inverse_basis[:,end - size(inverse_basis,2)+1:end] = inverse_basis
        ZZ = ZZ*augmented_inverse_basis
    else
        ZZ = zeros(1, n_states(m))
        ZZ[endo[:output][1]] = 1
        augmented_inverse_basis = zeros(n_states(m), freq*size(TTT, 1))
        augmented_inverse_basis[:,end - size(inverse_basis,2)+1:end] = inverse_basis
        ZZ = ZZ*augmented_inverse_basis
    end
    return Measurement(ZZ, DD, QQ, EE)

end

function measurement(m::KrusellSmithCT{T},
                     TT::Matrix{T},
                     RR::Matrix{T},
                     CC::Vector{T},
                     inverse_basis::Matrix{T},
                     dt::T,
                     track_lag::Bool = true;
                     abs_err::Float64 = 1e-8,
                     max_iter::Int64 = size(TT, 1)*2) where {T<:AbstractFloat}

    endo = m.endogenous_states

    freq = get_setting(m, :state_simulation_freq)
    track_last_state = get_setting(m, :track_lag)
    QQ = eye(freq, freq) * dt # shocks are Brownian motion -> variance is dt
    EE = zeros(1,1) * dt
    DD = [0.]
    dt = 1. / freq

    if track_lag
        ZZ = zeros(1, size(inverse_basis,2) * (freq + 1))

        # Numerically approximate integral of expm(T*dt) * s_0 using the formula
        # dt * (I + (TT * dt)/2! + (TT * dt)^2/3! + ...)
        TT_prev_power = TT * dt^2 # save previous power so can recursively multiply
        polynomial = eye(TT) * dt + TT_prev_power/2 # initialize polynomial
        for i = 2:max_iter+1
            TT_prev_power = TT_prev_power * TT * dt # compute next power

            # Check if error is small enough
            if maximum(abs.(TT_prev_power./factorial(i+1))) < abs_err
                break
            end

            # Otherwise add onto polynomial
            polynomial += TT_prev_power/factorial(i+1)
        end

        # Change basis from reduced basis to full basis on left side
        integ_output = inverse_basis * polynomial

        # grab row corresponding to output, fill in ZZ
        integ_output = integ_output[endo[:output], :]
        for i = 1:freq
            ZZ[1, 1+size(inverse_basis,2)*(i-1):size(inverse_basis,2)*i] = integ_output
        end

        # Change of basis into reduced basis
        # augment_inverse_basis = zeros(size(inverse_basis,1), size(inverse_basis,2) * (freq+1))
        # for i = 1:(freq+1)
        #     augment_inverse_basis[:, 1+size(inverse_basis,2)*(i-1):size(inverse_basis,2)*i] = inverse_basis
        # end
    else
        ZZ = zeros(1, size(inverse_basis,2) * freq)

        # Numerically approximate integral of expm(T*(-dt)) * s_t using the formula
        # -dt * (I + (TT * -dt)/2! + (TT * -dt)^2/3! + ...)
        # Note the negative dt is to conduct a backward motion
        TT_prev_power = TT * dt^2 # save previous power so can recursively multiply
        polynomial = eye(TT) * (-dt) + TT_prev_power/2
        for i = 2:max_iter+1
            TT_prev_power = TT_prev_power * TT * (-dt) # compute next power

            # Check if error is small enough
            if maximum(abs.(TT_prev_power./factorial(i+1))) < abs_err
                break
            end

            # Otherwise add onto polynomial
            polynomial += TT_prev_power/factorial(i+1)
        end

        # Change basis from reduced basis to full basis on left side
        integ_output = inverse_basis * polynomial

        # grab row corresponding to output, fill in ZZ
        integ_output = integ_output[endo[:output], :]
        for i = 1:freq
            ZZ[1, 1+size(inverse_basis,1)*(i-1):size(inverse_basis,1)*i] = integ_output
        end

        # Change of basis into reduced basis
        # augment_inverse_basis = zeros(size(inverse_basis,1), size(inverse_basis,2) * freq)
        # for i = 1:freq
        #     augment_inverse_basis[:, 1+size(inverse_basis,2)*(i-1):size(inverse_basis,2)*i] = inverse_basis
        # end
    end

    return Measurement(ZZ, DD, QQ, EE)
end


text:
#type definitions
abstract type Stem end

mutable struct cancercellM
    mutationsp::Array{Int64,1}
    mutationsd::Array{Int64,1}
    mutationsneg::Array{Int64,1}
    b::Float64
    d::Float64
    binitial::Float64
    dinitial::Float64
    fitness::Array{Float64, 1}
    fitnessneg::Array{Float64, 1}
    μp::Float64
    μd::Float64
    μneg::Float64
    timedrivers::Array{Float64, 1}
end

mutable struct SimResultM
    trueVAF::Array{Float64,1}
    trueVAFp::Array{Float64,1}
    trueVAFd::Array{Float64,1}
    trueVAFneg::Array{Float64,1}
    Ndrivers::Int64
    Npassengers::Int64
    Nnegative::Int64
    timedrivers::Array{Float64, 1}
    cells::Array{cancercellM, 1}
end

mutable struct InputParametersM
    Nmax::Int64
    detectionlimit::Float64
    ploidy::Int64
    read_depth::Float64
    clonalmutations::Int64
    μp::Float64
    μd::Float64
    μneg::Float64
    b::Float64
    d::Float64
    ρ::Float64
    cellularity::Float64
    s
    timefunction::Function
    fitnessfunc
end

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

function newmutations(cancercell::cancercellM, mutIDp, mutIDd, mutIDneg, Rmax, t, s, fitnessfunc)

    #function to add new mutations to cells based on μ
    numbermutationsp = rand(Poisson(cancercell.μp))
    numbermutationsd = rand(Poisson(cancercell.μd))
    numbermutationsneg = rand(Poisson(cancercell.μneg))
    cancercell.mutationsd = append!(cancercell.mutationsd, mutIDd:mutIDd+numbermutationsd-1)
    mutIDd = mutIDd + numbermutationsd

    cancercell.mutationsneg = append!(cancercell.mutationsneg, mutIDneg:mutIDneg+numbermutationsneg-1)
    mutIDneg = mutIDneg + numbermutationsneg

    b = cancercell.binitial
    #increase fitness due to driver mutations
    for i in 1:numbermutationsd
      stemp = s()
      push!(cancercell.fitness, stemp)
      push!(cancercell.timedrivers, t)
    end

    #decrease fitness due to driver mutations
    for i in 1:numbermutationsneg
      push!(cancercell.fitnessneg, -s())
    end

    b = fitnessfunc(cancercell, b)

    cancercell.b = b

    cancercell.mutationsp = append!(cancercell.mutationsp, mutIDp:mutIDp+numbermutationsp-1)
    mutIDp = mutIDp + numbermutationsp

    if cancercell.b + cancercell.d > Rmax
      Rmax = cancercell.b + cancercell.d
    end

    return cancercell, mutIDp, mutIDd, mutIDneg, Rmax
end

function multiplicativefitness(cancercell, b)
  for i in cancercell.fitness
    b = b * (1 + i)
  end
  for i in cancercell.fitnessneg
    b = b * (1 + i)
  end
  return b
end

function nonmultiplicativefitness(cancercell, b)
  for i in cancercell.fitness
    b = cancercell.binitial * (1 + i)
  end
  return b
end

function initializesim(mup, mud, muneg, b, d)

    #initialize time to zero
    t = 0.0
    tvec = Float64[]
    push!(tvec,t)

    #population starts with one cell
    N = 1
    Nvec = Int64[]
    push!(Nvec,N)

    #Initialize array of cell type that stores mutations for each cell and their fitness type
    #fitness type of 1 is the host population, lowest fitness
    cells = cancercellM[]
    push!(cells, cancercellM(Int64[], Int64[], Int64[], b, d, b, d, Float64[], Float64[], mup, mud, muneg, []))

    #need to keep track of mutations, assuming infinite sites, new mutations will be unique,
    #we assign each new muation with a unique integer by simply counting up from one
    mutIDp, mutIDd, mutIDneg = 1, 1, 1
    return t, tvec, N, Nvec, cells, mutIDp, mutIDd, mutIDneg
end

function initializesim(mup, mud, muneg, b, d, Ncells)

    #initialize time to zero
    t = 0.0
    tvec = Float64[]
    push!(tvec,t)

    #population starts with one cell
    N = 1
    Nvec = Int64[]
    push!(Nvec,N)

    #Initialize array of cell type that stores mutations for each cell and their fitness type
    #fitness type of 1 is the host population, lowest fitness
    cells = cancercellM[]
    push!(cells, cancercellM(Int64[], Int64[], Int64[], b, d, b, d, Float64[], Float64[], mup, mud, muneg, []))

    for i in 1:Ncells - 1
        push!(cells, cancercellM(Int64[], Int64[], Int64[], b, d, b, d, Float64[], Float64[], mup, mud, muneg, []))
    end

    #need to keep track of mutations, assuming infinite sites, new mutations will be unique,
    #we assign each new muation with a unique integer by simply counting up from one
    mutIDp, mutIDd, mutIDneg = 1, 1, 1
    return t, tvec, N, Nvec, cells, mutIDp, mutIDd, mutIDneg
end

function tumourmoran(b, d, Nmax, μp, μd, μneg, maxt; clonalmutations = μp, timefunction::Function = exptime, s = sfunc, fitnessfunc = nonmultiplicativefitness)

    t, tvec, N, Nvec, cells, mutIDp, mutIDd, mutIDneg = initializesim(μp, μd, μneg, b, d, Nmax)
    Rmax = b
    wts = b.*ones(Nmax)

    while t < maxt
      celldie = rand(1:Nmax)
      cellbirth = sample(1:Nmax, weights(wts))
      cells[celldie] = copycell(cells[cellbirth])
      cells[cellbirth], mutIDp, mutIDd, mutIDneg, Rmax = newmutations(cells[cellbirth], mutIDp, mutIDd, mutIDneg, Rmax, t, s, nonmultiplicativefitness)
      cells[celldie], mutIDp, mutIDd, mutIDneg, Rmax = newmutations(cells[celldie], mutIDp, mutIDd, mutIDneg, Rmax, t, s, nonmultiplicativefitness)
      Δt =  1/(Nmax) * timefunction()
      t = Δt + t
      push!(tvec, t)
      wts[celldie] = cells[celldie].b
      wts[cellbirth] = cells[cellbirth].b
    end

    return cells, tvec, Rmax
end

function copycell(cancercellold::cancercellM)
  newcancercell::cancercellM = cancercellM(copy(cancercellold.mutationsp),
  copy(cancercellold.mutationsd),
  copy(cancercellold.mutationsneg),
  copy(cancercellold.b),
  copy(cancercellold.d),
  copy(cancercellold.binitial),
  copy(cancercellold.dinitial),
  copy(cancercellold.fitness),
  copy(cancercellold.fitnessneg),
  copy(cancercellold.μp),
  copy(cancercellold.μd),
  copy(cancercellold.μneg),
  copy(cancercellold.timedrivers))
end

sfunc() = 0.0

function tumourgrow_birthdeath(b, d, Nmax, μp, μd, μneg; clonalmutations = μp, timefunction::Function = exptime, s = sfunc, fitnessfunc = nonmultiplicativefitness)

    #Rmax starts with b + d and changes once a fitter mutant is introduced, this ensures that
    # b and d have correct units
    Rmax = b + d

    #initialize arrays and parameters
    t, tvec, N, Nvec, cells, mutIDp, mutIDd, mutIDneg = initializesim(μp, μd, μneg, b, d)
    muts = Int64[]
    push!(muts, mutIDp)

    while N < Nmax

        #pick a random cell
        randcell = rand(1:N)
        r = rand(Uniform(0, Rmax))
	      Nt = N
        Rmaxt = Rmax

        #birth event if r<birthrate, access correct birthrate from cells array
        if r < cells[randcell].b

            #population increases by one
            N = N + 1
            #copy cell and mutations for cell that reproduces
            push!(cells,copycell(cells[randcell]))
            #add new mutations to both new cells
            cells[randcell], mutIDp, mutIDd, mutIDneg, Rmax = newmutations(cells[randcell], mutIDp, mutIDd, mutIDneg, Rmax, t, s,
            fitnessfunc)
            cells[end], mutIDp, mutIDd, mutIDneg, Rmax = newmutations(cells[end], mutIDp, mutIDd, mutIDneg, Rmax, t, s, fitnessfunc)
            push!(Nvec, N)
            Δt =  1/(Rmaxt * Nt) * timefunction()
            t = t + Δt
            push!(tvec,t)
        end

        #nothing if r > b+d
        if ((cells[randcell].b + cells[randcell].d) <= r )
          push!(Nvec, N)
          Δt =  1/(Rmax * Nt) * timefunction()
          t = t + Δt
          push!(tvec,t)
        end

        #death event if b<r<b+d
        if (cells[randcell].b <= r < (cells[randcell].b + cells[randcell].d))
            #population decreases by 1
            N = N - 1
            #frequency of cell type decreases
            #remove deleted cell
            deleteat!(cells,randcell)
            push!(Nvec,N)
            Δt =  1/(Rmax * Nt) * timefunction()
            t = t + Δt
            push!(tvec,t)
        end

        #every cell dies reinitialize simulation
        if (N == 0)
            t, tvec, N, Nvec, cells, mutIDp, mutIDd, mutIDneg = initializesim(μp, μd, μneg, b, d)
            muts = Int64[]
            push!(muts, mutIDp)
        end
    end

    return cells, tvec, Rmax
    #RawOutput(Nvec, tvec, muts, cells, birthrates, deathrates, clonetype, clonetime, subclonemutations, cloneN, Ndivisions, aveDivisions)
end

function cellsconvert(cells::Array{cancercellM, 1})
    #convert from array of cell types to one array with mutations and one array with cell fitness
    fitnessb = zeros(Float64,length(cells))
    mutationsp = Int64[]
    mutationsd = Int64[]
    mutationsneg = Int64[]
    drivertime = Float64[]
    fitness = Float64[]

    for i in 1:length(cells)
        append!(mutationsp,cells[i].mutationsp)
        append!(mutationsd,cells[i].mutationsd)
        append!(mutationsneg,cells[i].mutationsneg)
        fitnessb[i] = cells[i].b
        append!(drivertime, cells[i].timedrivers)
    end

    return mutationsp, mutationsd, mutationsneg, fitnessb, sort(unique(drivertime))
end

function getresults(b, d, μp, μd, μneg, Nmax; ploidy = 2, clonalmutations = 100, timefunction = exptime, s = sfunc, fitnessfunc = nonmultiplicativefitness)

    #Nvec,tvec,mvec,cells,br,dr,ct,clonetime
    cells, tvec, Rmax = tumourgrow_birthdeath(b, d, Nmax, μp, μd, μneg; clonalmutations = 0, timefunction = timefunction, s = s, fitnessfunc = fitnessfunc);
    Mp, Md, Mneg, fitness, timedrivers = cellsconvert(cells)

    return Mp, Md, Mneg, fitness, tvec, cells, timedrivers
end

function getresults(b, d, μp, μd, μneg, Nmax, tmax; ploidy = 2, clonalmutations = 100, timefunction = exptime, s = sfunc)
    #Nvec,tvec,mvec,cells,br,dr,ct,clonetime
    cells, tvec, Rmax = tumourmoran(b, d, Nmax, μp, μd, μneg, tmax; clonalmutations = 0, timefunction = timefunction, s = s);
    Mp, Md, Mneg, fitness, timedrivers = cellsconvert(cells)

    return Mp, Md, Mneg, fitness, tvec, cells, timedrivers
end

function run1simulation(IP::InputParametersM)

    Mp, Md, Mneg, fitness, tvec, cells, timedrivers = getresults(IP.b, IP.d, IP.μp, IP.μd, IP.μneg, IP.Nmax; ploidy = IP.ploidy, clonalmutations = IP.clonalmutations, timefunction = IP.timefunction, s = IP.s, fitnessfunc = IP.fitnessfunc)

    if length(Md) > 0
      AFd = counts(Md)
      filter!(x -> x>0, AFd)
    else
      AFd = []
    end

    if length(Mneg) > 0
      AFneg = counts(Mneg)
      filter!(x -> x>0, AFneg)
    else
      AFneg = []
    end

    if length(Mp) > 0
      AFp = counts(Mp)
      filter!(x -> x>0, AFp)
    else
      AFp = []
    end

    AFallmuts = prepend!([AFd; AFp; AFneg], repeat([Float64(IP.Nmax)], inner = IP.clonalmutations))

    return SimResultM(convert(Array{Float64, 1}, AFallmuts),
    convert(Array{Float64, 1}, AFp),
    convert(Array{Float64, 1}, AFd),
    convert(Array{Float64, 1}, AFneg),
    length(AFd),
    length(AFp),
    length(AFneg),
    timedrivers,
    cells), IP
end

function run1simulation(IP::InputParametersM, tmax)

    Mp, Md, Mneg, fitness, tvec, cells, timedrivers = getresults(IP.b, IP.d, IP.μp, IP.μd, IP.μneg, IP.Nmax, tmax; ploidy = IP.ploidy, clonalmutations = IP.clonalmutations, timefunction = IP.timefunction, s = IP.s)

    if length(Md) > 0
      AFd = counts(Md)
      filter!(x -> x>0, AFd)
    else
      AFd = []
    end

    if length(Mneg) > 0
      AFneg = counts(Mneg)
      filter!(x -> x>0, AFneg)
    else
      AFneg = []
    end

    if length(Mp) > 0
      AFp = counts(Mp)
      filter!(x -> x>0, AFp)
    else
      AFp = []
    end

    AFallmuts = prepend!([AFd; AFp; AFneg], repeat([Float64(IP.Nmax)], inner = IP.clonalmutations))

    return SimResultM(convert(Array{Float64, 1}, AFallmuts),
    convert(Array{Float64, 1}, AFp),
    convert(Array{Float64, 1}, AFd),
    convert(Array{Float64, 1}, AFneg),
    length(AFd),
    length(AFp),
    length(AFneg),
    timedrivers,
    cells), IP
end

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


text:
@adjoint function scatter_add!(ys::AbstractArray, us::AbstractArray, xs::AbstractArray)
    ys_ = copy(ys)
    scatter_add!(ys_, us, xs)
    ys_, Δ -> (Δ, gather(Δ, xs), nothing)
end

@adjoint function scatter_sub!(ys::AbstractArray, us::AbstractArray, xs::AbstractArray)
    ys_ = copy(ys)
    scatter_sub!(ys_, us, xs)
    ys_, Δ -> (Δ, -gather(Δ, xs), nothing)
end

@adjoint function scatter_mul!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where {T<:Real}
    ys_ = copy(ys)
    scatter_mul!(ys_, us, xs)
    ys_, function (Δ)
        Δy = Δ .+ zero(ys)
        scatter_mul!(Δy, us, xs)
        rev_xs = gather_indices(xs)
        Δu = gather(ys, xs) .* gather(Δ, xs)
        @inbounds for ind = CartesianIndices(xs)
            inds = filter(x -> x != ind, rev_xs[xs[ind]])
            for i = 1:size(us, 1)
                Δu[i, ind] *= prod(j -> us[i, j], inds)
            end
        end
        (Δy, Δu, nothing)
    end
end

@adjoint function scatter_div!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where {T<:Real}
    ys_ = copy(ys)
    scatter_div!(ys_, us, xs)
    ys_, function (Δ)
        Δy = Δ .+ zero(ys)
        scatter_div!(Δy, us, xs)
        rev_xs = gather_indices(xs)
        Δu = - gather(ys, xs) .* gather(Δ, xs) ./ us.^2
        @inbounds for ind = CartesianIndices(xs)
            inds = filter(x -> x != ind, rev_xs[xs[ind]])
            for i = 1:size(us, 1)
                Δu[i, ind] /= prod(j -> us[i, j], inds)
            end
        end
        (Δy, Δu, nothing)
    end
end

@adjoint function scatter_max!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where {T<:Real}
    max = copy(ys)
    scatter_max!(max, us, xs)
    max, function (Δ)
       Δy = (ys .== max) .* Δ
       Δu = (us .== gather(max, xs)) .* gather(Δ, xs)
       (Δy, Δu, nothing)
    end
end

@adjoint function scatter_min!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where {T<:Real}
    min = copy(ys)
    scatter_min!(min, us, xs)
    min, function (Δ)
       Δy = (ys .== min) .* Δ
       Δu = (us .== gather(min, xs)) .* gather(Δ, xs)
       (Δy, Δu, nothing)
    end
end

@adjoint function scatter_mean!(ys::AbstractArray, us::AbstractArray, xs::AbstractArray)
    ys_ = copy(ys)
    scatter_mean!(ys_, us, xs)
    ys_, function (Δ)
        Δu = gather(Δ, xs)
        counts = zero.(xs)
        @inbounds for i = 1:size(ys, 2)
            counts += sum(xs.==i) * (xs.==i)
        end
        @inbounds for ind = CartesianIndices(counts)
            view(Δu, :, ind) ./= counts[ind]
        end
        (Δ, Δu, nothing)
    end
end


text:
function gen_sys_ltc_mtk(wiring; T=Float32, name)
  @variables t
  D = Differential(t)

  n_neurons = wiring.n_total
  n_in = s_in = wiring.s_in
  w_sens = wiring.matrices[:w_sens]
  w_syns = wiring.matrices[:w_syns]

  sts = @variables begin
    (v[1:n_neurons](t) = rand_uniform(T, -0.3, 0.3, n_neurons)), [lower=T(-1), upper=T(1)]
    I_sens[1:n_neurons](t)
    I_syns[1:n_neurons](t)
    I_L[1:n_neurons](t)
  end

  ps = @parameters begin
    stim[1:n_in] = fill(T(1), n_in)
    C[1:n_neurons] = rand_uniform(T, 1, 1.01, n_neurons), [lower=T(0.9), upper=T(1.1)]
    G_L[1:n_neurons] = rand_uniform(T, 0.001, 0.2, n_neurons), [lower=T(0.001), upper=T(1)]
    E_L[1:n_neurons] = rand_uniform(T, -0.5, 0.5, n_neurons), [lower=T(-1), upper=T(1)]
  end

  sts = reduce(vcat, collect.(sts))
  ps = reduce(vcat, collect.(ps))

  eqs = Equation[]
  for i in 1:n_neurons
    push!(eqs, I_L[i] ~ G_L[i] * (v[i] - E_L[i]))
    push!(eqs, D(v[i]) ~ -(I_L[i] + I_sens[i] + I_syns[i]) / C[i])
  end

  systems = ODESystem[]

  I_senss = reshape(Num[0 for _ in 1:s_in*n_neurons], s_in, n_neurons)
  I_synss = reshape(Num[0 for _ in 1:n_neurons*n_neurons], n_neurons, n_neurons)

  for d in 1:n_neurons

    n_synapses = Int(sum(w_sens[:,d]))
    weights_exc_sens = Num[]
    weights_inh_sens = Num[]

    ss = Symbol("G_exc_sens$(d)")
    G_exc_sens = @parameters $ss = rand_uniform(T, 0.001, 0.2), [lower=T(0.001), upper=T(1)]
    push!(ps, G_exc_sens[1])
    ss = Symbol("G_inh_sens$(d)")
    G_inh_sens = @parameters $ss = rand_uniform(T, 0.001, 0.2), [lower=T(0.001), upper=T(1)]
    push!(ps, G_inh_sens[1])

    ss = Symbol("sum_weights_exc_sens$(d)")
    sum_weights_exc_sens = @variables $ss(t)
    push!(sts, sum_weights_exc_sens[1])
    ss = Symbol("sum_weights_inh_sens$(d)")
    sum_weights_inh_sens = @variables $ss(t)
    push!(sts, sum_weights_inh_sens[1])

    for s in 1:n_in
      w_sens[s,d] == 0 && continue

      E_def = wiring.matrices[:E_sens][s,d]

      ss = Symbol("weight_sensp$(s)_$(d)")
      weight_sensp = @parameters $ss = T(1/n_synapses), [lower=T(0.0001), upper=T(1)]
      ss = Symbol("weight_sens$(s)_$(d)")
      weight_sens = @variables $ss(t)

      push!(eqs, weight_sens[1] ~ weight_sensp[1])

      g_sens = 0

      if E_def < 0
        push!(weights_exc_sens, weight_sens[1])
        g_sens = (weight_sens[1] / sum_weights_exc_sens[1]) * G_exc_sens[1]
      else
        push!(weights_inh_sens, weight_sens[1])
        g_sens = (weight_sens[1] / sum_weights_inh_sens[1]) * G_inh_sens[1]
      end
      #g_sens = weight_sens[1] * G_sens[1]

      ss = Symbol("vh_sens$(s)_$(d)")
      vh_sens = @parameters $ss = rand_uniform(T, -0.4, -0.2) * E_def, [lower=T(-2), upper=T(2)]
      ss = Symbol("k_sens$(s)_$(d)")
      k_sens = @parameters $ss = rand_uniform(T, 2, 4), [lower=T(-20), upper=T(20)]

      ss = Symbol("E_sens$(s)_$(d)")
      #E_sens = @parameters $ss = rand_uniform(T, -0.5, 0.5), [lower=T(-1), upper=T(1)]
      E = @parameters $ss = E_def, [lower=T(-1), upper=T(1)]

      I_senss[s,d] = g_sens * mysigm(k_sens[1] * (fff(stim[s]) - vh_sens[1])) * (v[d] - E[1])

      push!(sts, weight_sens[1])
      push!(ps, vh_sens[1], k_sens[1], weight_sensp[1], E[1])

    end
    #push!(eqs, 0 ~ sum(ws->ws, weights_sens) - 1)
    push!(eqs, sum_weights_exc_sens[1] ~ sum(weights_exc_sens))
    push!(eqs, sum_weights_inh_sens[1] ~ sum(weights_inh_sens))

    n_synapses = Int(sum(w_syns[:,d]))
    weights_exc_syns = Num[]
    weights_inh_syns = Num[]

    ss = Symbol("G_exc_syns$(d)")
    G_exc_syns = @parameters $ss = rand_uniform(T, 0.001, 0.2), [lower=T(0.001), upper=T(1)]
    push!(ps, G_exc_syns[1])
    ss = Symbol("G_inh_syns$(d)")
    G_inh_syns = @parameters $ss = rand_uniform(T, 0.001, 0.2), [lower=T(0.001), upper=T(1)]
    push!(ps, G_inh_syns[1])

    ss = Symbol("sum_weights_exc_syns$(d)")
    sum_weights_exc_syns = @variables $ss(t)
    push!(sts, sum_weights_exc_syns[1])
    ss = Symbol("sum_weights_inh_syns$(d)")
    sum_weights_inh_syns = @variables $ss(t)
    push!(sts, sum_weights_inh_syns[1])

    for s in 1:n_neurons
      w_syns[s,d] == 0 && continue

      E_def = wiring.matrices[:E_syns][s,d]

      ss = Symbol("weight_synsp$(s)_$(d)")
      weight_synsp = @parameters $ss = T(1/n_synapses), [lower=T(0.0001), upper=T(1)]
      ss = Symbol("weight_syns$(s)_$(d)")
      weight_syns = @variables $ss(t)

      push!(eqs, weight_syns[1] ~ weight_synsp[1])

      g_syns = 0

      if E_def < 0
        push!(weights_exc_syns, weight_syns[1])
        g_syns = (weight_syns[1] / sum_weights_exc_syns[1]) * G_exc_syns[1]
      else
        push!(weights_inh_syns, weight_syns[1])
        g_syns = (weight_syns[1] / sum_weights_inh_syns[1]) * G_inh_syns[1]
      end
      #g_syns = weight_syns[1] * G_syns[1]

      ss = Symbol("vh_syns$(s)_$(d)")
      vh_syns = @parameters $ss = rand_uniform(T, -0.4, -0.2) * E_def, [lower=T(-2), upper=T(2)]
      ss = Symbol("k_syns$(s)_$(d)")
      k_syns = @parameters $ss = rand_uniform(T, 2, 4), [lower=T(-20), upper=T(20)]

      ss = Symbol("E_syns$(s)_$(d)")
      #E_syns = @parameters $ss = rand_uniform(T, -0.5, 0.5), [lower=T(-1), upper=T(1)]
      E = @parameters $ss = E_def, [lower=T(-1), upper=T(1)]

      I_synss[s,d] = g_syns * mysigm(k_syns[1] * (v[s] - vh_syns[1])) * (v[d] - E[1])

      push!(sts, weight_syns[1])
      push!(ps, vh_syns[1], k_syns[1], weight_synsp[1], E[1])
    end
    #push!(eqs, 0 ~ sum(ws->ws, weights_syns) - 1)
    push!(eqs, sum_weights_exc_syns[1] ~ sum(weights_exc_syns))
    push!(eqs, sum_weights_inh_syns[1] ~ sum(weights_inh_syns))
  end

  for i in 1:n_neurons
    push!(eqs, I_sens[i] ~ sum(I_senss[:,i]))
  end
  for i in 1:n_neurons
    push!(eqs, I_syns[i] ~ sum(I_synss[:,i]))
  end

  ODESystem(eqs, t, sts, ps; systems, name)
end

#Zygote.@nograd gen_sys_ltc_syn_state

function LTCMTK(wiring, solver, sensealg;
  T=Float32, tspan=T.((0, 1)),
  #wiring=nothing,
  tgrad=true, gen_jac=false, kwargs...)

  @named sys = gen_sys_ltc_mtk(wiring)
  ssys = structural_simplify(sys)
  defs = ModelingToolkit.get_defaults(ssys)

  prob = ODEProblem(ssys, defs, tspan, tgrad=tgrad, jac=gen_jac)

  lb, ub = BCTRNN.get_bounds(ssys, wiring.s_in; T)

  rnncell = BCTRNNCell(wiring, solver, sensealg, prob, lb, ub; kwargs...)
  MyRecur(rnncell), defs, ssys
end


text:
@enum LDAPVersion::Cint begin
    LDAP_VERSION1 = 1
    LDAP_VERSION2 = 2
    LDAP_VERSION3 = 3
end

@enum LDAPOption::UInt16 begin
    LDAP_OPT_PROTOCOL_VERSION = 0x0011
end

@enum LDAPMessageType::UInt8 begin
    LDAP_RES_BIND = 0x61
    LDAP_RES_SEARCH_ENTRY = 0x64
    LDAP_RES_SEARCH_REFERENCE = 0x73
    LDAP_RES_SEARCH_RESULT = 0x65
    LDAP_RES_MODIFY = 0x67
    LDAP_RES_ADD = 0x69
    LDAP_RES_DELETE = 0x6b
    LDAP_RES_MODDN = 0x6d
    LDAP_RES_COMPARE = 0x6f
    LDAP_RES_EXTENDED = 0x78
    LDAP_RES_INTERMEDIATE = 0x79
end

@enum LDAPScope::Cint begin
    LDAP_SCOPE_BASE = 0x0000
    LDAP_SCOPE_ONELEVEL = 0x0001
    LDAP_SCOPE_SUBTREE = 0x0002
    LDAP_SCOPE_CHILDREN = 0x0003
end


text:
function transform!(context::T) where {T<:Union{SHA2_224_CTX,SHA2_256_CTX}}
    pbuf = buffer_pointer(context)
    # Initialize registers with the previous intermediate values (our state)
    a = context.state[1]
    b = context.state[2]
    c = context.state[3]
    d = context.state[4]
    e = context.state[5]
    f = context.state[6]
    g = context.state[7]
    h = context.state[8]

    # Run initial rounds
    for j = 1:16
        @inbounds begin
            # We bitswap every input byte
            v = bswap(unsafe_load(pbuf, j))
            unsafe_store!(pbuf, v, j)

            # Apply the SHA-256 compression function to update a..h
            T1 = h + Sigma1_256(e) + Ch(e, f, g) + K256[j] + v
            T2 = Sigma0_256(a) + Maj(a, b, c)
            h = g
            g = f
            f = e
            e = UInt32(d + T1)
            d = c
            c = b
            b = a
            a = UInt32(T1 + T2)
        end
    end

    for j = 17:64
        @inbounds begin
            # Implicit message block expansion:
            s0 = unsafe_load(pbuf, mod1(j + 1, 16))
            s0 = sigma0_256(s0)
            s1 = unsafe_load(pbuf, mod1(j + 14, 16))
            s1 = sigma1_256(s1)

            # Apply the SHA-256 compression function to update a..h
            v = unsafe_load(pbuf, mod1(j, 16)) + s1 + unsafe_load(pbuf, mod1(j + 9, 16)) + s0
            unsafe_store!(pbuf, v, mod1(j, 16))
            T1 = h + Sigma1_256(e) + Ch(e, f, g) + K256[j] + v
            T2 = Sigma0_256(a) + Maj(a, b, c)
            h = g
            g = f
            f = e
            e = UInt32(d + T1)
            d = c
            c = b
            b = a
            a = UInt32(T1 + T2)
        end
    end

    # Compute the current intermediate hash value
    context.state[1] += a
    context.state[2] += b
    context.state[3] += c
    context.state[4] += d
    context.state[5] += e
    context.state[6] += f
    context.state[7] += g
    context.state[8] += h
end

function transform!(context::Union{SHA2_384_CTX,SHA2_512_CTX})
    pbuf = buffer_pointer(context)
    # Initialize registers with the prev. intermediate value
    a = context.state[1]
    b = context.state[2]
    c = context.state[3]
    d = context.state[4]
    e = context.state[5]
    f = context.state[6]
    g = context.state[7]
    h = context.state[8]

    for j = 1:16
        @inbounds begin
            v = bswap(unsafe_load(pbuf, j))
            unsafe_store!(pbuf, v, j)

            # Apply the SHA-512 compression function to update a..h
            T1 = h + Sigma1_512(e) + Ch(e, f, g) + K512[j] + v
            T2 = Sigma0_512(a) + Maj(a, b, c)
            h = g
            g = f
            f = e
            e = d + T1
            d = c
            c = b
            b = a
            a = T1 + T2
        end
    end

    for j = 17:80
        @inbounds begin
            # Implicit message block expansion:
            s0 = unsafe_load(pbuf, mod1(j + 1, 16))
            s0 = sigma0_512(s0)
            s1 = unsafe_load(pbuf, mod1(j + 14, 16))
            s1 = sigma1_512(s1)

            # Apply the SHA-512 compression function to update a..h
            v = unsafe_load(pbuf, mod1(j, 16)) + s1 + unsafe_load(pbuf, mod1(j + 9, 16)) + s0
            unsafe_store!(pbuf, v, mod1(j, 16))
            T1 = h + Sigma1_512(e) + Ch(e, f, g) + K512[j] + v
            T2 = Sigma0_512(a) + Maj(a, b, c)
            h = g
            g = f
            f = e
            e = d + T1
            d = c
            c = b
            b = a
            a = T1 + T2
        end
    end

    # Compute the current intermediate hash value
    context.state[1] += a
    context.state[2] += b
    context.state[3] += c
    context.state[4] += d
    context.state[5] += e
    context.state[6] += f
    context.state[7] += g
    context.state[8] += h
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
# # Prestressed simply-supported column modal analysis

# Source code: [`prestressed_column_modal_tut.jl`](prestressed_column_modal_tut.jl)

# ## Description

# Vibration analysis of a simply supported column loaded with axial force.
# The fundamental vibration frequency depends on the prestress force.

# ## Goals

# - Demonstrate the change in the fundamental vibration frequency due to the
#   prestress through an axial force.
# - Show how to evaluate the geometric stiffness from a static solution.
# - Solve the eigenvalue free vibration problem with the inclusion of the
#   geometric stiffness matrix.

# ## Definition of the basic inputs

# The finite element code realize on the basic functionality implemented in this
# package.
using FinEtools

# The material parameters may be defined with the specification of the units.
# The elastic properties are:
E = 30002.0 * phun("ksi")
nu = 0.0;

# The mass density is expressed in customary units as
g = 32.17*12 * phun("in/sec^2")
rho = 4.65 * phun("oz/in ^3") / g
# Here are the cross-sectional dimensions and the length of the beam between supports.
b = 1.8 * phun("in"); h = 1.8 * phun("in"); L = 300 * phun("in");

# ## Cross-section

# Cross-sectional properties are incorporated in the cross-section property. The
# three arguments supplied are functions. All are returning "constants". In
# particular the first two functions each return the dimension of the
# cross-section as a constant(the beam has a uniform cross-section); the third
# function defines the orientation of the cross-section in the global Cartesian
# coordinates. `[1.0, 0.0, 0.0]` is the vector that together with the tangent
# to the midline curve of the beam spans the $x_1x_2$ plane of the local
# coordinates for the beam.
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
cs = CrossSectionRectangle(s -> b, s -> h, s -> [1.0, 0.0, 0.0])

# Here we retrieve the cross-sectional properties at the arc length 0.0.
@show A, J, I1, I2, I3 = cs.parameters(0.0)

# ## Analytical frequencies

# The analytical frequencies were taken from table 8-1 of Formulas for
# natural frequency and mode shape, Robert D. Blevins, Krieger publishing
# company, Malabar Florida, reprint edition 2001.
#
# The beam has cylindrical supports at either end.
#
# Simply supported column without a pre-stressing force has a fundamental frequency of
@show analyt_freq = (1*pi)^2/(2*pi*L^2)*sqrt(E*I2/rho/A);

# The critical Euler buckling load (simple-support): (pi^2*E*I2/L^2).
@show PEul = (pi^2*E*I2/L^2);

# Ps=linspace(-1, 1, 100);

# The purpose of the numerical model is to calculate approximation to the fundamental
# analytical natural frequency.
neigvs = 1;

# Now we generate the mesh of the beam. The locations of its two endpoints are:
xyz = [[0 -L/2 0]; [0 L/2 0]]
# We will generate
n = 20
# beam elements along the member.
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member
fens, fes = frame_member(xyz, n, cs);

# ## Material

# Material properties can be now used to create a material: isotropic elasticity model of the `FinEtoolsDeforLinear` package is instantiated.
using FinEtoolsDeforLinear
material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

# ## Fields

# Now we start constructing the discrete finite element model.
# We begin by constructing the requisite fields, geometry and displacement.
# These are the so-called "configuration variables", all initialized to 0.
# This is that geometry field.
geom0 = NodalField(fens.xyz)
# This is the displacement field, three unknown displacements per node.
u0 = NodalField(zeros(size(fens.xyz, 1), 3))
# This is the rotation field, three unknown rotations per node are represented
# with a rotation matrix, in total nine numbers. The utility function
# `initial_Rfield`
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield
Rfield0 = initial_Rfield(fens)

# Finally, this is the displacement and rotation field for incremental changes,
# incremental displacements and incremental rotations. In total, 6 unknowns per
# node.
dchi = NodalField(zeros(size(fens.xyz, 1), 6))

# ## Support conditions

# Now we apply the essential boundary conditions (EBCs) to enforce the action of
# the supports at the ends of the beam.

# First we select the node at the location  `[0 -L/2 0]`. This is the immovable node.
immovable = selectnode(fens; box=[0 0 -L/2 -L/2 0 0], tolerance=L/n/1000)[1]

# The boundary condition at this point dictates zero displacements (degrees of
# freedom 1, 2, and 3) and zero rotations about the axis of the beam (degree of freedom 5).
for i in [1,2,3,5]
    setebc!(dchi, [immovable], true, i)
end

# Similarly, the node next to the other end of the beam is selected. This node
# is to be exposed to the external loading through the axial force.
movable = selectnode(fens; box=[0 0 L/2  L/2 0 0], tolerance=L/n/1000)[1]

# This time the transverse displacements and the axial rotation are suppressed.
for i in [1,3,5]
    setebc!(dchi, [movable], true, i)
end

# These boundary conditions now need to be "applied". This simply means that the
# prescribed values of the degrees of freedom are copied into the active degrees
# of freedom.
applyebc!(dchi)
# The essential boundary conditions will also reduce the number of free
# (unknown) degrees of freedom.
numberdofs!(dchi);
# Here we inspect the degrees of freedom in the incremental
# displacement/rotation field:
@show dchi.dofnums
# Note that the degrees of freedom are actually carried by the incremental
# field, not by the displacement or the rotation fields.

# ## Assemble the global discrete system

using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material);

# For disambiguation we will refer to the stiffness and mass functions by
# qualifying them with the corotational-beam module, `FEMMCorotBeamModule`.
using FinEtoolsFlexStructures.FEMMCorotBeamModule
CB = FEMMCorotBeamModule

# Thus we can construct the stiffness matrix as follows:
# Note that the finite element machine is the first argument. This provides
# access to the integration domain. The next argument is the geometry field,
# followed by the displacement, rotations, and incremental
# displacement/rotation fields.
K = CB.stiffness(femm, geom0, u0, Rfield0, dchi);

# Now we construct the means of applying the concentrated force of prestress at
# the movable node. The node is the "boundary" of the domain of the column.
loadbdry = FESetP1(reshape([movable], 1, 1))
# The concentrated force can be considered a distributed loading at a single
# point. This distributed loading can be integrated with a quadrature rule
# suitable for a single point domains.
lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))

# We assume that the magnitude of the force is unity. Note that the force is
# applied in the positive direction along the beam, meaning as tensile.
# Therefore, buckling will be caused by applying  a negative buckling factor.
# Numerically, the buckling factor corresponding to failure is then of the
# magnitude of the Euler force.
fi = ForceIntensity(FFlt[0, 1.0, 0, 0, 0, 0]);

# The distributed loading is now integrated over the "volume" of the integration
# domain.
F = CB.distribloads(lfemm, geom0, dchi, fi, 3);

# Solve for the displacement under @show the static load
scattersysvec!(dchi, K\F);

# Update deflections so that the initial stress can be computed. First the
# displacements:
u1 = deepcopy(u0)
u1.values .= dchi.values[:, 1:3]
# Then the rotations:
Rfield1 = deepcopy(Rfield0)
using FinEtoolsFlexStructures.RotUtilModule:  update_rotation_field!
update_rotation_field!(Rfield1, dchi)

# The static deflection is now used to compute the internal forces
# which in turn lead to the geometric stiffness matrix.
Kg = CB.geostiffness(femm, geom0, u1, Rfield1, dchi);

# Now we can evaluate the mass matrix,
M = CB.mass(femm, geom0, u0, Rfield0, dchi);
# and we have the complete discrete model for the solution of the free vibration problem.

# ## Solve the free-vibration problem

# The Arnoldi algorithm implemented in the well-known `Arpack` package is used
# to solve the generalized eigenvalue problem with the sparse matrices. As is
# common in structural dynamics, we request the smallest eigenvalues in
# absolute value (`:SM`).
using Arpack

# The prestress-modified natural frequencies are computed in the loop for
# a number of prestress force values.
Ps = collect(linearspace(-1.0, 1.0, 50)).*PEul
freqs = let freqs =[];
    for P in Ps
        # Note that we take the complete stiffness matrix: elastic plus prestress (initial stress).
        evals, evecs, nconv = eigs(K + P.*Kg, M; nev=neigvs, which=:SM, explicittransform = :none);
        # The fundamental frequency:
        f = sqrt(evals[1]) / (2 * pi);
        push!(freqs, f);
    end
    freqs
end

# Show the normalized force and the fundamental frequencies.
sigdig(n) = round(n * 1000) / 1000
@show sigdig.(Ps./PEul), sigdig.(freqs)

# ## Present a plot

using Gnuplot

@gp  "set terminal windows 0 "  :-

@gp  :- Ps./PEul freqs./analyt_freq " lw 2 lc rgb 'red' with p title 'Fundamental frequency' "  :-

@gp  :- "set xlabel 'P/P_{Euler}'" :-
@gp  :- "set ylabel 'Frequency(P)/Frequency(0) [Hz]'" :-
@gp  :- "set title 'Prestressed column'"

nothing


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using Test
using ACE, JuLIP, JuLIP.Testing, QuadGK, LinearAlgebra, ACE.JacobiPolys
using ACE: TransformedJacobi, transform, transform_d,
             alloc_B, alloc_temp

using JuLIP: evaluate!

##

@noinline get_IN(N) = collect((aceB.idx_Bll[N][1]+1):(aceB.idx_Bll[N][end]+aceB.len_Bll[N][end]))

# function barrier
@noinline gramian(N, aceB, Nsamples=100_000; normalise=false) =
   gramian(N, get_IN(N), aceB, Nsamples; normalise=normalise)

@noinline gramian(N, IN, aceB, Nsamples; normalise=false) =
   gramian(N, IN, alloc_temp(aceB), alloc_B(aceB), aceB, Nsamples, normalise)

function gramian(N, IN, tmp, B, aceB, Nsamples = 100_000, normalise = false)
   Zs = zeros(Int16, N)
   lenB = length(IN)
   G = zeros(Float64, lenB, lenB)
   for n = 1:Nsamples
      Rs = ACE.rand_vec(aceB.J, N)
      evaluate!(B, tmp, aceB, Rs, Zs, 0)
      for j = 1:lenB
         Bj = B[IN[j]]'
         @simd for i = 1:lenB
            @inbounds G[i,j] += Bj * B[IN[i]]
         end
      end
   end
   if normalise
      g = diag(G)
      for i = 1:lenB, j = 1:lenB
         G[i,j] /= sqrt(g[i]*g[j])
      end
   end
   return G / Nsamples
end

##
Nmax = 4
Nsamples = 1_000
rl, ru = 0.5, 3.0
fcut =  PolyCutoff2s(2, rl, ru)
trans = PolyTransform(2, 1.0)
aceB = SHIPBasis( SparseSHIP(6, 15), trans, fcut )

ctr = 0
nnz = 0
for N = 1:Nmax
   for νz in aceB.NuZ[N]
      kk, ll = ACE._kl(νz.ν, νz.izz, aceB.KL)
      for mm in ACE._mrange(ll)
         global ctr += 1
         global nnz += size( aceB.rotcoefs[N][ll], 2 )
      end
   end
end
length(aceB)
ctr
nnz

function fill_an_array!(A)
   for n = 1:length(A)
      A[n] = rand()
   end
   return A
end

A = zeros(100_000);
(@elapsed fill_an_array!(A)) * 1000

@show aceB.idx_Bll[end]
@show aceB.len_Bll[end]
@show length(aceB)

G = gramian(4, aceB, 10_000)
rank(G)
svdG = svd(G)
svdG.U

@info("Conditions numbers of gramians")
for N = 1:Nmax
   GN = gramian(N, aceB, Nsamples, normalise=false)
   @show N, cond(GN)
end

@info("Conditions numbers of normalised gramians")
for N = 1:Nmax
   GN = gramian(N, aceB, Nsamples, normalise=true)
   @show N, cond(GN)
end


text:
using MPI, MPIArrays

MPI.Init()
rank = MPI.Comm_rank(MPI.COMM_WORLD)
N = 8 # size of the matrix

# Create an uninitialized matrix and vector
x = MPIArray{Float64}(N)
A = MPIArray{Float64}(N,N)

# Set random values
forlocalpart!(m -> fill!(m,1+rank),x)
forlocalpart!(m -> fill!(m,1+rank),A)

# Make sure every process finished initializing the coefficients
sync(A, x)

b = A*x

# This will print on the first process, using slow element-by-element communication, but that's OK to print to screen
rank == 0 && println("Matvec result: $b")

# Clean up
free(A)
free(x)
MPI.Finalize()


text:
abstract type AbstractFileInfo end

# type declarations

mutable struct FileCabinet{ V <: Vector{<:AbstractFileInfo} }
  file_infos::V
  boot_module::Module
  export_list::Vector{Symbol}
  load_order::Vector{AbstractString}
end

mutable struct FileInfo <: AbstractFileInfo
  name::AbstractString
  parent::FileCabinet{Vector{FileInfo}}
  undef::Union{Any, Symbol}
  time::AbstractFloat
  loaded_shards::AbstractArray{Expr}
  unloaded_shards::AbstractArray{Expr}
end

# constructor functions

function FileCabinet(cur_package::Module)
  cur_cabinet = FileCabinet(
    FileInfo[],
    cur_package.eval(parse("module BootDummy_$( camelize(string(Base.Random.uuid4())) ) end")),
    Symbol[],
    AbstractString[]
  )

  cur_cabinet
end

FileInfo{U <: AbstractFileInfo}(name::AbstractString, parent::FileCabinet{Vector{U}}) =
  FileInfo(
    name,
    parent,
    nothing,
    0.0,
    Expr[],
    Expr[]
  )


text:
# utility/z-test.jl

using Test

@test downsample(:test)
@test map_many(:test)

@test mask_or(:test)
@test mask_outline(:test)
@test embed(:test)
@test maskit(:test)
@test ndgrid(:test)


text:
function Laplacian{T<:AbstractFloat}(nz::Int,nx::Int,dz::T,dx::T)

# This function computes the 2D Laplacian operator/matrix that acts on a
# vectorized 2D field. The vectorization is assumed to be done using column
# major order.
#
# INPUTS:     nz        - Number of grid points in the z-direction
#             nx        - Number of grid points in the x-direction
#             dz        - Grid spacing in z-direction
#             dx        - Grid spacing in x-direction
#
# OUTPUTS:    L         - 2D Laplace operator

    Dzz = (1/dz^2)*spdiagm((ones(nz-1),-2*ones(nz),ones(nz-1)),(-1,0,1))
    Dxx = (1/dx^2)*spdiagm((ones(nx-1),-2*ones(nx),ones(nx-1)),(-1,0,1))

    Mzz = kron(speye(nx),Dzz)
    Mxx = kron(Dxx,speye(nz))

    L = Mzz + Mxx

    return L

end


text:
"""
    install_deps([env]; [force=false])

Uses Conda.jl to install HUMAnN and MetaPhlAn.
In order to use the commandline tools,
you must have the conda environment bin directory in `ENV["PATH"]`.
See "[Using Conda](@ref using-conda)" for more information.
"""
function install_deps(env=:BiobakeryUtils; force=false)
    if isdir(Conda.bin_dir(env)) && !force
        @warn """
        You already seem to have an environment, '$env'.
        If you've already installed the bioBakery packages, try

        `ENV["PATH"] = ENV["PATH"] * ":" * "$(Conda.bin_dir(env))"`

        Use `force=true` to install anyway"""
        return nothing
    end

    Conda.add_channel("bioconda", env)
    Conda.add_channel("conda-forge", env)
    Conda.add("humann", env)
    Conda.add("tbb=2020.2", env) # https://www.biostars.org/p/494922/

    @warn """
    Don't forget to add $(Conda.bin_dir(env)) to your PATH!

    This can be done in a julia session with:

    `ENV["PATH"] = ENV["PATH"] * ":" * "$(Conda.bin_dir(env))"`,
    or you can set it in your shell environment.
    """
    return nothing
end

function add_cli_kwargs!(cmd, kwargs; optunderscores=true)
    for (key,val) in pairs(kwargs)
        if val isa Bool
            val && push!(cmd, string("--", key))
        elseif val isa AbstractVector
            append!(cmd, [string("--", key), string.(val)...])
        else
            append!(cmd, [string("--", key), string(val)])
        end
    end
    !optunderscores && map!(c-> startswith(c, "--") ? replace(c, "_"=>"-") : c, cmd, cmd)
    return cmd
end

function check_for_install(tool)
    try run(pipeline(`which $tool`, stdout=devnull))
        return nothing
    catch e
        @error """
        Can not find `$tool`! If you think it should be installed,
        try running:

        ```
        ENV["PATH"] = ENV["PATH"] * ":" * Conda.bin_dir(env)
        ```

        Where `env` is something like `:BiobakeryUtils`.
        """
        rethrow()
    end
    return nothing
end


text:
"""
    LK(Args...)

A differential method for optical flow estimation developed by Bruce D. Lucas
and Takeo Kanade. It assumes that the flow is essentially constant in a local
neighbourhood of the pixel under consideration, and solves the basic optical flow
equations for all the pixels in that neighbourhood, by the least squares criterion.

The different arguments are:

 -  prev_points       =  Vector of SVector{2} for which the flow needs to be found
 -  next_points       =  Vector of SVector{2} containing initial estimates of new positions of
                         input features in next image
 -  window_size       =  Size of the search window at each pyramid level; the total size of the
                         window used is 2*window_size + 1
 -  max_level         =  0-based maximal pyramid level number; if set to 0, pyramids are not used
                         (single level), if set to 1, two levels are used, and so on
 -  estimate_flag     =  true -> Use next_points as initial estimate
                         false -> Copy prev_points to next_points and use as estimate
 -  term_condition    =  The termination criteria of the iterative search algorithm i.e the number of iterations
 -  min_eigen_thresh  =  The algorithm calculates the minimum eigenvalue of a (2 x 2) normal matrix of optical
                         flow equations, divided by number of pixels in a window; if this value is less than
                         min_eigen_thresh, then a corresponding feature is filtered out and its flow is not processed
                         (Default value is 10^-6)

## References

B. D. Lucas, & Kanade. "An Interative Image Registration Technique with an Application to Stereo Vision,"
DARPA Image Understanding Workshop, pp 121-130, 1981.

J.-Y. Bouguet, “Pyramidal implementation of the afﬁne lucas kanadefeature tracker description of the
algorithm,” Intel Corporation, vol. 5,no. 1-10, p. 4, 2001.
"""

struct LK{T <: Int, F <: Float64, V <: Bool} <: OpticalFlowAlgo
    prev_points::Array{SVector{2, T}, 1}
    next_points::Array{SVector{2, F}, 1}
    window_size::T
    max_level::T
    estimate_flag::V
    term_condition::T
    min_eigen_thresh::F
end

LK(prev_points::Array{SVector{2, T}, 1}, next_points::Array{SVector{2, F}, 1}, window_size::T, max_level::T, estimate_flag::V, term_condition::T) where {T <: Int, F <: Float64, V <: Bool} = LK{T, F, V}(prev_points, next_points, window_size, max_level, estimate_flag, term_condition, 0.000001)

function optflow(first_img::AbstractArray{T, 2}, second_img::AbstractArray{T,2}, algo::LK{}) where T <: Gray
    if algo.estimate_flag
        @assert size(algo.prev_points) == size(algo.next_points)
    end

    # Construct gaussian pyramid for both the images
    first_pyr = gaussian_pyramid(first_img, algo.max_level, 2, 1.0)
    second_pyr = gaussian_pyramid(second_img, algo.max_level, 2, 1.0)

    # Initilisation of the output arrays
    guess_flow = Array{SVector{2, Float64}, 1}(length(algo.prev_points))
    final_flow = Array{SVector{2, Float64}, 1}(length(algo.prev_points))
    status_array = trues(size(algo.prev_points))
    error_array = zeros(Float64, size(algo.prev_points))

    for i = 1:length(algo.prev_points)
        if algo.estimate_flag
            guess_flow[i] = SVector{2}((algo.next_points[i][1] - algo.prev_points[i][1])/2^algo.max_level,(algo.next_points[i][2] - algo.prev_points[i][2])/2^algo.max_level)
        else
            guess_flow[i] = SVector{2}(0.0,0.0)
        end
        final_flow[i] = SVector{2}(0.0,0.0)
    end
    min_eigen = 0.0
    grid_size = 0

    # Iterating over each level of the pyramids
    for i = algo.max_level+1:-1:1
        # Find x and y image gradients with scharr kernel (As written in Refernce [2])
        Ix, Iy = imgradients(first_pyr[i], KernelFactors.scharr, Fill(zero(eltype(first_pyr[i]))))

        # Extrapolate the second image so as to get intensity at non integral points
        itp = interpolate(second_pyr[i], BSpline(Quadratic(Flat())), OnGrid())
        etp = extrapolate(itp, zero(eltype(second_pyr[i])))

        # Calculate Ixx, Iyy and Ixy taking Gaussian weights for pixels in the grid
        Ixx = imfilter(Ix.*Ix, Kernel.gaussian(4))
        Iyy = imfilter(Iy.*Iy, Kernel.gaussian(4))
        Ixy = imfilter(Ix.*Iy, Kernel.gaussian(4))

        # Running the algorithm for each input point
        for j = 1:length(algo.prev_points)

            # Check if point is lost
            if status_array[j]
                # # Position of point in current pyramid level
                point = SVector{2}(floor(Int,(algo.prev_points[j][1])/2^i),floor(Int,(algo.prev_points[j][2])/2^i))
                # Bounds for the search window
                if all(x->isa(x, Base.OneTo), indices(first_pyr[i]))
                    rs = max(1,point[1] - algo.window_size):min(indices(first_pyr[i])[1].stop, point[1] + algo.window_size)
                    cs = max(1,point[2] - algo.window_size):min(indices(first_pyr[i])[2].stop,point[2] + algo.window_size)
                    grid = SVector{2}( rs, cs )
                else
                    rs = max(indices(first_pyr[i])[1].start,point[1] - algo.window_size):min(indices(first_pyr[i])[1].stop, point[1] + algo.window_size)
                    cs = max(indices(first_pyr[i])[2].start,point[2] - algo.window_size):min(indices(first_pyr[i])[2].stop, point[2] + algo.window_size)
                    grid = SVector{2}( rs, cs )
                end

                # Spatial Gradient Matrix
                G = SMatrix{2,2,Float64}(sum(Ixx[grid[1],grid[2]]),sum(Ixy[grid[1],grid[2]]), sum(Ixy[grid[1],grid[2]]), sum(Iyy[grid[1],grid[2]]) )

                temp_flow = SVector{2}(0.0,0.0)

                # Iterating till terminating condition
                for k = 1:algo.term_condition
                    # Diff flow between images to calculate the image difference
                    diff_flow = guess_flow[j] + temp_flow

                    # Check if the search window is present completely in the image
                    grid_1, grid_2, flag = get_grid(first_pyr[i], grid, SVector{2}(point), diff_flow, algo.window_size)
                    # Calculate the flow
                    A = @view first_pyr[i][grid_1[1],grid_1[2]]
                    B = @view etp[grid_2[1],grid_2[2]]
                    I_x = @view Ix[grid_1[1],grid_1[2]]
                    I_y = @view Iy[grid_1[1],grid_1[2]]
                    bx = 0.0
                    by = 0.0
                    P, Q = size(A)
                    # Evaluates sum(δI.*I_x) and sum(δI.*I_y) where δI = A .- B
                    @inbounds begin
                        for q = 1:Q
                           for p = 1:P
                                bx = bx + (A[p,q] - B[p,q])*I_x[p,q]
                                by = by + (A[p,q] - B[p,q])*I_y[p,q]
                            end
                        end
                    end
                    bk = SVector{2}(bx, by)

                    if !flag
                        G = SMatrix{2,2,Float64}(sum(Ixx[grid_1[1],grid_1[2]]),
                                         sum(Ixy[grid_1[1],grid_1[2]]),
                                         sum(Ixy[grid_1[1],grid_1[2]]),
                                         sum(Iyy[grid_1[1],grid_1[2]]))
                    end

                    # Use a custom Moore-Pensore inverse for a 2-by-2 matrix
                    # because the default  pinv(Array(G)) allocates memory.
                    # TODO: Add the pinv2x2 method to the StaticArrays package.
                    G_inv = pinv2x2(G)
                    ηk = G_inv*bk
                    temp_flow = ηk + temp_flow

                    # Minimum eigenvalue in the matrix is used as the error function
                    D, V =  eig(G)
                    min_eigen = D[1] < D[2] ? D[1] : D[2]
                    grid_size = (grid_1[1].stop - grid_1[1].start + 1) * (grid_1[2].stop - grid_1[2].start + 1)
                    error_array[j] = min_eigen/grid_size

                    # Check whether point is lost
                    if is_lost(first_pyr[i], SVector{2}(point + temp_flow), min_eigen, algo.min_eigen_thresh, grid_size)
                        status_array[j] = false
                        final_flow[j] = SVector{2}(0.0,0.0)
                        break
                    end
                end

                if status_array[j]
                    final_flow[j] = temp_flow
                end
                # Guess for the next level of pyramid
                guess_flow[j] = 2*(guess_flow[j] + final_flow[j])

                # Declare point lost if very high displacement
                if is_lost(SVector(0.5*guess_flow[j]), 2*algo.window_size)
                    status_array[j] = false
                    final_flow[j] = SVector{2}(0.0,0.0)
                end
            end
        end
    end
    # Final output flow
    output_flow = 0.5*guess_flow
    return output_flow, status_array, error_array
end

function in_image(img::AbstractArray{T, 2}, point::SVector{2, U}, window::Int) where {T <: Gray, U <: Union{Int, Float64}}
    if all(x->isa(x, Base.OneTo), indices(img))
        if point[1] < 1 || point[2] < 1 || point[1] > indices(img)[1].stop || point[2] > indices(img)[2].stop
            return false
        end
    else
        if point[1] < indices(img)[1].start || point[2] < indices(img)[2].start || point[1] > indices(img)[1].stop || point[2] > indices(img)[2].stop
            return false
        end
    end
    return true
end

function is_lost(img::AbstractArray{T, 2}, point::SVector{2, U}, min_eigen::Float64, min_eigen_thresh::Float64, window::Int) where {T <: Gray, U <: Union{Int, Float64}}
    if !(in_image(img, point, window))
        return true
    else
        val = min_eigen/window

        if val < min_eigen_thresh
            return true
        else
            return false
        end
    end
end

function is_lost(point::SVector{2, Float64}, window_size::Int)
    if point[1] > window_size || point[2] > window_size
        return true
    else
        return false
    end
end

function lies_in(area::Tuple{UnitRange{Int},UnitRange{Int}}, point::SVector{2, T}) where T <: Union{Int,Float64}
    if area[1].start <= point[1] && area[1].stop >= point[1] && area[2].start <= point[2] && area[2].stop >= point[2]
        return true
    else
        return false
    end
end

function get_grid(img::AbstractArray{T, 2}, grid::AbstractArray{UnitRange{Int}, 1}, point::SVector{2, U}, diff_flow::SVector{2, Float64}, window_size::Int) where {T <: Gray, U <: Union{Int, Float64}}

    if all(x->isa(x, Base.OneTo), indices(img))
        allowed_area = map(i -> 1+window_size:i.stop-window_size, indices(img))
    else
        allowed_area = map(i -> i.start+window_size:i.stop-window_size, indices(img))
    end

    if !lies_in(allowed_area, point + diff_flow)
        new_point = point + diff_flow

        if all(x->isa(x, Base.OneTo), indices(img))
            x21 = clamp(new_point[2]-window_size, 1, indices(img)[2].stop)
            x22 = clamp(new_point[2]+window_size, 1, indices(img)[2].stop)
        else
            x21 = clamp(new_point[2]-window_size, indices(img)[2].start, indices(img)[2].stop)
            x22 = clamp(new_point[2]+window_size, indices(img)[2].start, indices(img)[2].stop)
        end

        #TODO: Handle case for upper bound
        if x21 == 1
            x22 = ceil(x22)
            x11 = convert(Int, x21)
            x12 = convert(Int, x22)
        else
            x11 = point[2] - window_size
            x12 = point[2] + window_size
        end

        if x11 < 1
            x21 += (1 - x11)
            x11 = 1
        end

        if all(x->isa(x, Base.OneTo), indices(img))
            y21 = clamp(new_point[1]-window_size, 1, indices(img)[1].stop)
            y22 = clamp(new_point[1]+window_size, 1, indices(img)[1].stop)
        else
            y21 = clamp(new_point[1]-window_size, indices(img)[1].start, indices(img)[1].stop)
            y22 = clamp(new_point[1]+window_size, indices(img)[1].start, indices(img)[1].stop)
        end

        if y21 == 1
            y22 = ceil(y22)
            y11 = Int(y21)
            y12 = Int(y22)
        else
            y11 = point[1] - window_size
            y12 = point[1] + window_size
        end

        if y11 < 1
            y21 += (1 - y11)
            y11 = 1
        end

        grid_1 = SVector{2}(y11:y12, x11:x12)
        grid_2 = SVector{2}(y21:y22, x21:x22)
        flag = false

    else
         grid_1 = grid
         grid_2 = grid + diff_flow
         flag = true
    end
    return grid_1, grid_2, flag
end