text:
# GeneBracket
# = reduced Bracket for genetic programming
# no symbols
# no ESC, VAL, META, PRINT
# builtin for side_effects
# no need for GC during single run
# numbers as int8

# new design of bracket for julia
# closely follow the go-branch

const CELLS     = 24*1024*1024
const GCMARGIN  = CELLS - 24
const STACKSIZE = 1024*1024

# Tagbits (from right  to left)
# three bits are used (from Bit 1 to Bit 3), Bit 4 is free and can be used for gc for tree traversals
# local pointer, global pointer, Int, Prim, Symbol, Float
# global pointers not used yet, will become meta gene pool of gene bracket
# Bit 1 = 0 ->  Cell
#    Bit 2 = 0 --> pointer to cell on local heap
#    Bit 2 = 1 --> pointer to cell on global heap  (not used in the moment)
#    Bit 3 = 0 --> cons (ie, list or quotation)
#    Bit 3 = 1 --> closure
# Bit 1 = 1 ->  Number or Symb
# Bit 2 = 0 --> Symb
#    Bit 3 = 0 --> assignable symbol
#    Bit 3 = 1 --> primitive
# Bit 2 = 1 ->  Number
#    Bit 3 = 0 --> Int
#    Bit 3 = 1 --> Float

const tagType    = 7  # mask with bits 111
const tagGlobal  = 2  # bits 010    cell on global heap
const tagCell    = 1  # bits 001
const tagCons    = 5  # bits 101
const tagClosure = 4  # bits 100
const tagPrim    = 5  # bits 101
const tagSymb    = 1  # bits 001
const tagNumb    = 2  # bits 010
const tagInt     = 3  # bits 011
const tagFloat   = 7  # bits 111

boxCons(x) = x<<4    # create a new local cons
boxClosure(x) = x<<4 | tagClosure   # create a new local closure
#func boxGlobal(x) = x<<4 | tagGlobal
boxPrim(x) = x<<4 | tagPrim  # create a local primitive
boxSymb(x) = x<<4 | tagSymb
boxInt(x)  = Int64(x)<<4 | tagInt

unbox(x) = x>>4   # remove all tags
unbox8(x) = Int8(x>>4)   # remove all tags
# in contrast to C, here the pointer is
# just the heap index, that is, a number
#func ptr(x value) int    {return int(x)>>4}

isInt(x)    = x & tagType == tagInt
isFloat(x)  = x & tagType == tagFloat
isPrim(x)   = x & tagType == tagPrim
isSymb(x)   = x & tagType == tagSymb
isLocal(x)  = x & tagGlobal == 0
isGlobal(x) = x & tagGlobal == tagGlobal
isCell(x)   = x & tagCell == 0
isAtom(x)   = x & tagCell == tagCell
isCons(x)   = x & tagCons == 0
isClosure(x) = x & tagCons == tagClosure
isNumb(x)   = x & tagNumb == tagNumb
isAbstractSymb(x) = x & tagNumb == 0  # symbol or primitive

isNil(x) = x == NIL
isDef(x) = x != NIL
isCell2(vm,x) = isCell(x) && isCell(cdr(vm,x))
#isCons2(vm,x) = isCons(x) && isCons(cdr(vm,x))
#isCons3(x,vm) = isCons(x) && isCons(cdr(x,vm)) && isCons(cddr(x,vm))

struct Cell
    car :: Int
    cdr :: Int
end

# use a closure to generate primitives
function make_primitives()
   counter = 0
   function closure()
      counter += 1
      boxPrim(counter)
   end
end

newprimitive = make_primitives()
const NIL  = newprimitive()
const DUP  = newprimitive()
const DROP = newprimitive()
const SWAP = newprimitive()
const ROT  = newprimitive()
const CONS = newprimitive()
const CAR  = newprimitive()
const CDR  = newprimitive()
const ADD  = newprimitive()
const SUB  = newprimitive()
const MUL  = newprimitive()
const DIV  = newprimitive()
const LT   = newprimitive()
const GT   = newprimitive()
const RND  = newprimitive()
const EQ   = newprimitive()
const IF   = newprimitive()
const EVAL = newprimitive()
#const COND = newprimitive()
const DIP  = newprimitive()
const VAL  = newprimitive()
const REC  = newprimitive()
const DEF  = newprimitive()
const LAMBDA = newprimitive()
#const ESC   = newprimitive()
#const VESC  = newprimitive()
#const META  = newprimitive()
const TRACE = newprimitive()
const TYP   = newprimitive()
const SE    = newprimitive()  # side effects
#const PRINT = newprimitive()
# probably not really needed SET, WHL,COND, RTo, RIS
const UNBOUND = 0

const symboltable = Dict(
   "dup"=>DUP, "drop"=>DROP, "swap"=>SWAP, "rot"=>ROT, "cons"=>CONS, "car"=>CAR, "cdr"=>CDR,
   "add"=>ADD, "sub"=>SUB, "mul"=>MUL, "div"=>DIV, "+"=>ADD, "-"=>SUB,
   "*"=>MUL, "/"=>DIV, "lt"=>LT, "gt"=>GT, "eq"=>EQ,
   "rnd" =>RND, "if"=>IF, "eval"=>EVAL, "dip"=>DIP, "val"=>VAL,
   "rec"=>REC, "def"=>DEF, "lambda"=>LAMBDA, "\\"=>LAMBDA,
   "trace"=>TRACE, "typ"=>TYP, "se"=>SE)

const symboltable1 = Dict(
   "u"=>DUP, "o"=>DROP, "s"=>SWAP, "t"=>ROT,
   "n"=>CONS, "a"=>CAR, "d"=>CDR,
   "+"=>ADD, "-"=>SUB, "*"=>MUL, "/"=>DIV, "lt"=>LT, "gt"=>GT, "eq"=>EQ,
   "rnd" =>RND, "?"=>IF, "e"=>EVAL, "i"=>DIP, "v"=>VAL,
   "r"=>REC, "d"=>DEF, "l"=>LAMBDA, "\\"=>LAMBDA,
   "trace"=>TRACE, "typ"=>TYP, "se"=>SE)

#   "dip"=>DIP, "cond"=>COND,
#   "`"=>VESC, "vesc"=>VESC,

mutable struct Stats   # some statistics about the running program
    nInst :: Int    # number of Instructions
    nRecur :: Int   # recursion depth
    nSteps :: Int   # no of performed programming steps
    extent :: Int   # exent of genome at birth (vm size)
end

# virtual machine
mutable struct Vm
    arena    :: Vector{Cell}  # memory arena to hold the cells
    brena    :: Vector{Cell}  # second arena, needed for copying gc
    next     :: Int    # index into current cell in arena
    bra      :: Int    # global program stack
    ket      :: Int    # global data stack
    env      :: Int    # environment
    stack    :: Vector{Int}
    stackindex :: Int
    stats    :: Stats # some statistics about the running program
    trace    :: Int   # trace mode e: 0=no trace, 1=trace non-verbose, 3=verbose
    depth    :: Int   # current recursion depth
    need_gc  :: Bool  # flag to indicate that heap space gets rare
    side_effect
end

function Vm(side_eff)
    arena   = Vector{Cell}(undef,CELLS)
    brena   = Vector{Cell}(undef,CELLS)
    next    = 0
    stack   = Vector{Int}(undef,STACKSIZE)
    stats = Stats(0,0,0,0)
    vm = Vm(arena,brena,next,NIL,NIL,NIL,stack,0,stats,0,0,false,side_eff)
    vm.env  = cons(vm,NIL,NIL)
    vm
end

function reset!(vm)
    vm.next = 0
    vm.stats = Stats(0,0,0,0)
    vm.bra = NIL
    vm.ket = NIL
    vm.env = cons(vm,NIL,NIL)
    vm.stackindex = 0
    vm.depth = 0
    vm.trace = 0
    vm.need_gc = false
end

#  garbage collector  *********************************
#  implement Cheney copying algorithm
#    Cheney :  non-recursive traversal of live-objects
function relocate!(vm, c)
    if !isCell(c)
        return c
    end
    indb = unbox(c)    # index into brena
    @inbounds bcell = vm.brena[indb]
    if bcell.car == UNBOUND
        return bcell.cdr
    end
    inda = vm.next     # index into arena
    if isCons(c)
        c1 = boxCons(inda)
    else
        c1 = boxClosure(inda)
    end
    @inbounds begin
         vm.arena[inda] = bcell
         vm.brena[indb] = Cell(UNBOUND, c1)
    end
    vm.next += 1
    c1
end

function gc(vm)
    println("starting gc ************************************************")

    (vm.brena, vm.arena) = (vm.arena, vm.brena)
    finger =  1
    vm.next = 1

    # scan root of every live object
    vm.bra = relocate!(vm, vm.bra)
    vm.ket = relocate!(vm, vm.ket)
    vm.env = relocate!(vm, vm.env)
    @inbounds for i = 1 : vm.stackindex
      vm.stack[i] = relocate!(vm,vm.stack[i])
    end

    # scan remaining objects in heap (including objects added by this loop)
    @inbounds while finger < vm.next
       c = vm.arena[finger]
       vm.arena[finger] = Cell(relocate!(vm,c.car), relocate!(vm,c.cdr))
       finger += 1
    end

    #println("GC: live objects found: ", vm.next-1)
    #println("GC: stack ", vm.stackindex, " ", vm.depth)

    if vm.next >= GCMARGIN
        error("Bracket GC, heap too small")
    end
    vm.need_gc = false
    #println("GC finished")
end

# **********************

function new_cons(vm,pcar, pcdr)
    vm.next += 1
    if vm.next > GCMARGIN
      vm.need_gc = true
    end
    @inbounds vm.arena[vm.next] = Cell(pcar,pcdr)
    vm.next  # return a boxed index
end

cons(vm,pcar,pcdr) = boxCons(new_cons(vm,pcar,pcdr))
closure(vm,pcar,pcdr) = boxClosure(new_cons(vm,pcar,pcdr))

strip_closure(vm,cl) = isClosure(cl) ? car(vm,cl) : cl

# ------- careful that we do not leak mutablilty
#         should be used only for environments
# modify car or cdr of a cell without allocating a new cell
# should only be used for bindings
@inline function setcar!(vm, cl, newcar)
    ind = unbox(cl)
    @inbounds  c = vm.arena[ind]
    @inbounds  vm.arena[ind] = Cell(newcar, c.cdr)
end

@inline function setcdr!(vm, cl, newcdr)
    ind = unbox(cl)
    @inbounds  c = vm.arena[ind]
    @inbounds  vm.arena[ind] = Cell(c.car, newcdr)
end
# --------------------

#unsafe, assumes c is a Cell
@inbounds car(vm,c) = vm.arena[unbox(c)].car
@inbounds cdr(vm,c) = vm.arena[unbox(c)].cdr
caar(vm,c) = car(vm,car(vm,c))
cadr(vm,c) = car(vm,cdr(vm,c))
cddr(vm,c) = cdr(vm,cdr(vm,c))
pop(vm,list) = (car(vm,list), cdr(vm,list))   # unsafe
pop2(vm,list) = (car(vm,list), car(vm,cdr(vm,list)), cdr(vm,cdr(vm,list)))  # unsafe
#popsafe(vm,elem) = isCons(elem) ? pop(vm,elem) : (elem,elem)

# just count the number of conses, ie dotted pair has length 1
function length_list(vm,l)
    n=0
    while isCell(l)
        n += 1
        l = cdr(vm,l)
    end
    n
 end

 #list length, without quoted values (but also including dotted pairs)
 function length_nonquoted(vm, list)
    n = 0
    while isCell(list)
        elem = car(vm,list)
        #if elem == ESC     # a quoted element
        #   if isDef(cdr(vm,list))
        #      list = cdr(vm,list)
        #   end
        #else
        #if elem != VESC && elem != LAMBDA    #no quoted element
        if elem != VAL && elem != LAMBDA    #no quoted element
           n += 1
        end
        list = cdr(vm,list)
    end
    if isDef(list)   #count last element in dotted pair
        n += 1
    end
    n
 end

# reverse a list
# if list contained a dotted pair, reverse returns normal list
# but also a flag
# reverse only to first occurence of a closure, because
#   closures can occur only at end of a list
#   to avoid infinite loop when printing environments
function reverse_list_long(vm,list)
    l = NIL
    while isCons(list)    # take care not to pop from a closure
      p, list = pop(vm,list)
      l = cons(vm,p,l)
      if vm.need_gc
        #pushstack(vm,l)
        #pushstack(vm,list)
        #gc(vm)
        throw(BracketException("need gc"))
        #list = popstack(vm)
        #l = popstack(vm)
      end
    end
    if isClosure(list)    # take only quotation from closure, not the environment
         # seems to be wrong, no reversal!
      l = cons(vm,car(vm,list),l)
      return l, true
    elseif isDef(list)    # list contained a dotted pair
      l = cons(vm,list,l)
      return l, true
    else
      return l, false     # list did not contain a dotted pair
    end
end

# simple version if list does not contain a dotted pair
# still might be issues of list is a closure
function reverse_list(vm,list)
    l = NIL
    #list = strip_closure(vm,list)  # no need because list is not a closure
    while isCons(list)    # take care not to pop from a closure
      p, list = pop(vm,list)
      l = cons(vm,p,l)
      if vm.need_gc
        #pushstack(vm,l)
        #pushstack(vm,list)
        #gc(vm)
        throw(BracketException("need gc"))
        #list = popstack(vm)
        #l = popstack(vm)
      end
    end
    l
end

function isEqual(vm,p1,p2)
    p1 = strip_closure(vm,p1)
    p2 = strip_closure(vm,p2)
    if isCell(p1) && isCell(p2)
        isEqual(vm,car(vm,p1),car(vm,p2)) &&
        isEqual(vm,cdr(vm,p1),cdr(vm,p2))
    else
        p1 == p2
    end
end

# stack functions #############
@inline function pushstack(vm, x)
    vm.stackindex += 1
    if vm.stackindex > STACKSIZE
        error("VM stack overflow")
    end
    @inbounds vm.stack[vm.stackindex] = x
end

@inline function popstack(vm)
    #println("pop")
    #if vm.stackindex == 0
    #  error("VM stack underflow")
    #end
    @inbounds x = vm.stack[vm.stackindex]
    vm.stackindex -= 1
    return x
end

@inline function getstack(vm)
    #if vm.stackindex == 0
    #  error("VM stack underflow")
    #end
    @inbounds x = vm.stack[vm.stackindex]
    x
end

@inline function replacestack(vm, x)
    @inbounds vm.stack[vm.stackindex] = x
end

function print_stack(vm)
    println("Stack: ")
    for i = 1 : vm.stackindex
        printElem(vm, vm.stack[i])
        println()
    end
    println()
end

# creates new empty environment
@inline newenv(vm,env) = cons(vm, NIL, env)

# ************ bindings ************************************************

function find_localkey(vm, key, env)
# search binding with key in current (= top of env) frame
    if isNil(key)
        return NIL
    end
    bnds = car(vm,env)  #current frame (list of bindings) is on top of env
    while isCell(bnds)
       bnd = car(vm,bnds)
       if car(vm,bnd) == key
           return bnd
       end
       bnds = cdr(vm,bnds)
    end
    return NIL
end

function findkey(vm, key)
# search binding with key in whole environment
    if isNil(key)
        return NIL
    end
    env = vm.env
    while isDef(env)
       bnd = find_localkey(vm, key, env)
       if isDef(bnd)
          return bnd
       end
       env = cdr(vm,env)
    end
    return NIL
end

function boundvalue(vm, key) # lookup symbol..
    bnd = findkey(vm, key)
    if bnd == NIL
        return NIL
    else
        return cdr(vm,bnd)
    end
end

function bindkey(vm, key, val)
# search for key in top frame, if key found override
# otherwise make new binding in top frame
    env = vm.env
    bnd = find_localkey(vm,key,env)
    if isNil(bnd)  # key does not yet exist
       bnd = cons(vm, key, val)
       setcar!(vm, env, cons(vm,bnd,car(vm,env)))
    else           # key exists, just override val
       setcdr!(vm,bnd,val)
    end
end

function setkey(vm, key, val)
# search for key in full environment, if key found override
# otherwise make new binding in top frame
    env = vm.env
    bnd = findkey(vm,key)
    if isNil(bnd)  # key does not yet exist
       bnd = cons(vm, key, val)
       setcar!(vm, env, cons(vm,bnd,car(vm,env)))
    else           # key exists, just override val
       setcdr!(vm,bnd,val)
    end
end

#istrue(l) = isDef(l) ?  (unbox(l) != 0) : false

@inline istrue(l) = l != NIL && unbox(l) != 0
@inline isfalse(l) = l == NIL || unbox(l) == 0

# ************ io   ****************************************************+

# compared to Base64 we place the digits at the beginning
# and use 'minus' and 'underscore' as additional chars
# 0..9 have position 0..9, 'A' .. 'Z' have position 10..35
# 'a' .. 'z' have position 36..61, '-' has position 62 and '_' has 63*/
const base64_enc_table =
	"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz-_"

# position in table is ASCII index, value is the index in Base64
const base64_dec_dict = Dict(
'0'=>0, '1'=>1, '2'=>2, '3'=>3, '4'=>4, '5'=>5, '6'=>6, '7'=>7, '8'=>8, '9'=>9,
'A'=>10,'B'=>11,'C'=>12,'D'=>13,'E'=>14,'F'=>15,'G'=>16,'H'=>17,'I'=>18,'J'=>19,
'K'=>20,'L'=>21,'M'=>22,'N'=>23,'O'=>24,'P'=>25,'Q'=>26,'R'=>27,'S'=>28,'T'=>29,
'U'=>30,'V'=>31,'W'=>32,'X'=>33,'Y'=>34,'Z'=>35,'a'=>36,'b'=>37,'c'=>38,'d'=>39,
'e'=>40,'f'=>41,'g'=>42,'h'=>43,'i'=>44,'j'=>45,'k'=>46,'l'=>47,'m'=>48,'n'=>49,
'o'=>50,'p'=>51,'q'=>52,'r'=>53,'s'=>54,'t'=>55,'u'=>56,'v'=>57,'w'=>58,'x'=>59,
'y'=>60,'z'=>61,'-'=>62,'_'=>63 )

#=
function string2symbol(str)
# encode each character into 6 bits Base64-value
# only 10*6=60 bits are used,
# the remaining 4 bits can be used either as flags bits (at the right side)
   x = 0
   nok = 0
   for j=1:lastindex(str)
       if haskey(base64_dec_dict, str[j])
         x  = ( x << 6 ) | base64_dec_dict[str[j]]
         nok += 1
         if nok == 10; break; end
       end
   end
   boxSymb(x)
end
=#

#=
function symbol2string(symb)
# decode symbol back to string for output
   x = unbox(symb)  # remove the flag bits
   s = base64_enc_table[(x & 0x3F)+1]
   for i = 2 : 10
       x = x>>6
       s = s * base64_enc_table[(x & 0x3F)+1]
   end
   s = reverse(s) # revert ..
   for i=1:10     # ..  and cut-off trailing 0's
      if s[i] != '0'
        return s[i:10]
      end
   end
   return ""
end
=#

function printElem(vm, q)
    if isInt(q)
        print(unbox(q))
    #elseif isFloat(q)
    #    print(unbox_float(q))
    elseif isNil(q)
         print("[]")
    elseif isPrim(q)
         print(first(keys(filter(p->p.second == q , symboltable))))
    #elseif isSymb(q)
    #     print(symbol2string(q))
    else
        printList(vm,q)
    end
end

function printInnerList(vm, list, invert)
    isdotted = false
    list = strip_closure(vm,list)
    if isCell(list)
       if invert
          list, isdotted = reverse_list_long(vm,list)
       end
       p, list = pop(vm,list)
       printElem(vm,p)
       if isdotted   # dotted list that was reversed
          print(" .")
       end
       while isCell(list)
          p, list = pop(vm,list)
          print(" ")
          printElem(vm,p)
       end
       if isDef(list)   # dotted list
          print(" *")
          printElem(vm,list)
       end
    end
end

function printList(vm,l)
    print("[")
    printInnerList(vm,l,true)
    print("]")
end

function printKet(vm,l)
    print("[")
    printInnerList(vm,l,false)
    println(">")
end

function printBra(vm,l)
    print("<")
    printInnerList(vm,l,true)
    println("]")
end

function atom(token)
     if (length(token) == 1) && haskey(symboltable,token)   # + - * /
         return symboltable[token]
     end
     sign = 1
     if token[1] == '+'; token = token[2:end]; end
     if token[1] == '-'
       token = token[2:end]
       sign = -1
     end

     #ind = findfirst(x -> x=='.', token)
     #if ind != nothing && all(isdigit,token[1:ind-1]) && all(isdigit,token[ind+1:end])
     #      f = Float32(sign*Meta.parse(token))
     #      a = box_float(f)
     #else
     if all(isdigit,token)   # Ints
           a = boxInt(sign*Meta.parse(token))
     elseif haskey(symboltable,token)
           a = symboltable[token]
     else
         error("Symbol in code")
          # a = string2symbol(token)
     end
     a
end

function nextchar(io)
    c = '\0'
    while true
       if eof(io); return '\0'; end
       c = read(io, Char)
       if c == ';'  # single line comment
          while true
             if eof(io); return '\0'; end
             c = read(io, Char)
             if c == '\n'; break; end
          end
       end
       if !isspace(c); break; end
    end
    return c
end

function read_token(io)
    token = ""
    ndots = 0  # number of dots in token
    while true
       if eof(io); break; end
       c = read(io, Char)
       if c == '.'
         ndots += 1
       end
       if ndots > 1 # can have at most 1 dot in number
          skip(io,-1)
          break
       end
       if c == '\n' || isspace(c)
          break
       elseif !(isdigit(c) || c=='+' || c=='-' || c=='.' || c=='_'
                          || c == '*' || c == '/' || isletter(c))
          #write(io,c)
          skip(io,-1)
          break
       elseif c=='`' || c == '\'' || c == '\\' || c == ';'
          #write(io,c)
          skip(io,-1)
          break
       else
         #push!(buf,c)
         token = token * c
       end
    end
    token
end

function read_tokens!(vm,io)
    val = NIL
    while true
      c = nextchar(io)
      if c == '\0'       # end of stream
        return val
      elseif c == ']'    # end of list
        return val
      elseif c == '['    # begin of list
        newval = read_tokens!(vm,io)
        val = cons(vm,newval, val)
      #elseif c == '\''   # escape
      #  val = cons(vm,ESC, val)
      #elseif c == '`'   # Vesc (escape value)
      #  val = cons(vm, VESC, val)
      elseif c == '\\'   # Backslash = lambda
        val = cons(vm, LAMBDA, val)
      else               # read new atom
        #write(io, c)
        skip(io,-1)
        token = read_token(io)
        newval = atom(token)
        val = cons(vm, newval, val)
      end
    end
end

function make_bra(vm,prog::String)
    io = IOBuffer(prog)
    read_tokens!(vm, io)
end

function load_file(vm, prog)
     io = open(prog)
     val = read_tokens!(vm, io)
     close(io)
     val
end

# **************************** builtins *********************************
function f_dup(vm)
    if isCell(vm.ket)
        vm.ket = cons(vm,car(vm,vm.ket), vm.ket)
    end
 end

function f_drop(vm)
    if isCell(vm.ket)
        vm.ket = cdr(vm,vm.ket)
    end
end

function f_swap(vm)
    if isCell2(vm,vm.ket)
        a,b,vm.ket  = pop2(vm,vm.ket)
        vm.ket = cons(vm,a,vm.ket)
        vm.ket = cons(vm,b,vm.ket)
    end
end

function f_rot(vm)
    if isCell(vm.ket)
        a, vm.ket = pop(vm,vm.ket)
        if isCell(vm.ket)
            b,vm.ket = pop(vm,vm.ket)
            if isCell(vm.ket)
              c,vm.ket = pop(vm,vm.ket)
              vm.ket = cons(vm,b,vm.ket)
              vm.ket = cons(vm,a,vm.ket)
              vm.ket = cons(vm,c,vm.ket)
            end
        end
    end
end

function f_cons(vm)
    if isCell2(vm,vm.ket)
       p1, p2, vm.ket = pop2(vm,vm.ket)
       p2 = strip_closure(vm, p2)  # cons to a closure strips the closure
                  # instead here we could cons to the quotation of the closure
       vm.ket = cons(vm,cons(vm,p1,p2),vm.ket)
    end
end

function f_car(vm)
    if isCell(vm.ket)
        p, vm.ket = pop(vm,vm.ket)
        p = strip_closure(vm, p)
        if isCell(p) # car a list
            head,p = pop(vm,p)
            # !!!!!!
            #vm.ket = cons(vm,p,vm.ket)  # leave rest of list on the ket
            vm.ket = cons(vm,head,vm.ket)
        else          # car a symbol
            ####if isNil(p); return; end
            val = boundvalue(vm,p)  # lookup symbol
            val = strip_closure(vm, val)
            if isCons(val)
              vm.ket = cons(vm,car(vm,val), vm.ket)
            end
        end
    end
end

function f_cdr(vm)
    if isCell(vm.ket)
        p, vm.ket = pop(vm,vm.ket)
        p = strip_closure(vm, p)
        if isCell(p)   # cdr a list
            head, p = pop(vm,p)
            vm.ket = cons(vm,p,vm.ket)
        else
            val = boundvalue(vm,p)  # look up symbole
            val = strip_closure(vm, val)
            if isCons(val)
               vm.ket = cons(vm,cdr(vm,val), vm.ket)
            else
               vm.ket = cons(vm,NIL,vm.ket) # at least leave a nill on ket
            end
        end
    end
end

# some math functions (we can easily extend to more..)
rnd(x) = x < 1 ? 0 : rand(1:x)
lt(x1,x2) = x1 < x2 ? 1 : 0
gt(x1,x2) = x1 > x2 ? 1 : 0
my_div(x1,x2) = x2 == 0 ? 0 : div(x1,x2)

function f_math1(vm,op)
    if isCell2(vm,vm.ket)
       n1,n2, vm.ket = pop2(vm,vm.ket)
       #if isSymb(n1)
       #  n1 = boundvalue(vm,n1)
       #end
       #if isSymb(n2)
       #  n2 = boundvalue(vm,n2)
       #end
       # Here we should check that n1 and n2 are numbers now !!!!!!
       n3 = boxInt(op(unbox8(n1), unbox8(n2)))
       vm.ket = cons(vm,n3,vm.ket)
    end
end

function f_math(vm,op)
    if isCell2(vm,vm.ket)
       n1,n2, vm.ket = pop2(vm,vm.ket)
       #if isSymb(n1)
       #  n1 = boundvalue(vm,n1)
       #end
       #if isSymb(n2)
       #  n2 = boundvalue(vm,n2)
       #end
       if isNumb(n1) && isNumb(n2)
          n3 = boxInt(op(unbox8(n1), unbox8(n2)))
          vm.ket = cons(vm, n3, vm.ket)
       elseif isCell(n1) && isCell(n2)
          n1 = strip_closure(vm,n1)
          n2 = strip_closure(vm,n2)
          c = NIL
          while isCell(n1) && isCell(n2)
             c1, n1 = pop(vm,n1)
             c2, n2 = pop(vm,n2)
             #if isSymb(c1)
             #    c1 = boundvalue(vm,c1)
             #end
             #if isSymb(c2)
             #    c2 = boundvalue(vm,c2)
             #end
             if isNumb(c1) && isNumb(c2)
                 n3 = boxInt(op(unbox8(c1), unbox8(c2)))
                 c = cons(vm, n3, c)
             end
             if vm.need_gc
                #pushstack(vm,c)
                #pushstack(vm,n1)
                #pushstack(vm,n2)
                #gc(vm)
                throw(BracketException("need gc"))
                #n2 = popstack(vm)
                #n1 = popstack(vm)
                #c  = popstack(vm)
             end
          end
          #c,_ =reverse_list(vm,c)
          c =reverse_list(vm,c)
          vm.ket = cons(vm,c,vm.ket)
       elseif isCell(n1)
          n1 = strip_closure(vm,n1)
          c = NIL
          while isCell(n1)
             c1, n1 = pop(vm,n1)
             #if isSymb(c1)
             #    c1 = boundvalue(vm,c1)
             #end
             if isNumb(c1) && isNumb(n2)
                 n3 = boxInt(op(unbox8(c1), unbox8(n2)))
                 c = cons(vm, n3, c)
             end
             if vm.need_gc
                #pushstack(vm,c)
                #pushstack(vm,n1)
                #gc(vm)
                throw(BracketException("need gc"))
                #n1 = popstack(vm)
                #c  = popstack(vm)
             end
          end
          #c,_ =reverse_list(vm,c)
          c =reverse_list(vm,c)
          vm.ket = cons(vm,c,vm.ket)
       elseif isCell(n2)
          n2 = strip_closure(vm,n2)
          c = NIL
          while isCell(n2)
             c2, n2 = pop(vm,n2)
             #if isSymb(c2)
             #    c2 = boundvalue(vm,c2)
             #end
             if isNumb(n1) && isNumb(c2)
                 n3 = boxInt(op(unbox8(n1), unbox8(c2)))
                 c = cons(vm, n3, c)
             end
             if vm.need_gc
                #pushstack(vm,c)
                #pushstack(vm,n2)
                #gc(vm)
                throw(BracketException("need gc"))
                #n2 = popstack(vm)
                #c  = popstack(vm)
             end
          end
          #c,_ =reverse_list(vm,c)
          c =reverse_list(vm,c)
          vm.ket = cons(vm,c,vm.ket)
       end
    end
end

function f_rnd(vm)
    if isCell(vm.ket)
        p, vm.ket = pop(vm,vm.ket)
        if isInt(p)
            if p >= 1
               p=boxInt(rand(1:unbox8(p)))
            else
               p=boxInt(0)
            end
        elseif isCell(p)
            p = strip_closure(vm,p)
            n=length_list(vm,p)-1
            n1=rand(0:n)
            for i=1:n1
                p = cdr(vm,p)
            end
            p = car(vm,p)
        end
        vm.ket=cons(vm,p,vm.ket)
    end
end

function f_eq(vm)
    if isCell2(vm,vm.ket)
       p1, p2, vm.ket = pop2(vm,vm.ket)
       b = isEqual(vm,p1,p2)
       # <<4 of Boolean is automatically converted to Int64
       vm.ket = cons(vm,boxInt(b),vm.ket)
    end
end

function f_if(vm)
    if isCell(vm.ket)
        b, vm.ket = pop(vm,vm.ket)
        if isCell(vm.ket)
            e1,vm.ket = pop(vm,vm.ket)
            if isCell(vm.ket)
              e2,vm.ket = pop(vm,vm.ket)
              vm.ket = istrue(b) ? cons(vm,e1,vm.ket) : cons(vm,e2,vm.ket)
            end
        end
    end
end

function f_dip(vm)
    if isCell2(vm,vm.ket)
        q1,q2,vm.ket = pop2(vm,vm.ket)
        vm.bra = cons(vm,q2,vm.bra)
        vm.bra = cons(vm,EVAL,vm.bra)
        vm.ket = cons(vm,q1,vm.ket)
    end
end

function f_esc(vm)
    if isCell(vm.bra)
        val, vm.bra = pop(vm,vm.bra)
        vm.ket = cons(vm,val,vm.ket)
     end
end

#=
function f_vesc(vm)
    if isCell(vm.bra)
        val, vm.bra = pop(vm,vm.bra)
        vm.ket = cons(vm,val,vm.ket)
        f_val(vm)
     end
end
=#

function f_val(vm)
    if isCell(vm.ket)
        key, vm.ket = pop(vm,vm.ket)
        if isCell(key)
           vm.ket = cons(vm,key,vm.ket)
        else
           val = boundvalue(vm,key)       # lookup symbol ..
           vm.ket = cons(vm,val,vm.ket)   # .. and place on ket
        end
     end
end

function f_typ(vm)
    if isCell(vm.ket)
       p,vm.ket = pop(vm,vm.ket)
       if isInt(p)
          t = 1
       elseif isPrim(p)
          t = 2
       #elseif isSymb(p)
       #   t=3
       elseif isCons(p)
          t=4
       elseif isClosure(p)
          t=5
       else
          t=0
       end
       vm.ket = cons(vm,boxInt(t),vm.ket)
    end
end

function f_trace(vm)
    if isCell(vm.ket)
       p,vm.ket = pop(vm,vm.ket)
       vm.trace = unbox(p)
    end
end

function f_print(vm)
    if isCell(vm.ket)
        p,vm.ket = pop(vm,vm.ket)
        #print_bra(p,vm)
        #println(" here goes the next")
        printElem(vm,p)
        print(" ")
    end
end

function f_rec(vm)
    #anonymous recursion: replace bra of this scope by original value
    if isCell(vm.ket)
        b, vm.ket = pop(vm,vm.ket)
        if istrue(b)
            vm.bra = getstack(vm)
        end
    end
end

function f_lambda(vm)
   if isCell2(vm,vm.ket)
       keys,q,vm.ket = pop2(vm,vm.ket)
       if isAtom(q)
          q = boundvalue(vm, q)
       end
       if isCons(q) # make a closure (only of not yet)
           if isDef(keys)           # if arguments are not NIL ..
               q = cons(vm,DEF, q)  # .. push a definition on q
               q = cons(vm,keys, q)
               #if isAtom(keys)
               #    q = cons(vm,ESC, q)
               #end
           end
          env  = newenv(vm, vm.env)
          clos = closure(vm,q,env)
        elseif isClosure(q) # new closure with keys as quotation
          env  = cdr(vm,q)
          clos = closure(vm,keys,env)
        else  # isAtom(q) # we need a quotation to do lambda
           return
        end
        vm.ket = cons(vm, clos, vm.ket) # push the new closure on ket
    end
end

function deepbind(vm, keys, val)
# recursively bind all values of list keys to atom val
# keys must be a list, val an atom
    while isCell(keys)
        key, keys = pop(vm,keys)
        if isAtom(key)
            bindkey(vm,key,val)
            if vm.need_gc
                #pushstack(vm,keys)
                #gc(vm)
                throw(BracketException("need gc"))
                #keys = popstack(vm)
            end
        else   # key itself is a list
            #pushstack(vm,keys)
            deepbind(vm,key, val)
            #keys = popstack(vm)
        end
    end
end

function match(vm, keys, vals)
# bind elements from keys to elements from vals with pattern matching
# keys must be a list
    if isAtom(vals)
        deepbind(vm,keys,vals)
        return
    end                                # Q: do we need an else here??
    while isCell(keys)
       key, keys = pop(vm,keys)
       if isNil(keys)
           bindkey(vm,key,vals)
       else
           val, vals = pop(vm, vals)    # Q: do we need to check for list ??
           bindkey(vm,key,val)
       end
       if isAtom(vals)
          deepbind(vm, keys, vals)
          return
       end
    end
end

function f_def(vm)
   if isCell(vm.ket)
       key, vm.ket = pop(vm,vm.ket)
       if isAtom(key)
           if isCell(vm.ket)
              val, vm.ket = pop(vm, vm.ket)
              bindkey(vm,key,val)  #bind key to val in top env-frame
           end
       elseif isDef(key)       # binding a list of keys
           n = length_nonquoted(vm, key)
           n1 = 0    # push max n values from ket to stack
           for i = 1 : n
               if isCell(vm.ket)
                   val, vm.ket = pop(vm,vm.ket)
                   n1 += 1
                   pushstack(vm,val)
               else
                  break
               end
           end
           for i = 1 : n1 # make the bindings
               k, key = pop(vm, key)  # we should have enough elements in key
               #if k == VESC
               if k == VAL
                  k, key = pop(vm,key)   # this is interpreted as set
                  setkey(vm, k, popstack(vm))
               elseif isAtom(k)
                  bindkey(vm, k, popstack(vm))
               else
                  elem = popstack(vm)
                  pushstack(vm,key)  # safe key in case of gc
                  match(vm,k,elem)
                  key = popstack(vm)
               end
           end
       end
   end
end

# we remove set (replaced by backtick)
#=
function f_set(vm)
   if isCell2(vm, vm.ket)
       key, val, vm.ket = pop2(vm, vm.ket)
       if isAtom(key)
           setkey(vm,key,val)  # bind key to val in top env-frame
       end
    end
end
=#

function f_side(vm)
# perform a side effect
   e,vm.ket = pop(vm,vm.ket)
   vm.side_effect(unbox(e))
   #vm.ket = cons(vm, res, vm.ket)
end

function f_eval(vm)
    if isCell(vm.ket)
        op,vm.ket = pop(vm,vm.ket)
        if isCons(op)
            eval_cons(vm,op)
        elseif isClosure(op)
            eval_closure(vm,op)
        elseif isNil(op)
            return
        elseif isPrim(op)
              eval_prim(vm,op)
        #elseif isSymb(op)
        #      eval_symb(vm,op)
        else    # eval a number
              eval_numb(vm,op)
        end
    end
end

function eval_cons(vm, op)
    if isCell(vm.bra)
        vm.depth += 1
        pushstack(vm,vm.env)
        pushstack(vm,vm.bra)
        pushstack(vm,op)   # 2nd
        vm.env = newenv(vm,vm.env)
    else      # tail position
        replacestack(vm,op)
    end
    vm.bra = op
end

function eval_closure(vm, clos)
    op = car(vm,clos)
    env = newenv(vm, cdr(vm,clos))
    if isCell(vm.bra)
        vm.depth += 1
        pushstack(vm,vm.env)
        pushstack(vm,vm.bra)
        pushstack(vm,op)   # 2nd
    else      # tail position
        replacestack(vm,op)
    end
    vm.env = env
    vm.bra = op
end

function eval_numb(vm, n)
    val = boundvalue(vm,n)
    if isCons(val)
        eval_cons(vm,val)
    elseif isClosure(val)
        eval_closure(vm,val)
    else
        vm.ket = cons(vm,val,vm.ket)
    end
    #vm.ket = cons(vm,n,vm.ket)
end

#=
function eval_symb(vm,sym)
    val = boundvalue(vm,sym)
    if isCons(val)
        eval_cons(vm,val)
    elseif isClosure(val)
        eval_closure(vm,val)
    else
        vm.ket = cons(vm,val,vm.ket)
    end
end
=#

@inline function eval_prim(vm,x)
    #println("eval prim")
    if x == DUP
        f_dup(vm)
    elseif x == DROP
        f_drop(vm)
    elseif x == SWAP
        f_swap(vm)
    elseif x == ROT
        f_rot(vm)
    elseif x == CONS
        f_cons(vm)
    elseif x == CAR
        f_car(vm)
    elseif x == CDR
        f_cdr(vm)
    elseif x == ADD
        f_math(vm,+)
    elseif x == SUB
        f_math(vm,-)
    elseif x == MUL
        f_math(vm,*)
    elseif x == DIV
        f_math(vm,my_div)
    elseif x == LT
        f_math(vm,lt)
    elseif x == GT
        f_math(vm,gt)
    elseif x == RND
        f_rnd(vm)
    elseif x == EQ
        f_eq(vm)
    elseif x == IF
        f_if(vm)
    #elseif x == COND
    #    f_cond(vm)
    elseif x ==TYP
        f_typ(vm)
    elseif x == EVAL
        f_eval(vm)
    elseif x == DEF
        f_def(vm)
    elseif x == LAMBDA
        f_lambda(vm)
    elseif x == REC
        f_rec(vm)
    elseif x == VAL
        f_val(vm)
   # elseif x == ESC
   #     f_esc(vm)
   # elseif x == VESC
   #     f_vesc(vm)
    elseif x == DIP
        f_dip(vm)
    elseif x == SE
        f_side(vm)
    elseif x == TRACE
        f_trace(vm)
    #elseif x == PRINT
    #    f_print(vm)
    #elseif x == META
    #    f_meta(vm)
    end
    nothing

#      @match s begin
#         2 => f_dup!(vm)
#       #  3 => f_drop!(vm)
#       #  4 => f_swap!(vm)
#       #  5 => f_cons!(vm)
#       #  6 => f_car!(vm)
#       #  7 => f_cdr!(vm)
#         8 => f_add!(vm)
#       #  9 => f_lt!(vm)
#         10 => f_gt!(vm)
#       #  11 => f_eq!(vm)
#       #  12 => f_if!(vm)
#       #  13 => f_eval!(vm)
#       #  14 => f_each!(vm)
#         15 => f_whl!(vm)
#         _  => vm.ket = cons(e,vm.ket,vm)
#      end

end

function eval_bra(vm)
    #try

    #println("eval bra")
    if isAtom(vm.bra)
        return
    end
    starting_depth = vm.depth
    pushstack(vm,vm.bra)

    while true
      if vm.trace > 0
         #println("trace ")
         printBra(vm,vm.bra)
         printKet(vm,vm.ket)
         #println("stackdepth ", vm.stackindex)
         #printElem(vm.stack[vm.stackindex],false, vm)
         #printBra(vm.ket,vm)
         #println("depth ", vm.depth)
         #if vm.trace > 1
         #    dump_env(vm.env[vm.depth],vm)
         #end
         println()
      end

      e, vm.bra = pop(vm,vm.bra)
      if isNil(e)
          vm.ket = cons(vm,e,vm.ket)
      elseif isPrim(e)
          eval_prim(vm,e)
      #elseif isSymb(e)
      #     eval_symb(vm,e)
      else
           vm.ket = cons(vm,e,vm.ket)
      end

      if vm.need_gc
         #gc(vm)
         throw(BracketException("need gc"))

      end

      if isAtom(vm.bra)   # exit scope
         if vm.depth == starting_depth
               break
         end
         vm.depth -= 1
         popstack(vm) # 2nd
         vm.bra = popstack(vm)
         vm.env = popstack(vm)
      end
    end
    vm.bra = popstack(vm)

    #catch
    #  println("error")
    #end
 end

 # questions
 #
 # def: should this leave the result on the ket?
 #
 # car quotion: keep the cdr or not?
 #
 # do we need rec?

 # break (leave quotation if false)
 # call-cc, continuations


text:
module SQL

include("sql_exports.jl")
include("sql_imports.jl") # DBMS Structured SubQuery _to_sql _placeholder _placeholders

const DatabaseID = DBMS.SQL

"""
    to_sql(query::Structured)::String
"""
to_sql(query::Structured)::String = _to_sql(DatabaseID(), query)

"""
    to_sql(subquery::SubQuery)::String
"""
to_sql(subquery::SubQuery)::String = _to_sql(DatabaseID(), subquery)

placeholder(nth::Int) = _placeholder(DatabaseID(), nth)
placeholders(dims::Int) = _placeholders(DatabaseID(), dims)

using ...AdapterBase: AdapterBase, SQLElement

function Base.show(io::IO, mime::MIME"text/plain", element::Union{E,Structured} where E<:SQLElement)
    dbms = DatabaseID()
    AdapterBase._show(io, mime, dbms, element)
end

end # Octo.Adapters.SQL

ext: jl


text:
"""
    Dubins3D

Implementation of 2D/3D Dubins maneuver
"""
module Dubins3D

include("dubinsmaneuver3d.jl")

export DubinsManeuver2D
export DubinsManeuver3D, getLowerBound, getUpperBound, compute_sampling, sample_path, get_configuration

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# Since there is no minimal working example (MWE), I could not figure out what you really wanted to do.
# (Please read https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757)
#
# However, I thought that possibly you would like to do something like this. Please take a look at the MWE below.

# %%
module SCI

export Problem, Algorithm_A, Algorithm_B, solve

abstract type AbstractProblem end
struct Problem{T} <: AbstractProblem x::T end

abstract type AbstractAlgorithm end
Base.@kwdef struct Algorithm_A{T} <: AbstractAlgorithm a::T = 2.0 end
Base.@kwdef struct Algorithm_B{T} <: AbstractAlgorithm a::T = 3.0 end
default_algorithm(prob::Problem) = Algorithm_A()

struct Solution{R, P<:AbstractProblem, A<:AbstractAlgorithm} result::R; prob::P; alg::A end
solve(prob::AbstractProblem) = solve(prob, default_algorithm(prob))
solve(prob::AbstractProblem, alg::AbstractAlgorithm) = Solution(alg.a * prob.x, prob, alg)

"""
Here, `Base` module plays the role of `COM` package.
SCI module can define the method of
`Base.show` = `COM.show` function for `SCI.Solution` type.
"""
function Base.show(io::IO, sol::Solution)
    result = """
    Problem:   $(sol.prob)
    Algorithm: $(sol.alg)
    Result:    $(sol.result)
    """
    print(io, result)
end

end

using .SCI
@show prob = Problem(1.5)
println()
println(solve(prob))
println(solve(prob, Algorithm_B()))

# %%
module SCI2 # extension of SCI module

export Problem2, Algorithm_C

using ..SCI: SCI, AbstractProblem, AbstractAlgorithm, Solution

struct Problem2{T} <: AbstractProblem x::T; y::T end
Base.@kwdef struct Algorithm_C{T} <: AbstractAlgorithm a::T = 2.0; b = 3.0 end
SCI.default_algorithm(prob::Problem2) = Algorithm_C()
SCI.solve(prob::Problem2, alg::Algorithm_C) = Solution(alg.a * prob.x + alg.b * prob.y, prob, alg)

end

using .SCI2
@show prob2 = Problem2(10.0, 1.0)
println()
println(solve(prob2))

# %% [markdown]
# In the above example I don't define `AbstractSolution` type, nor did I make `Solution` a subtype of it. But you can define it, and then define another Solution type as a subtype of it in `SCI2` module, and also define how to display objects of another Solution type.   Many other extensions are also possible.
#
# See also https://discourse.julialang.org/t/function-depending-on-the-global-variable-inside-module/64322/10

# %%

max_stars_repo_path: 0011/Handling of package inter-dependency.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using BlockArrays, SparseArrays, Setfield
####################################################################################################
# method using the Trapezoidal rule (Order 2 in time) and discretisation of the periodic orbit.

"""
	pb = PeriodicOrbitTrapProblem(F, J, ϕ, xπ, M::Int)
This composite type implements Finite Differences based on a Trapezoidal rule to locate periodic orbits. The arguments are as follows
- `F` vector field
- `J` jacobian of `F`
- `d2F = nothing` Hessian of F (optional), useful for continuation of Fold of periodic orbits
- `ϕ` used to set a section for the phase constraint equation
- `xπ` used in the section for the phase constraint equation
- `M::Int` number of time slices
- `linsolver: = DefaultLS()` linear solver for each time slice, i.e. to solve `J⋅sol = rhs`. This is only used for the computation of the Floquet multipliers.
- `isinplace::Bool` whether `F` and `J` are inplace functions (Experimental). In this case, the functions `F` and `J` must have the following definitions `(o, x) ->  F(o, x)` and `(o, x, dx) -> J(o, x, dx)`.
- `ongpu::Bool` whether the computation takes place on the gpu (Experimental)

You can then call `pb(orbitguess)` to compute the functional on a `orbitguess`. Note that `orbitguess` must be of size M * N + 1 where N is the number of unknowns in the state space and `orbitguess[M*N+1]` is an estimate of the period of the limit cycle.

The scheme is as follows, one look for `T = x[end]` and

 ``\\left(x_{i} - x_{i-1}\\right) - \\frac{h}{2} \\left(F(x_{i}) + F(x_{i-1})\\right) = 0,\\ i=1,\\cdots,m-1``

with `u_{0} = u_{m-1}` and the periodicity condition `u_{m} - u_{1} = 0` and

where `h = T/M`. Finally, the phase of the periodic orbit is constrained by using a section

 ``\\langle x[1] - x_\\pi, \\phi\\rangle=0.``

 A functional, hereby called `G`, encodes this problem. The following methods are available

- `pb(orbitguess)` evaluates the functional G on `orbitguess`
- `pb(orbitguess, du)` evaluates the jacobian `dG(orbitguess).du` functional at `orbitguess` on `du`
- `pb(Val(:JacFullSparse), orbitguess)` return the sparse matrix of the jacobian `dG(orbitguess)` at `orbitguess` without the constraints. It is called `A_γ` in the docs.
- `pb(Val(:JacCyclicSparse), orbitguess)` return the sparse cyclic matrix Jc (see the docs) of the jacobian `dG(orbitguess)` at `orbitguess`
- `pb(Val(:BlockDiagSparse), orbitguess)` return the diagonal of the sparse matrix of the jacobian `dG(orbitguess)` at `orbitguess`. This allows to design Jacobi preconditioner. Use `blockdiag`.

!!! note "GPU call"
    For these methods to work on the GPU, for example with `CuArrays` in mode `allowscalar(false)`, we face the issue that the function `extractPeriodFDTrap` won't be well defined because it is a scalar operation. One may have to redefine it like `extractPeriodFDTrap(x::CuArray) = x[end:end]` or something else. Also, note that you must pass the option `ongpu = true` for the functional to be evaluated efficiently on the gpu.

"""
@with_kw struct PeriodicOrbitTrapProblem{TF, TJ, TJt, Td2F, vectype, Tls <: AbstractLinearSolver}
	# Function F(x, p) = 0
	F::TF

	# Jacobian of F wrt x
	J::TJ

	# Jacobian transpose of F wrt x. This is mainly used for matrix-free computation of Folds of limit cycles
	Jt::TJt = nothing

	# Hessian of F wrt x, useful for continuation of Fold of periodic orbits
	d2F::Td2F = nothing

	# variables to define a Section for the phase constraint equation
	ϕ::vectype
	xπ::vectype

	# discretisation of the time interval
	M::Int = 1

	# dimension of the problem in case of an AbstractVector
	N::Int = 0

	# linear solver for each slice, i.e. to solve J⋅sol = rhs. This is mainly used for the computation of the Floquet coefficients
	linsolver::Tls = DefaultLS()

	# whether F and J are inplace functions
	isinplace::Bool = false

	# whether the computation takes place on the gpu
	ongpu::Bool = false
end

isinplace(pb::PeriodicOrbitTrapProblem) = pb.isinplace
ongpu(pb::PeriodicOrbitTrapProblem) = pb.ongpu

function applyF(pb::PeriodicOrbitTrapProblem, dest, x)
	if isinplace(pb)
		pb.F(dest, x)
	else
		dest .= pb.F(x)
	end
	dest
end

function applyJ(pb::PeriodicOrbitTrapProblem, dest, x, dx)
	if isinplace(pb)
		pb.J(dest, x, dx)
	else
		dest .= apply(pb.J(x), dx)
	end
	dest
end

function PeriodicOrbitTrapProblem(F::TF, J::TJ, ϕ::vectype, xπ::vectype, M, ls::AbstractLinearSolver = DefaultLS(); isinplace = false, ongpu = false) where {TF, TJ, vectype}
	if ϕ isa AbstractVector
		return PeriodicOrbitTrapProblem{TF, TJ, Nothing, Nothing, vectype, typeof(ls)}(F = F, J = J, ϕ = ϕ, xπ = xπ, M = M, N = length(ϕ), linsolver = ls, isinplace = isinplace, ongpu = ongpu)
	else
		return PeriodicOrbitTrapProblem{TF, TJ, Nothing, Nothing, vectype, typeof(ls)}(F = F, J = J, ϕ = ϕ, xπ = xπ, M = M, 0, linsolver = ls, isinplace = isinplace, ongpu = ongpu)
	end
end

function PeriodicOrbitTrapProblem(F::TF, J::TJ, d2F::Td2f, ϕ::vectype, xπ::vectype, M, ls::AbstractLinearSolver = DefaultLS(); isinplace = false, ongpu = false) where {TF, TJ, Td2f, vectype}
	if ϕ isa AbstractVector
		return PeriodicOrbitTrapProblem{TF, TJ, Nothing, Td2f, vectype, typeof(ls)}(F = F, J = J, d2F = d2F, ϕ = ϕ, xπ = xπ, M = M, length(ϕ), linsolver = ls, isinplace = isinplace, ongpu = ongpu)
	else
		return PeriodicOrbitTrapProblem{TF, TJ, Nothing, Td2f, vectype, typeof(ls)}(F = F, J = J, d2F = d2F, ϕ = ϕ, xπ = xπ, M = M, 0, linsolver = ls, isinplace = isinplace, ongpu = ongpu)
	end
end

# these functions extract the last component of the periodic orbit guess
extractPeriodFDTrap(x::AbstractVector) = x[end]
extractPeriodFDTrap(x::BorderedArray)  = x.T

# these functions extract the time slices components
extractTimeSlice(x::AbstractVector, N, M) = @views reshape(x[1:end-1], N, M)
extractTimeSlice(x::BorderedArray,  N, M) = x.u

function POTrapScheme!(pb::PeriodicOrbitTrapProblem, dest, u1, u2, h, tmp, linear::Bool = true)
	# this function implements the basic implicit scheme used for the time integration
	# because this function is called in a cyclic manner, we save in tmp the value F(u2) in order to avoid recomputing it in a subsequent call
	# basically tmp is F(u2)
	if linear
		dest .= tmp
		# tmp <- pb.F(u1)
		applyF(pb, tmp, u1) #TODO this line does not almost seem to be type stable in code_wartype, gives @_11::Union{Nothing, Tuple{Int64,Int64}}
		dest .= u1 .- u2 .- h .* (dest .+ tmp)
	else
		dest .-= h .* tmp
		# tmp <- pb.F(u1)
		applyF(pb, tmp, u1)
		dest .-= h .* tmp
	end
end

function POTrapSchemeJac!(pb::PeriodicOrbitTrapProblem, dest, u1, u2, du1, du2, h, tmp)
	# this function implements the basic implicit scheme used for the time integration
	# useful for the matrix-free jacobian
	# basically tmp is dF(u2).du2 (see above for explanation)
	dest .= tmp
	# tmp <- apply(pb.J(u1), du1)
	applyJ(pb, tmp, u1, du1)
	dest .= du1 .- du2 .- h .* (dest .+ tmp)
end

"""
This function implements the functional for finding periodic orbits based on finite differences using the Trapezoidal rule. It works for inplace / out of place vector fields pb.F
"""
function POTrapFunctional!(pb::PeriodicOrbitTrapProblem, out::AbstractVector, u0::AbstractVector)
	M = pb.M
	N = pb.N
	T = extractPeriodFDTrap(u0)
	h = T / M

	u0c  = extractTimeSlice(u0, N, M)
	outc = extractTimeSlice(out, N, M)

	# outc[:, M] plays the role of tmp until it is used just after the for-loop
	@views applyF(pb, outc[:, M], u0c[:, M-1])

	@views POTrapScheme!(pb, outc[:, 1], u0c[:, 1], u0c[:, M-1], h/2, outc[:, M])

	for ii = 2:M-1
		# this function avoids computing F(u0c[:, ii]) twice
		@views POTrapScheme!(pb, outc[:, ii], u0c[:, ii], u0c[:, ii-1], h/2, outc[:, M])
	end

	# closure condition ensuring a periodic orbit
	outc[:, M] .= @views u0c[:, M] .- u0c[:, 1]

	# this is for CuArrays.jl to work in the mode allowscalar(false)
	if ongpu(pb)
		return @views vcat(out[1:end-1], dot(u0c[:, 1], pb.ϕ) - dot(pb.xπ, pb.ϕ)) # this is the phase condition
	else
		out[end] = @views dot(u0c[:, 1], pb.ϕ) - dot(pb.xπ, pb.ϕ) #dot(u0c[:, 1] .- pb.xπ, pb.ϕ)
		return out
	end
end

"""
Matrix free expression of the Jacobian of the problem for computing periodic obits when evaluated at `u0` and applied to `du`.
"""
function POTrapFunctionalJac!(pb::PeriodicOrbitTrapProblem, out, u0::AbstractVector, du)
	M = pb.M
	N = pb.N
	T  = extractPeriodFDTrap(u0)
	dT = extractPeriodFDTrap(du)
	h = T / M

	u0c = extractTimeSlice(u0, N, M)
	outc = extractTimeSlice(out, N, M)
	duc = extractTimeSlice(du, N, M)

	# compute the cyclic part
	@views Jc(pb, outc, u0[1:end-1-N], h, du[1:end-N-1], outc[:, M])

	# outc[:, M] plays the role of tmp until it is used just after the for-loop
	tmp = @view outc[:, M]

	# compute the partial derivative w.r.t. the period T
	# the .+ is for the GPU
	# out .+= @views (pb(vcat(u0[1:end-1], T .+ δ)) .- pb(u0)) ./ δ .* dT
	@views applyF(pb, tmp, u0c[:, M-1])

	@views POTrapScheme!(pb, outc[:, 1], u0c[:, 1], u0c[:, M-1], dT / (2M), tmp, false)
	for ii = 2:M-1
		@views POTrapScheme!(pb, outc[:, ii], u0c[:, ii], u0c[:, ii-1], dT / (2M), tmp, false)
	end

	# closure condition ensuring a periodic orbit
	outc[:, M] .= @views duc[:, M] .- duc[:, 1]

	# this is for CuArrays.jl to work in the mode allowscalar(false)
	if ongpu(pb)
		return @views vcat(out[1:end-1], dot(duc[:, 1], pb.ϕ))
	else
		out[end] = @views dot(duc[:, 1], pb.ϕ)
		return out
	end
end

function (pb::PeriodicOrbitTrapProblem)(u0::AbstractVector)
	out = similar(u0)
	POTrapFunctional!(pb, out, u0)
end

function (pb::PeriodicOrbitTrapProblem)(u0::AbstractVector, du)
	out = similar(du)
	POTrapFunctionalJac!(pb, out, u0, du)
end

# function (pb::PeriodicOrbitTrapProblem)(out, u0::AbstractVector, du)
# 	# out = similar(du)
# 	POTrapFunctionalJac!(pb, out, u0, du)
# end
####################################################################################################
# Matrix free expressions of matrices related to the Jacobian Matrix of the PO functional
"""
Function to compute the Matrix-Free version of Aγ, see docs for its expression.
"""
function Agamma!(pb::PeriodicOrbitTrapProblem, outc, u0::vectype, du) where {vectype <: AbstractVector}
	# u0 of size N * M + 1
	# du of size N * M
	M = pb.M
	N = pb.N
	T = extractPeriodFDTrap(u0)
	h = T / M

	u0c = extractTimeSlice(u0, N, M)

	# compute the cyclic part
	@views Jc(pb, outc, u0[1:end-1-N], h, du[1:end-N], outc[:, M])

	# closure condition ensuring a periodic orbit
	duc = reshape(du, N, M)
	outc[:, M] .= @views duc[:, M] .- duc[:, 1]
	return nothing
end

"""
Function to compute the Matrix-Free version of the cyclic matrix Jc, see docs for its expression.
"""
function Jc(pb::PeriodicOrbitTrapProblem, outc::AbstractMatrix, u0::vectype, h, du, tmp)  where {vectype <: AbstractVector}
	# tmp plays the role of buffer array
	# u0 of size N * (M - 1)
	# du of size N * (M - 1)
	# outc of size N * M
	# h = T / M
	M = pb.M
	N = pb.N

	u0c = reshape(u0, N, M-1)
	duc = reshape(du, N, M-1)

	@views applyJ(pb, tmp, u0c[:, M-1], duc[:, M-1])

	@views POTrapSchemeJac!(pb, outc[:, 1], u0c[:, 1], u0c[:, M-1],
											duc[:, 1], duc[:, M-1], h/2, tmp)

	for ii = 2:M-1
		@views POTrapSchemeJac!(pb, outc[:, ii], u0c[:, ii], u0c[:, ii-1],
												 duc[:, ii], duc[:, ii-1], h/2, tmp)
	end

	return vec(outc)
end

function Jc(pb::PeriodicOrbitTrapProblem, u0, du)
	M = pb.M
	N = pb.N

	T = extractPeriodFDTrap(u0)
	h = T / M

	out  = similar(du)
	outc = reshape(out, N, M-1)
	tmp = similar(view(outc, :, 1))
	return @views Jc(pb, outc, u0[1:end-1-N], h, du, tmp)
end
####################################################################################################
"""
Matrix by blocks expression of the Jacobian for the PO functional computed at the space-time guess: `u0`
"""
function jacobianPOFD_block(pb::PeriodicOrbitTrapProblem, u0::vectype, γ = 1.0) where {vectype <: AbstractVector}
	# extraction of various constants
	M = pb.M
	N = pb.N
	T = extractPeriodFDTrap(u0)
	h = T / M

	Aγ = BlockArray(spzeros(M * N, M * N), N * ones(Int64, M),  N * ones(Int64, M))
	cylicPOFD_block!(pb, u0, Aγ)

	In = spdiagm( 0 => ones(N))
	setblock!(Aγ, -γ * In, M, 1)
	setblock!(Aγ,  In,     M, M)
	return Aγ
end

"""
This function populates Jc with the cyclic matrix using the different Jacobians
"""
function cylicPOFD_block!(pb::PeriodicOrbitTrapProblem, u0::vectype, Jc::BlockArray) where {vectype <: AbstractVector}
	# extraction of various constants
	M = pb.M
	N = pb.N
	T = extractPeriodFDTrap(u0)
	h = T / M

	In = spdiagm( 0 => ones(N))
	On = spzeros(N, N)

	u0c = extractTimeSlice(u0, N, M)
	outc = similar(u0c)

	tmpJ = @views pb.J(u0c[:, 1])

	@views Jn = In - h/2 .* tmpJ
	setblock!(Jc, Jn, 1, 1)

	@views Jn = -In - h/2 .* pb.J(u0c[:, M-1])
	setblock!(Jc, Jn, 1, M-1)

	for ii=2:M-1
		@views Jn = -In - h/2 .* tmpJ
		setblock!(Jc, Jn, ii, ii-1)

		tmpJ .= @views pb.J(u0c[:, ii])

		@views Jn = In - h/2 .* tmpJ
		setblock!(Jc, Jn, ii, ii)
	end
	return Jc
end

function cylicPOFD_block(pb::PeriodicOrbitTrapProblem, u0::vectype) where {vectype <: AbstractVector}
	# extraction of various constants
	M = pb.M
	N = pb.N
	Jc = BlockArray(spzeros((M - 1) * N, (M - 1) * N), N * ones(Int64, M-1),  N * ones(Int64, M-1))
	cylicPOFD_block!(pb, u0, Jc)
end

cylicPOFD_sparse(pb::PeriodicOrbitTrapProblem, orbitguess0) = blockToSparse(cylicPOFD_block(pb, orbitguess0))

"""
This method returns the jacobian of the functional G encoded in PeriodicOrbitTrapProblem using a Sparse representation.
"""
function (pb::PeriodicOrbitTrapProblem{TF, TJ, TJt, Td2F, vectype, Tls})(::Val{:JacFullSparse}, u0::vectype, γ = 1.0) where {TF, TJ, TJt, Td2F, vectype <: AbstractVector, Tls}
	# extraction of various constants
	M = pb.M
	N = pb.N
	T = extractPeriodFDTrap(u0)
	h = T / M
	Aγ_block = jacobianPOFD_block(pb, u0, γ)

	# we now set up the last line / column
	δ = 1e-9
	@views ∂TGpo = (pb(vcat(u0[1:end-1], T + δ)) .- pb(u0)) ./ δ

	# this is bad for performance. Get converted to SparseMatrix at the next line
	Aγ = blockToSparse(Aγ_block) # most of the computing time is here!!
	@views Aγ = hcat(Aγ, ∂TGpo[1:end-1])
	Aγ = vcat(Aγ, spzeros(1, N * M + 1))

	Aγ[N*M+1, 1:N] .=  pb.ϕ
	Aγ[N*M+1, N*M+1] = ∂TGpo[end]
	return Aγ
end

function (pb::PeriodicOrbitTrapProblem{TF, TJ, TJt, Td2F, vectype, Tls})(::Val{:JacCyclicSparse}, u0::vectype, γ = 1.0) where {TF, TJ, TJt, Td2F, vectype <: AbstractVector, Tls}
	# extraction of various constants
	M = pb.M
	N = pb.N
	T = extractPeriodFDTrap(u0)
	h = T / M
	Aγ_block = jacobianPOFD_block(pb, u0, γ)

	# this is bad for performance. Get converted to SparseMatrix at the next line
	Aγ = blockToSparse(Aγ_block) # most of the computing time is here!!
	return Aγ[1:end-N, 1:end-N]
end

function (pb::PeriodicOrbitTrapProblem{TF, TJ, TJt, Td2F, vectype, Tls})(::Val{:BlockDiagSparse}, u0::vectype) where {TF, TJ, TJt, Td2F, vectype <: AbstractVector, Tls}
	# extraction of various constants
	M = pb.M
	N = pb.N
	T = extractPeriodFDTrap(u0)
	h = T / M
	A_diag_block = BlockArray(spzeros(M * N, M * N), N * ones(Int64, M),  N * ones(Int64, M))

	In = spdiagm( 0 => ones(N))

	u0c = reshape(u0[1:end-1], N, M)
	outc = similar(u0c)

	@views Jn = In - h/2 .* pb.J(u0c[:, 1])
	setblock!(A_diag_block, Jn, 1, 1)

	for ii=2:M-1
		@views Jn = In - h/2 .* pb.J(u0c[:, ii])
		setblock!(A_diag_block, Jn, ii, ii)
	end
	setblock!(A_diag_block, In, M, M)

	A_diag_sp = blockToSparse(A_diag_block) # most of the computing time is here!!
	return A_diag_sp
end
####################################################################################################
# The following struct encodes a jacobian of PeriodicOrbitTrapProblem which is a convenient composite type for the computation of Floquet multipliers
mutable struct PeriodicOrbitTrapJacobianFull{Tpb, Tj, vectype}
	pb::Tpb								# PeriodicOrbitTrapProblem
	J::Tj								# jacobian of the problem
	orbitguess0::vectype				# point at which the jacobian is computed
end

# computation of the jacobian, nothing to be done
(pojacfull::PeriodicOrbitTrapJacobianFull)(x) = return pojacfull

# linear solver for the PO functional, akin to a bordered linear solver
@with_kw mutable struct PeriodicOrbitTrapLS{Tl} <: AbstractLinearSolver
	linsolver::Tl = DefaultLS()	# linear solver
end

# linear solver for the jacobian
(pols::PeriodicOrbitTrapLS)(pojacfull::PeriodicOrbitTrapJacobianFull, rhs) = pols.linsolver(pojacfull.J, rhs)

(pols::PeriodicOrbitTrapLS)(pojacfull::PeriodicOrbitTrapJacobianFull, rhs1, rhs2) = pols.linsolver(pojacfull.J, rhs1, rhs2)

####################################################################################################
# Linear solvers of the linearized version of the functional G implemented by PeriodicOrbitTrapProblem

# composite type to encode the Aγ Operator and its associated cyclic matrix
@with_kw mutable struct AγOperator{Tvec, Tjc, T, Tpb}
	N::Int64 = 0				    		# dimension of a time slice
	orbitguess::Tvec = zeros(1)				# point at which Aγ is evaluated, of size N * M + 1
	Jc::Tjc	= lu(spdiagm(0 => ones(1)))	    # lu factorisation of the cyclic matrix
	is_matrix_free::Bool = false	    	# whether we consider a sparse matrix representation or a Matrix Free one
	γ::T = 1.0				    			# factor γ can be used to compute Floquet multipliers
	prob::Tpb = nothing						# PO functional, used when is_matrix_free = true
end

ismatrixfree(A::AγOperator) = A.is_matrix_free

# linear solver designed specifically to deal with AγOperator
@with_kw struct AγLinearSolver{Tls} <: AbstractLinearSolver
	# Linear solver to invert the cyclic matrix Jc contained in Aγ
	linsolver::Tls = DefaultLS()
end

# linear solver for the PO functional, akin to a bordered linear solver
@with_kw mutable struct PeriodicOrbitTrapBLS{Tl} <: AbstractLinearSolver
	linsolverbls::Tl = BorderingBLS(AγLinearSolver())	# linear solver
end

# The following struct encodes a jacobian of PeriodicOrbitTrapProblem which eases the use of PeriodicOrbitTrapBLS. It is made so that accessing to the cyclic matrix Jc or Aγ is easier. It is combined with a specific linear solver. It is also a convenient structure for the computation of Floquet multipliers
@with_kw mutable struct PeriodicOrbitTrapJacobianBordered{Tpb, T∂, vectype}
	pb::Tpb								# PeriodicOrbitTrapProblem
	∂TGpo::T∂	= nothing				# derivative of the PO functional G wrt T
	Aγ::AγOperator = AγOperator()		# Aγ Operator involved in the Jacobian of the PO functional
	orbitguess0::vectype = nothing		# point at which the jacobian is computed
end

# this function is called whenever the jacobian of G has to be updated
function (J::PeriodicOrbitTrapJacobianBordered)(orbitguess0::AbstractVector)
	# u0 must be an orbit guess
	@views J.orbitguess0 .= orbitguess0[1:length(J.orbitguess0)]

	# we compute the derivative of the problem wrt the period TODO: remove this or improve!!
	δ = 1e-9
	T = extractPeriodFDTrap(orbitguess0)
	# TODO REMOVE CE vcat!
	@views J.∂TGpo .= (J.pb(vcat(orbitguess0[1:end-1], T + δ)) .- J.pb(orbitguess0)) ./ δ

	# update Aγ
	J.Aγ(J.pb, orbitguess0)

	# return J, needed to properly call the linear solver.
	return J
end

# function to update the cyclic matrix
function (Aγ::AγOperator)(pb::PeriodicOrbitTrapProblem, orbitguess::AbstractVector)
	if Aγ.is_matrix_free == false
		# we store the lu decomposition of the newly computed cyclic matrix
		Aγ.Jc = SparseArrays.lu(cylicPOFD_sparse(pb, orbitguess))
	else
		copyto!(Aγ.orbitguess, orbitguess)
		Aγ.prob = pb
	end
end

# this function is called whenever one wants to invert Aγ
function (ls::AγLinearSolver)(A::AγOperator, rhs)
	N = A.N
	if ismatrixfree(A) == false
		# we invert the cyclic part Jc of Aγ
		xbar, flag, numiter = @views ls.linsolver(A.Jc, rhs[1:end - N])
		!flag && @warn "Sparse solver for Aγ did not converge"
	else
		# we invert the cyclic part Jc of Aγ
		xbar, flag, numiter = @views ls.linsolver(dx -> Jc(A.prob, A.orbitguess, dx), rhs[1:end - N])
		!flag && @warn "Matrix Free solver for Aγ did not converge"
	end
	x = similar(rhs)
	x[1:end-N] .= xbar
	x[end-N+1:end] .= @views A.γ .* x[1:N] .+ rhs[end-N+1:end]
	return x ,flag , numiter
end

# Linear solver associated to PeriodicOrbitTrapJacobianBordered
function (ls::PeriodicOrbitTrapBLS)(J::PeriodicOrbitTrapJacobianBordered, rhs)
	N = J.pb.N

	# TODO REMOVE THIS HACK
	ϕ = zeros(length(rhs)-1)
	ϕ[1:N] .= J.pb.ϕ

	# we solve the bordered linear system as follows
	dX, dl, flag, liniter = @views ls.linsolverbls(J.Aγ, J.∂TGpo[1:end-1],
	 										 		  ϕ, J.∂TGpo[end],
													  rhs[1:end-1], rhs[end])

	return vcat(dX, dl), flag, sum(liniter)
end
####################################################################################################
# newton wrapper
"""
	newton(prob::PeriodicOrbitTrapProblem, orbitguess, options::NewtonPar, linearPO = :BorderedLU; kwargs...)

This is the Newton Solver for computing a periodic orbit using a functional G based on Finite Differences and a Trapezoidal rule.

# Arguments:
- `prob` a problem of type `PeriodicOrbitTrapProblem` encoding the functional G
- `orbitguess` a guess for the periodic orbit where `orbitguess[end]` is an estimate of the period of the orbit. It should be a vector of size `N * M + 1` where `M` is the number of time slices, `N` is the dimension of the phase space. This must be compatible with the numbers `N,M` in `prob`.
- `options` same as for the regular `newton` method
- `linearPO = :BorderedLU`. Specify the choice of the linear algorithm, which must belong to `[:FullLU, :BorderedLU, :FullMatrixFree, :BorderedMatrixFree]`. This is used to select a way of inverting the jacobian `dG` of the functional G. For `:FullLU`, we use the default linear solver on a sparse matrix representation of `dG`. For `:BorderedLU`, we take advantage of the bordered shape of the linear solver and use LU decomposition to invert `dG` using a bordered linear solver. This is the default algorithm. For `:FullMatrixFree`, a matrix free linear solver is used for `dG`: note that a preconditioner is very likely required here because of the cyclic shape of `dG` which affects the convergence properties of GMRES. Finally, for `:BorderedMatrixFree`, a matrix free linear solver is used but for `Jc` only (see docs): it means that `options.linsolver` is used to invert `Jc`. These two Matrix-Free options thus expose different part of the jacobian `dG` in order to use specific preconditioners. For example, an ILU preconditioner on `Jc` could remove the constraints in `dG` and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required.

# Output:
- solution
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(probPO::PeriodicOrbitTrapProblem, orbitguess, options::NewtonPar, linearPO::Symbol = :BorderedLU; kwargs...)
	@assert linearPO in [:FullLU, :BorderedLU, :FullMatrixFree, :BorderedMatrixFree]
	N = probPO.N
	M = probPO.M

	if linearPO in [:FullLU, :FullMatrixFree]
		@assert orbitguess isa AbstractVector
		@assert length(orbitguess) == N * M + 1 "Error with size of the orbitguess"

		jac = linearPO == :FullLU ? x -> probPO(Val(:JacFullSparse), x) :
		 								x -> ( dx -> probPO(x, dx))

		return newton(
			x -> probPO(x),
			jac,
			orbitguess, options;
			kwargs...)

	else
		@assert orbitguess isa AbstractVector
		@assert length(orbitguess) == N * M + 1 "Error with size of the orbitguess"

		if linearPO == :BorderedLU
			Aγ = AγOperator(is_matrix_free = false, N = probPO.N)
			Aγ.Jc = lu(spdiagm( 0 => ones(N * (M - 1)) ))
			# linear solver
			lspo = PeriodicOrbitTrapBLS()
		else	# :BorderedMatrixFree
			Aγ = AγOperator(is_matrix_free = true, prob = probPO)
			Aγ.N = probPO.N
			Aγ.orbitguess = zeros(N * M + 1)
			Aγ.Jc = lu(spdiagm( 0 => ones(N * (M - 1)) ))
			# linear solver
			lspo = PeriodicOrbitTrapBLS(BorderingBLS(AγLinearSolver(options.linsolver)))
		end

		# create the jacobian
		JacPO = PeriodicOrbitTrapJacobianBordered(probPO, zeros(N * M + 1), Aγ, zeros(N * M + 1))

		return newton(
			x -> probPO(x),
			x ->  JacPO(x),
			orbitguess, @set options.linsolver = lspo;
			kwargs...)
	end
end

"""
	newton(probPO::PeriodicOrbitTrapProblem, orbitguess, options::NewtonPar, defOp::DeflationOperator{T, Tf, vectype}, linearPO = :BorderedLU; kwargs...) where {T, Tf, vectype}

This function is similar to `newton(probPO, orbitguess, options, linearPO; kwargs...)` except that it uses deflation in order to find periodic orbits different from the one in `defOp`. For example, it can be used in the vicinity of a Hopf bifurcation to prevent the Newton algorithm from converging to the equilibrium point.
"""
function newton(probPO::PeriodicOrbitTrapProblem, orbitguess, options::NewtonPar, defOp::DeflationOperator{T, Tf, vectype}, linearPO::Symbol; kwargs...) where {T, Tf, vectype}
	@assert linearPO in [:FullLU, :FullMatrixFree, :BorderedLU, :BorderedMatrixFree]
	N = probPO.N
	M = probPO.M

	if linearPO in [:FullLU, :FullMatrixFree]
		@assert orbitguess isa AbstractVector
		@assert length(orbitguess) == N * M + 1 "Error with size of the orbitguess"

		jac = linearPO == :FullLU ? x -> probPO(Val(:JacFullSparse), x) :
		 								x -> ( dx -> probPO(x, dx))

		return newton(
			x -> probPO(x),
			jac,
			orbitguess, options, defOp;
			kwargs...)

	else
		@assert orbitguess isa AbstractVector
		@assert length(orbitguess) == N * M + 1 "Error with size of the orbitguess"

		if linearPO == :BorderedLU
			Aγ = AγOperator(is_matrix_free = false, N = probPO.N)
			Aγ.Jc = lu(spdiagm( 0 => ones(N * (M - 1)) ))
			# linear solver
			lspo = PeriodicOrbitTrapBLS()
		else 	#:BorderedMatrixFree
			Aγ = AγOperator(is_matrix_free = true, prob = probPO)
			Aγ.N = probPO.N
			Aγ.orbitguess = zeros(N * M + 1)
			Aγ.Jc = lu(spdiagm( 0 => ones(N * (M - 1)) ))
			# linear solver
			lspo = PeriodicOrbitTrapBLS(BorderingBLS(AγLinearSolver(options.linsolver)))
		end

		# create the jacobian
		JacPO = PeriodicOrbitTrapJacobianBordered(probPO, zeros(N * M + 1), Aγ, zeros(N * M + 1))

		return newton(
			x -> probPO(x),
			x ->  JacPO(x),
			orbitguess, (@set options.linsolver = lspo), defOp;
			kwargs...)
	end

end

####################################################################################################
# continuation wrapper

"""
	continuationPOTrap(prob, orbitguess, p0::Real, contParams::ContinuationPar, linearPO = :BorderedLU; printSolution = (u,p) -> u[end], kwargs...)

This is the continuation routine for computing a periodic orbit using a functional G based on Finite Differences and a Trapezoidal rule.

# Arguments
- `p -> prob(p)` is a family such that `prob(p)::PeriodicOrbitTrapProblem` encodes the functional G
- `orbitguess` a guess for the periodic orbit where `orbitguess[end]` is an estimate of the period of the orbit. It could be a vector of size `N * M + 1` where `M` is the number of time slices, `N` is the dimension of the phase space. This must be compatible with the numbers `N, M` in `prob`.
- `p0` initial parameter, must be a real number
- `contParams` same as for the regular `continuation` method
- `linearPO = :BorderedLU`. Same as `newton` when applied to `PeriodicOrbitTrapProblem`.

Note that by default, the methods prints the period of the periodic orbit as function of the parameter. This can be changed by providing your `printSolution` argument.
"""
function continuationPOTrap(probPO, orbitguess, p0::Real, _contParams::ContinuationPar, linearPO = :BorderedLU; printSolution = (u,p) -> u[end], kwargs...)
	@assert linearPO in [:FullLU, :FullMatrixFree, :BorderedLU, :BorderedMatrixFree]
	contParams = check(_contParams)

	N = probPO(p0).N
	M = probPO(p0).M
	options = contParams.newtonOptions

	if contParams.computeEigenValues
		contParams = @set contParams.newtonOptions.eigsolver = FloquetQaDTrap(contParams.newtonOptions.eigsolver)
	end

	if linearPO in [:FullLU, :FullMatrixFree]
		@assert length(orbitguess) == N * M + 1 "Error with size of the orbitguess"

		jac = linearPO == :FullLU ? (x, p) -> probPO(p)(Val(:JacFullSparse), x) :
		 								(x, p) ->  ( dx -> probPO(p)(x, dx))
		lspo = PeriodicOrbitTrapLS(options.linsolver)
		# contParams.computeEigenValues && @warn "Floquet computation is disabled, choose other linear solver"
		# contParams.detectBifurcation = false
		# contParams.computeEigenValues = false

		return continuation(
			(x, p) -> probPO(p)(x),
			(x, p) -> PeriodicOrbitTrapJacobianFull(probPO(p), jac(x, p), x),
			orbitguess, p0,
			(@set contParams.newtonOptions.linsolver = lspo);
			printSolution = printSolution,
			kwargs...)
	else
		@assert orbitguess isa AbstractVector
		@assert length(orbitguess) == N * M + 1 "Error with size of the orbitguess"

		if linearPO == :BorderedLU
			Aγ = AγOperator(is_matrix_free = false, N = probPO(p0).N)
			Aγ.Jc = lu(spdiagm( 0 => ones(N * (M - 1)) ))
			# linear solver
			lspo = PeriodicOrbitTrapBLS()
		else
			Aγ = AγOperator(is_matrix_free = true, prob = probPO(p0))
			Aγ.N = probPO(p0).N
			Aγ.orbitguess = zeros(N * M + 1)
			Aγ.Jc = lu(spdiagm( 0 => ones(N * (M - 1)) ))
			# linear solver
			lspo = PeriodicOrbitTrapBLS(BorderingBLS(AγLinearSolver(options.linsolver)))
		end

		# create the jacobian
		JacPO = p -> PeriodicOrbitTrapJacobianBordered(probPO(p), zeros(N * M + 1), Aγ, zeros(N * M + 1))

		return continuation(
			(x, p) -> probPO(p)(x),
			(x, p) ->  JacPO(p)(x),
			orbitguess, p0,
			(@set contParams.newtonOptions.linsolver = lspo);
			printSolution = printSolution,
			kwargs...)
	end
end


text:
# Test if we can pass Dual types to Distributions

using Turing
using Distributions
using ForwardDiff: Dual
using Test

float1 = 1.1
float2 = 2.3
d1 = Dual(1.1)
d2 = Dual(2.3)

@test logpdf(Normal(0, 1), d1).value ≈ logpdf(Normal(0, 1), float1) atol=0.001
@test logpdf(Gamma(2, 3), d2).value ≈ logpdf(Gamma(2, 3), float2) atol=0.001
@test logpdf(Beta(2, 3), (d2 - d1) / 2).value ≈ logpdf(Beta(2, 3), (float2 - float1) / 2) atol=0.001

@test pdf(Normal(0, 1), d1).value ≈ pdf(Normal(0, 1), float1) atol=0.001
@test pdf(Gamma(2, 3), d2).value ≈ pdf(Gamma(2, 3), float2) atol=0.001
@test pdf(Beta(2, 3), (d2 - d1) / 2).value ≈ pdf(Beta(2, 3), (float2 - float1) / 2) atol=0.001


text:
using Documenter, FourierFilterFlux
ENV["GKSwstype"] = "nul"
using Pkg
Pkg.build("GR")
makedocs(sitename="FourierFilterFlux.jl",
         pages=[
             "Install" => "installation.md",
             "ConvFFT" => [
                 "Core Type" => "coreType.md",
                 "Initialization" => "init.md",
                 "Boundary Conditions" => "bound.md"
             ],
             "Built-in Constructors" => "constructors.md"
         ])

deploydocs(
    repo="github.com/dsweber2/FourierFilterFlux.jl.git",
)


text:
# GeneticVariation.jl
# ===================
#
# A julia package for the representation and analysis of genetic variation.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE

__precompile__()

module GeneticVariation

export
    # Site types
    Conserved,
    Mutated,
    #Transition,
    #Transversion,
    Segregating,

    # Distances
    Proportion,
    Jaccard,
    MASH,
    distance,
    pdistance,
    mash,
    jaccard,

    # Allele frequencies
    gene_frequencies,

    # Diversity measures
    avg_mut,
    NL79,

    # VCF and BCF
    VCF,
    BCF,
    header,
    metainfotag,
    metainfoval,
    isfilled,
    MissingFieldException

import BioSequences:
    BioSequences,
    Alphabet,
    AA_Term,
    BioSequence,
    bp_chunk_count,
    Certain,
    Composition,
    DNAAlphabet,
    GeneticCode,
    ispurine,
    Kmer,
    Match,
    Mismatch,
    MinHashSketch,
    NucAlphs,
    Position,
    RNAAlphabet,
    Sequence

import BioCore:
    metainfotag,
    metainfoval,
    header

import Combinatorics.permutations
import IntervalTrees: Interval, IntervalValue
import Twiddle:
    enumerate_nibbles,
    nibble_mask,
    count_zero_nibbles,
    count_nonzero_nibbles,
    count_one_nibbles,
    count_zero_bitpairs,
    count_nonzero_bitpairs

include("vcf/vcf.jl")
include("bcf/bcf.jl")
include("site_counting.jl")
include("seg_sites.jl")
include("distances/minhash.jl")
include("distances/proportion.jl")
include("allele_freq.jl")
include("diversity_measures.jl")

end # Module GeneticVariation


text:
#=
    Figure S1
=#

using Model

using CSV
using DataFrames
using DataFramesMeta
using Plots
using StatsPlots

include("../figure_defaults.jl")

##############################################################
## Run/load analysis
##############################################################

include("../../Results/EstimateBinding.jl")

##############################################################
## Fig S1
##############################################################

figS1 = plot(log_γ, LLopt .- maximum(LLopt),
    ylim=(-3.0,0.5),
    frange = -4.0,
    fα = 0.2,
    lw = 2.0,
    c  = :black,
    box = :on,
    widen=false,
    legend=:none,
    xlabel="log(γ)",
    ylabel="Profile Loglikelihood"
)
hline!(figS1,[-1.92],ls=:dash,lw=2.0,c=:red)

##############################################################
## Save...
##############################################################

savefig("$(@__DIR__)/FigS1.svg")

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
#TODO: not used. Plan to remove Plots.jl dependency and add a Plots recipe

@recipe function f(dt::stfitxy,
            framestyle=:axes,
            yminorgrid=true,
            titefont  = font(6, "Computer Modern"),
            ticfont  = font(6, "Computer Modern"),
            guiefont = font(7, "Computer Modern"),
            legndfont = font(6, "Computer Modern"),
            marershape = :circle,
            marersize = 2,
            marerstrokewidth = 0,
            fg_color_legend = nothing,
            size = (600,400),
            dpi = 300,
            margin = 0mm)

    # Default values
    titlefont  := titlefont
    tickfont  := tickfont
    guidefont := guidefont
    legendfont := legendfont
    markershape := markershape
    markersize := markersize
    markerstrokewidth := markerstrokewidth
    fg_color_legend := fg_color_legend
    size := size,
    margin := margin,
    dpi := dpi)

    legend := isnothing(legend) ? (size(x,2) <= 4 ? :outertop : :outerright) : legend

    @series begin
        label := permutedims(names_x)
        color_palette := :tab10
        seriescolor := reshape(collect(1:size(x,2)),1,size(x,2))
        seriestype := :line

        ts,x
    end
end


text:
include(joinpath(@__DIR__, "types.jl"))

function get_yaml_path(platform::Platform)
    lib_dir = @__DIR__
    buildkite_dir = dirname(lib_dir)
    yaml_path = joinpath(buildkite_dir, "test-$(platform.arch).yml")
    return yaml_path
end

max_stars_repo_name: mozilla/rr

max_issues_repo_name: rocallahan/rr

max_forks_repo_name: rocallahan/rr


text:
# Define functions for the Student's t model:

function log_students_t_pdf(x::T,sig::T,nu::T) where {T <: Real}
# -2 log PDF of the Student's t distribution, without normalization:
# fac = gamma(alp)/sqrt(pi*nu)/gamma(0.5*nu)/sig
 # pdf[i] = fac/(1+(x[i]/sig)^2/nu)^alp
return (nu+1)*log(1+(x/sig)^2/nu)
#return pdf/sig
end

# Define -2 log normalization of Student's t distribution:
function log_norm_students_t(x::Array{T,1}) where {T <: Real}
nu,lnV1 = x
# Normalization is: Gamma((\nu+1)/2)/\sqrt{\nu \pi}/Gamma(\nu/2)/\sigma
# (sigma is factor by which the error bars are inflated - it equals sqrt(exp(lnV1)))
return -2*(log(gamma(0.5*(nu+1)))-log(gamma(0.5*nu)))+log(pi*nu)+lnV1
end

# Full log probabilistic model for a single datapoint (missing prior - should add this in [ ]):
function lnprob(p::Array{T,1}) where {T <: Real}
  # Returns -2 log(likelihood) for Student's t distribution.
  # y = observed value; model = model prediction; yerr = uncertainty (sqrt variance)
  # ndof = number of degrees of freedom (\nu) for Student's t distribution;
  # lnV1 = log of factor inflating variance foreground = log(sig^2)
  y, model, yerr, ndof, lnV1 = p

  ll = log_students_t_pdf((y-model)/yerr,exp(0.5*lnV1),ndof)

# return negative log likelihood (so optimize can minimize this,
# or mazimize likelihood).  If the likelihood were exp(-\chi^2/2), which
# it is for ndof = infinity and lnV1 = 0, then this would return \chi^2:
  return ll #+log_norm_students_t([ndof,lnV1])
end

# Compute gradients of this function:
grad_point = p -> ForwardDiff.gradient(lnprob,p)
hess_point = p -> ForwardDiff.hessian(lnprob,p)

# Compute gradients of normalization:
grad_norm = x -> ForwardDiff.gradient(log_norm_students_t,x)
hess_norm = x -> ForwardDiff.hessian(log_norm_students_t,x)


text:
@doc raw"""
`HODE`: Hamiltonian Ordinary Differential Equation

Defines a Hamiltonian ordinary differential initial value problem, that is
a canonical Hamiltonian system of equations,
```math
\begin{aligned}
\dot{q} (t) &= v(t, q(t), p(t)) , & q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f(t, q(t), p(t)) , & p(t_{0}) &= p_{0} ,
\end{aligned}
```
with vector fields ``v`` and ``f``, given by
```math
\begin{aligned}
v &=   \frac{\partial H}{\partial p} , &
f &= - \frac{\partial H}{\partial q} ,
\end{aligned}
```
initial conditions ``(q_{0}, p_{0})`` and the dynamical variables ``(q,p)``
taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}``.

### Parameters

* `DT <: Number`: data type
* `TT <: Real`: time step type
* `AT <: AbstractArray{DT}`: array type
* `vType <: Function`: type of `v`
* `fType <: Function`: type of `f`
* `PType <: Function`: type of `P`
* `hamType <: Function`: Hamiltonian type
* `invType <: OptionalNamedTuple`: invariants type
* `parType <: OptionalNamedTuple`: parameters type
* `perType <: OptionalArray{AT}`: periodicity type

### Fields

* `d`: dimension of dynamical variables ``q`` and ``p`` as well as the vector fields ``v`` and ``f``
* `v`: function computing the vector field ``v``
* `f`: function computing the vector field ``f``
* `poisson`: function computing the Poisson matrix ``P``
* `hamiltonian`: function computing the Hamiltonian ``H``
* `invariants`: either a `NamedTuple` containing the equation's invariants or `nothing`
* `parameters`: either a `NamedTuple` containing the equation's parameters or `nothing`
* `periodicity`: determines the periodicity of the state vector `q` for cutting periodic solutions

### Constructors

```julia
HODE(v, f, poisson, t₀, q₀, p₀, hamiltonian, invariants, parameters, periodicity)

HODE(v, f, h, t₀, q₀::StateVector, p₀::StateVector; kwargs...)
HODE(v, f, h, q₀::StateVector, p₀::StateVector; kwargs...)
HODE(v, f, h, t₀, q₀::State, p₀::State; kwargs...)
HODE(v, f, h, q₀::State, p₀::State; kwargs...)
```

### Keyword arguments:

* `poisson = symplectic_matrix`
* `invariants = nothing`
* `parameters = nothing`
* `periodicity = nothing`

"""
struct HODE{vType <: Callable, fType <: Callable,
            poiType <: Callable,
            hamType <: Callable,
            invType <: OptionalInvariants,
            parType <: OptionalParameters,
            perType <: OptionalPeriodicity} <: AbstractEquationPODE{invType,parType,perType}

    v::vType
    f::fType

    poisson::poiType
    hamiltonian::hamType
    invariants::invType
    parameters::parType
    periodicity::perType

    function HODE(v, f, poisson, hamiltonian, invariants, parameters, periodicity)
        new{typeof(v), typeof(f), typeof(poisson), typeof(hamiltonian),
            typeof(invariants), typeof(parameters), typeof(periodicity)}(
                v, f, poisson, hamiltonian, invariants, parameters, periodicity)
    end
end

HODE(v, f, poisson, hamiltonian; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity()) = HODE(v, f, poisson, hamiltonian, invariants, parameters, periodicity)

GeometricBase.invariants(equation::HODE) = equation.invariants
GeometricBase.parameters(equation::HODE) = equation.parameters
GeometricBase.periodicity(equation::HODE) = equation.periodicity

hasvectorfield(::HODE) = true
hashamiltonian(::HODE) = true

function check_initial_conditions(::HODE, ics::NamedTuple)
    haskey(ics, :q) || return false
    haskey(ics, :p) || return false
    eltype(ics.q) == eltype(ics.p) || return false
    typeof(ics.q) == typeof(ics.p) || return false
    axes(ics.q) == axes(ics.p) || return false
    return true
end

function check_methods(equ::HODE, tspan, ics, params)
    applicable(equ.v, tspan[begin], ics.q, ics.p, zero(ics.q), params) || return false
    applicable(equ.f, tspan[begin], ics.q, ics.p, zero(ics.p), params) || return false
    applicable(equ.hamiltonian, tspan[begin], ics.q, ics.p, params) || return false
    return true
end

function datatype(equ::HODE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    return eltype(ics.q)
end

function arrtype(equ::HODE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    typeof(ics.q)
end

_get_v(equ::HODE, params) = (t,q,p,v) -> equ.v(t, q, p, v, params)
_get_f(equ::HODE, params) = (t,q,p,f) -> equ.f(t, q, p, f, params)
_get_v̄(equ::HODE, params) = _get_v(equ, params)
_get_f̄(equ::HODE, params) = _get_f(equ, params)
_get_h(equ::HODE, params) = (t,q,p) -> equ.hamiltonian(t, q, p, params)
_get_poisson(equ::HODE, params) = (t,q,p,ω) -> equ.poisson(t, q, p, ω, params)

_functions(equ::HODE) = (v = equ.v, f = equ.f, poisson = equ.poisson, h = equ.hamiltonian)
_functions(equ::HODE, params::OptionalParameters) = (v = _get_v(equ, params), f = _get_f(equ, params), poisson = _get_poisson(equ, params), h = _get_h(equ, params))


text:
using Documenter, REDCap
using Literate

makedocs()

deploydocs(
	deps = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4"),
    repo = "github.com/bcbi/REDCap.jl.git",
    julia = "1.0",
    osname = "linux"
)


text:
# By replacing the 1st digit of the 2-digit number *3, it turns out that six of
# the nine possible values: 13, 23, 43, 53, 73, and 83, are all prime.
#
# By replacing the 3rd and 4th digits of 56**3 with the same digit, this
# 5-digit number is the first example having seven primes among the ten
# generated numbers, yielding the family: 56003, 56113, 56333, 56443, 56663,
# 56773, and 56993. Consequently 56003, being the first member of this family,
# is the smallest prime with this property.
#
# Find the smallest prime which, by replacing part of the number (not
# necessarily adjacent digits) with the same digit, is part of an eight prime
# value family.

using ProjectEulerSolutions

# Use Seive to get set of primes, then make a list of possible indexes to
# substitute (we only go up to 6 digits).  Sustitute with new digits and cycle
# through testing for primality until a set of 8 occurs.
function p051solution(n::Integer=100)::Integer

    primes = sieve_eratosthenes(n)
    primeset = Set(primes)
    maxnumprimes = 0
    minprime = n
    checked = Set{Integer}()
    for p in primes
        # Don't need to check primes that have already been tested
        if p in checked
            continue
        end
        numdigs = max(2, ndigits(p))
        mult = (10 .^ (0:numdigs-1))
        # Can't change an even number of digits
        if numdigs == 2
            ss = [[2]]
        elseif numdigs == 3
            ss = [[2], [3]]
        elseif numdigs == 4
            ss = [[2], [3], [4], [2,3,4]]
        elseif numdigs == 5
            ss = [[2], [3], [4], [5], [2,3,4], [2,3,5], [2,4,5]]
        elseif numdigs == 6
            ss = [[2], [3], [4], [5], [6], [2,3,4], [2,3,5], [2,3,6], [2,4,5], [2,4,6], [2,5,6], [2,3,4,5,6]]
        end
        if p > 10
            # Don't use the "last" digit since it would have to include even
            # numbers, which aren't prime
            for s in ss
                d = digits(p)
                numprimes = 0
                minp = 0
                for j = 0:9
                    d[s] .= j
                    newp = sum(d .* mult)
                    if isprime(newp)
                        push!(checked, newp)
                        numprimes += 1
                        if minp == 0
                            minp = newp
                        end
                    end
                end
                if numprimes == 8
                    return minp
                end
            end
        end
    end
    return -1
end

p051 = Problems.Problem(p051solution)

Problems.benchmark(p051, 1_000_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using Gridap
using GridapGmsh

using GridapMHD
using GridapMHD: analytical_shercliff_u
using GridapMHD: analytical_shercliff_j
using GridapMHD: transient_driver_inductionless_MHD
using GridapMHD: driver_inductionless_MHD

# Problem setting
ρ = 1.0
ν = 1.0
σ = 1.0

U0 = 10.0
B0 = 10.0
L = 1.0
Re = U0 * L / ν
Ha = B0 * L * sqrt(σ/(ρ*ν))

# Mesh load
partition = (4,4,3)
domain = (-1,1,-1,1,0,0.1)
model = CartesianDiscreteModel(domain,partition,isperiodic=(false,false,true))

# Boundary conditions location
function hartmann_walls(coord)
  y = coord[2]
  tol = 1e-6
  if (1.0-tol < y ) | ( y < -1.0+tol )
    return true
  end
  false
end

function side_walls(coord)
  x = coord[1]
  tol = 1e-6

  if (1.0-tol < x ) | ( x < -1.0+tol )
    return true
  end
  false
end

function inlet(coord)
  z = coord[3]
  tol = 1e-6
  if (z < tol )
    return true
  end
  false
end

function outlet(coord)
  z = coord[3]
  tol = 1e-6

  if (0.1 - tol < z )
    return true
  end
  false
end

# Name inlet, outlet, and walls element sets.
add_entity!(model,hartmann_walls,"hartmann_walls")
add_entity!(model,side_walls,"side_walls")
add_entity!(model,inlet,"inlet")
add_entity!(model,outlet,"outlet")

K = Ha / (1-0.825*Ha^(-1/2)-Ha^(-1))
∂p∂z = -Re * K / L^3

# Boundary conditions values
g_u(coord,t) = VectorValue(0.0, 0.0, 0.0)
g_u(t::Real) = x -> g_u(x,t)
g_u(x::VectorValue) =  g_u(x,0.0)

g_j(x,t) = VectorValue(0.0,0.0,0.0)
g_j(t::Real) = x -> g_j(x,t)
g_j(x::VectorValue) = g_j(x,0.0)

# Body force
f_u(x) = VectorValue(0.0,0.0, -∂p∂z) * L/U0^2

xh_t, trian, dΩ = transient_driver_inductionless_MHD(;
  t0=0.0,
  tF=1.0,
  Δt=0.25,
  θ=0.5,
  Re=Re,
  Ha=Ha,
  model=model,
  fluid_dirichlet_tags = ["side_walls","hartmann_walls"],
  fluid_neumann_tags = [],
  magnetic_dirichlet_tags = ["side_walls"],
  magnetic_non_perfectly_conducting_walls_tag = ["hartmann_walls",],
  magnetic_neumann_tags = [],
  c_w = 0.0,
  α = 10.0,
  fluid_dirichlet_conditions = g_u ,
  magnetic_dirichlet_conditions = g_j,
  fluid_body_force = f_u,
  constraint_presures = (false,false),
  usegmres = true,
  verbosity = Verbosity(1),
  resultsfile = "output_test"
)


text:
mutable struct TriangleInfo
    total::Int
    perVertex::Vector{Int}
    perEdge::Vector{Int}
end
function TriangleInfo(total::Int, nnodes::Int, nedges::Int)
    return TriangleInfo(total, zeros(Int, nnodes), zeros(Int, nedges))
end

"""
    enumerate_wedges(G::SimpleDiGraph)

Using the degree oriented DAG, enumerate wedges in order count triangles
"""
function enumerate_wedges(DG::SimpleDiGraph)

    n = nv(DG)
    m = ne(DG)

    outinfo = TriangleInfo(0, n, m)

    # # node i
    for i in vertices(DG)
        # # neighbor of i, j
        for (jind, j) in enumerate(neighbors(DG, i))
            # # neighbor of i that isn't j, i.e., the wedge j-i-k.
            # # every neighbor before the jind^th index has already been considered as an endpoint of the wedge. So, start from the j+1st index of neighbors
            for k in neighbors(DG, i)[jind+1:end]

                # # need to check julia docs more carefully, but should be fast since has_edge uses insorted, which assumes what it's looking through is sorted
                if has_edge(DG, j, k) || has_edge(DG, k, j)
                    outinfo.total += 1

                    outinfo.perVertex[i] += 1
                    outinfo.perVertex[j] += 1
                    outinfo.perVertex[k] += 1

                    println("@@@ figure out edge indexing")
                end
                println("wedge $j $i $k")

            end
        end
    end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
Convention used to reconstruct the faces of a polytope
from a list of vertices.
"""
abstract type OrderingConvention end

"""
    connectivities(polytope, rank, convention)

Return the list of `Connectivity` elements used to construct all
faces of rank `rank` of a polytope, following an ordering convention.
"""
function connectivities(polytope::Type{<:Polytope}, rank::Integer,
                        convention::Type{<:OrderingConvention})
  connectivities(polytope, Val(rank), convention)
end

# ---------------------------
# CONVENTION IMPLEMENTATIONS
# ---------------------------

include("conventions/gmsh.jl")


text:
using Dates;
using Printf;
using Distributions;

root_path = raw"D:\AppliRTE\repo\scopf-quanti";
push!(LOAD_PATH, root_path);
cd(root_path);
include(joinpath(root_path, "AmplTxt.jl"));
include(joinpath(root_path, "Workflow.jl"));
include(joinpath(root_path, "DataManager.jl"));

# test_name = "5buses_wind";
test_name = "2buses"
dir_path = joinpath(root_path, test_name);

flexibility = read_flexibility(dir_path);

N_SCENARIO = 5;
df = Dates.DateFormat("Y-m-d-H:M");
h = Dates.DateTime("2015-01-01-11:00", df);
h_1 = Dates.DateTime("2015-01-01-10:00", df);
h_2 = Dates.DateTime("2015-01-01-09:00", df);

# TIME_STEP =  [Minute(0), Minute(15), Minute(30), Minute(45)];
# TIME_STEP =  [Minute(0), Minute(15)];
TIME_STEP =  [Minute(0)];
HORIZON = [h_2, h_1, h];

SCENARIO = [@sprintf("S%d", x) for x in 1:N_SCENARIO];

println("TIME_STEP : ", TIME_STEP);
println("SCENARIO : ", SCENARIO);
println("HORIZON : ", HORIZON);

id_fuel = read_id_fuel(dir_path);
amplTxt = AmplTxt.read(test_name);
num_to_name, bus_load = read_bus_load(amplTxt);
uncertainties, pmax_uncertain = read_uncertainties(dir_path);

println(bus_load);

println(id_fuel);
println(flexibility);

output_file_path = joinpath(dir_path, "pscopf_uncertainties.txt");
write_pscopf_uncertainties(output_file_path, uncertainties, SCENARIO, HORIZON, TIME_STEP);

output_file_path = joinpath(dir_path, "pscopf_previsions.txt");
write_pscopf_previsions(output_file_path, uncertainties, SCENARIO, HORIZON, TIME_STEP);

output_file_path = joinpath(dir_path, "pscopf_units.txt");
write_pscopf_units(output_file_path, amplTxt, id_fuel, flexibility);

output_file_path = joinpath(dir_path, "pscopf_gen_type_bus.txt");
write_pscopf_gen_type_bus(output_file_path, amplTxt, id_fuel, num_to_name);

println(id_fuel);

max_stars_repo_name: azized/pscopf

max_issues_repo_name: azized/pscopf

max_forks_repo_name: azized/pscopf


text:
"""
Stick breaking representation of the simplex.
"""
struct Simplex{P,T,Pm1} <: AbstractVector{T}
    θ::SVector{Pm1,T}
end
Base.length(::Simplex{P}) where P = P
Base.size(::Simplex{P}) where P = (P,)
logit(p) = log(p/(1-p))
inv_logit(y) = 1/(1+exp(-y))

@generated function Simplex(x::NTuple{P,T}) where {P,T}
    q = quote
        y_1 = logit(x[1]) - $(log(1/(P-1)))
        x_sum = zero(T)
    end
    for p ∈ 2:P-1
        push!(q.args, :(x_sum += x[$p]))
        push!(q.args, :( $(Symbol(:y_, p)) = logit(x[$p]/(1-x_sum)) - $(log(1/(P-p))) ) )
    end
    push!(q.args, :(Simplex{$P,$T,$(P-1)}(SVector(@ntuple $(P-1) y))))
    q
end

@generated function Base.getindex(θ::Simplex{P,T,Pm1}, i) where {P,T,Pm1}
    q = quote
        @boundscheck (i > $P || i < 0) && throw(BoundsError())
        x_sum = x = z_1 = inv_logit(θ.θ[1] + $(log(1/(P-1))) )
    end
    for p ∈ 2:min(i,Pm1)
        zsym = Symbol(:z_, p)
        push!(q.args, :( $zsym = inv_logit(θ.θ[$p] + $(log(1/(P-p))) )) )
        push!(q.args, :( x = (1- x_sum)*$zsym ))
        push!(q.args, :( x_sum += x))
    end
    push!(q.args, :( ifelse(i == $P, 1 - x_sum : x) ) )
    q
end
@generated function probabilities(θ::Simplex{P,T,Pm1}) where {P,T,Pm1}
    q = quote
        x_sum = x_1 = z_1 = inv_logit(θ.θ[1] + $(T(log(1/(P-1))) ))
    end
    for p ∈ 2:Pm1
        zsym, xsym = Symbol(:z_, p), Symbol(:x_, p)
        push!(q.args, :( $zsym = inv_logit(θ.θ[$p] + $(T(log(1/(P-p)))) )) )
        push!(q.args, :( $xsym = (1- x_sum)*$zsym ))
        push!(q.args, :( x_sum += $xsym))
    end
    push!(q.args, :( $(Symbol(:x_, P)) = 1 - x_sum ))
    push!(q.args, :( SVector(@ntuple $P p -> x_p) ) )
    q
end
@generated function log_jac(θ::Simplex{P,T,Pm1}) where {P,T,Pm1}
    q = quote
        x_sum = x_1 = z_1 = inv_logit(θ.θ[1] + $(T(log(1/(P-1))) ))
        log_det = log(z_1) + log(1 - z_1)
    end
    for p ∈ 2:Pm1
        zsym, xsym = Symbol(:z_, p), Symbol(:x_, p)
        push!(q.args, :( $zsym = inv_logit(θ.θ[$p] + $(T(log(1/(P-p))) )) ))
        push!(q.args, :( $xsym = (1- x_sum)*$zsym ))
        push!(q.args, :( log_det += log($zsym) + log(1 - $zsym) + log(1 - x_sum) ))
        push!(q.args, :( x_sum += $xsym))
    end
    # push!(q.args, :( $(Symbol(:x_, P)) = 1 - x_sum ))
    push!(q.args, :( SVector(@ntuple $P p -> x_p), log_det ) )
    q
end
@generated function log_jac_probs(θ::Simplex{P,T,Pm1}) where {P,T,Pm1}
    q = quote
        x_sum = x_1 = z_1 = inv_logit(θ.θ[1] + $(T(log(1/(P-1))) ))
        log_det = log(z_1) + log(1 - z_1)
    end
    for p ∈ 2:Pm1
        zsym, xsym = Symbol(:z_, p), Symbol(:x_, p)
        push!(q.args, :( $zsym = inv_logit(θ.θ[$p] + $(T(log(1/(P-p))) )) ))
        push!(q.args, :( $xsym = (1- x_sum)*$zsym ))
        push!(q.args, :( log_det += log($zsym) + log(1 - $zsym) + log(1 - x_sum) ))
        push!(q.args, :( x_sum += $xsym))
    end
    push!(q.args, :( $(Symbol(:x_, P)) = 1 - x_sum ))
    push!(q.args, :( SVector(@ntuple $P p -> x_p), log_det ) )
    q
end


text:
module Krylov

using LinearOperators, LinearAlgebra, SparseArrays, Printf

"Abstract type for statistics returned by a solver"
abstract type KrylovStats end;

"Type for statistics returned by non-Lanczos solvers"
mutable struct SimpleStats <: KrylovStats
  solved :: Bool
  inconsistent :: Bool
  residuals :: Array{Float64,1}
  Aresiduals :: Array{Float64,1}
  status :: String
end

mutable struct LanczosStats <: KrylovStats
  solved :: Bool
  residuals :: Array{Float64}
  flagged :: Union{Bool, Array{Bool,1}, BitArray{1}}
  Anorm :: Float64
  Acond :: Float64
  status :: String
end

mutable struct SymmlqStats <: KrylovStats
  solved :: Bool
  residuals :: Array{Float64}
  residualscg :: Array{Float64}
  errors :: Array{Float64}
  errorscg :: Array{Float64}
  Anorm :: Float64
  Acond :: Float64
  status :: String
end

import Base.show

function show(io :: IO, stats :: SimpleStats)
  s  = "\nSimple stats\n"
  s *= @sprintf("  solved: %s\n", stats.solved)
  s *= @sprintf("  inconsistent: %s\n", stats.inconsistent)
  s *= @sprintf("  residuals:  %s\n", vec2str(stats.residuals))
  s *= @sprintf("  Aresiduals: %s\n", vec2str(stats.Aresiduals))
  s *= @sprintf("  status: %s\n", stats.status)
  print(io, s)
end

function show(io :: IO, stats :: LanczosStats)
  s  = "\nCG Lanczos stats\n"
  s *= @sprintf("  solved: %s\n", stats.solved)
  s *= @sprintf("  residuals: %s\n", typeof(stats.residuals))
  s *= @sprintf("  flagged: %s\n", stats.flagged)
  s *= @sprintf("  ‖A‖F: %7.1e\n", stats.Anorm)
  s *= @sprintf("  κ₂(A): %7.1e\n", stats.Acond)
  s *= @sprintf("  status: %s\n", stats.status)
  print(io, s)
end

function show(io :: IO, stats :: SymmlqStats)
  s  = "\nSYMMLQ stats\n"
  s *= @sprintf("  solved: %s\n", stats.solved)
  s *= @sprintf("  residuals: %s\n", vec2str(stats.residuals))
  s *= @sprintf("  residuals (cg): %s\n", vec2str(stats.residualscg))
  s *= @sprintf("  errors: %s\n", vec2str(stats.errors))
  s *= @sprintf("  errors(cg): %s\n", vec2str(stats.errorscg))
  s *= @sprintf("  ‖A‖F: %7.1e\n", stats.Anorm)
  s *= @sprintf("  κ₂(A): %7.1e\n", stats.Acond)
  s *= @sprintf("  status: %s\n", stats.status)
  print(io, s)
end

include("krylov_aux.jl")
include("krylov_utils.jl")

include("cgs.jl")
include("cg.jl")
include("cg_lanczos.jl")
include("minres.jl")
include("dqgmres.jl")
include("diom.jl")
include("symmlq.jl")
include("cr.jl")

include("cgls.jl")
include("crls.jl")
include("cgne.jl")
include("crmr.jl")

include("lsqr.jl")
include("craig.jl")
include("lsmr.jl")
include("craigmr.jl")
include("lslq.jl")

include("variants.jl")

end


text:
module EfficientNorm.jl

!isdefined(Float) && typealias Float AbstractFloat;

import Base: abs, sqrt, norm

export norm, normsqr, netlib_norm

using ErrorfreeArithmetic
using FloatFloat

include("Constants.jl")
include("NormSqr.jl")
include("Norm.jl")

end # module