text:
"""
  Push an integer or bool.

  The argument is a newline-terminated decimal literal string.
  The intent may have been that this always fit in a short Python int,
  but INT can be generated in pickles written on a 64-bit box that
  require a Python long on a 32-bit box.  The difference between this
  and LONG then is that INT skips a trailing 'L', and produces a short
  int whenever possible.

  Another difference is due to that, when bool was introduced as a
  distinct type in 2.3, builtin names True and False were also added to
  2.2.2, mapping to ints 1 and 0.  For compatibility in both directions,
  True gets pickled as INT + "I01\\n", and False as INT + "I00\\n".
  Leading zeroes are never produced for a genuine integer.  The 2.3
  (and later) unpicklers special-case these and return bool instead;
  earlier unpicklers ignore the leading "0" and return the int.
"""
INT

"""
  Push a four-byte signed integer.

  This handles the full range of Python (short) integers on a 32-bit
  box, directly as binary bytes (1 for the opcode and 4 for the integer).
  If the integer is non-negative and fits in 1 or 2 bytes, pickling via
  BININT1 or BININT2 saves space.
"""
BININT

"""
  Push a one-byte unsigned integer.

  This is a space optimization for pickling very small non-negative ints,
  in range(256).
"""
BININT1

"""
  Push a two-byte unsigned integer.

  This is a space optimization for pickling small positive ints, in
  range(256, 2**16).  Integers in range(256) can also be pickled via
  BININT2, but BININT1 instead saves a byte.
"""
BININT2

"""
  Push a long integer.

  The same as INT, except that the literal ends with 'L', and always
  unpickles to a Python long.  There doesn't seem a real purpose to the
  trailing 'L'.

  Note that LONG takes time quadratic in the number of digits when
  unpickling (this is simply due to the nature of decimal->binary
  conversion).  Proto 2 added linear-time (in C; still quadratic-time
  in Python) LONG1 and LONG4 opcodes.
"""
LONG

"""
  Long integer using one-byte length.

  A more efficient encoding of a Python long; the long1 encoding
  says it all.
"""
LONG1

"""
  Long integer using found-byte length.

  A more efficient encoding of a Python long; the long4 encoding
  says it all.
"""
LONG4

"""
  Push a Python string object.

  The argument is a repr-style string, with bracketing quote characters,
  and perhaps embedded escapes.  The argument extends until the next
  newline character.  These are usually decoded into a str instance
  using the encoding given to the Unpickler constructor. or the default,
  'ASCII'.  If the encoding given was 'bytes' however, they will be
  decoded as bytes object instead.
"""
STRING

"""
  Push a Python string object.

  There are two arguments: the first is a 4-byte little-endian
  signed int giving the number of bytes in the string, and the
  second is that many bytes, which are taken literally as the string
  content.  These are usually decoded into a str instance using the
  encoding given to the Unpickler constructor. or the default,
  'ASCII'.  If the encoding given was 'bytes' however, they will be
  decoded as bytes object instead.
"""
BINSTRING

"""
  Push a Python string object.

  There are two arguments: the first is a 1-byte unsigned int giving
  the number of bytes in the string, and the second is that many
  bytes, which are taken literally as the string content.  These are
  usually decoded into a str instance using the encoding given to
  the Unpickler constructor. or the default, 'ASCII'.  If the
  encoding given was 'bytes' however, they will be decoded as bytes
  object instead.
"""
SHORT_BINSTRING

"""
  Push a Python bytes object.

  There are two arguments:  the first is a 4-byte little-endian unsigned int
  giving the number of bytes, and the second is that many bytes, which are
  taken literally as the bytes content.
"""
BINBYTES

"""
  Push a Python bytes object.

  There are two arguments:  the first is a 1-byte unsigned int giving
  the number of bytes, and the second is that many bytes, which are taken
  literally as the string content.
"""
SHORT_BINBYTES

"""
  Push a Python bytes object.

  There are two arguments:  the first is an 8-byte unsigned int giving
  the number of bytes in the string, and the second is that many bytes,
  which are taken literally as the string content.
"""
BINBYTES8

"""
  Push a Python bytearray object.

  There are two arguments:  the first is an 8-byte unsigned int giving
  the number of bytes in the bytearray, and the second is that many bytes,
  which are taken literally as the bytearray content.
"""
BYTEARRAY8

"Push an out-of-band buffer object."
NEXT_BUFFER

"Make an out-of-band buffer object read-only."
READONLY_BUFFER

"Push None on the stack."
NONE

"Push True onto the stack."
NEWTRUE

"Push False onto the stack."
NEWFALSE

"""
  Push a Python Unicode string object.

  The argument is a raw-unicode-escape encoding of a Unicode string,
  and so may contain embedded escape sequences.  The argument extends
  until the next newline character.
"""
UNICODE

"""
  Push a Python Unicode string object.

  There are two arguments:  the first is a 1-byte little-endian signed int
  giving the number of bytes in the string.  The second is that many
  bytes, and is the UTF-8 encoding of the Unicode string.
"""
SHORT_BINUNICODE

"""
  Push a Python Unicode string object.

  There are two arguments:  the first is a 4-byte little-endian unsigned int
  giving the number of bytes in the string.  The second is that many
  bytes, and is the UTF-8 encoding of the Unicode string.
"""
BINUNICODE

"""
  Push a Python Unicode string object.

  There are two arguments:  the first is an 8-byte little-endian signed int
  giving the number of bytes in the string.  The second is that many
  bytes, and is the UTF-8 encoding of the Unicode string.
"""
BINUNICODE8

"""
  Newline-terminated decimal float literal.

  The argument is repr(a_float), and in general requires 17 significant
  digits for roundtrip conversion to be an identity (this is so for
  IEEE-754 double precision values, which is what Python float maps to
  on most boxes).

  In general, FLOAT cannot be used to transport infinities, NaNs, or
  minus zero across boxes (or even on a single box, if the platform C
  library can't read the strings it produces for such things -- Windows
  is like that), but may do less damage than BINFLOAT on boxes with
  greater precision or dynamic range than IEEE-754 double.
"""
FLOAT

"""
  Float stored in binary form, with 8 bytes of data.

  This generally requires less than half the space of FLOAT encoding.
  In general, BINFLOAT cannot be used to transport infinities, NaNs, or
  minus zero, raises an exception if the exponent exceeds the range of
  an IEEE-754 double, and retains no more than 53 bits of precision (if
  there are more than that, "add a half and chop" rounding is used to
  cut it back to 53 significant bits).
"""
BINFLOAT

"Push an empty list."
EMPTY_LIST

"""
  Append an object to a list.

    Stack before:  ... pylist anyobject
    Stack after:   ... pylist+[anyobject]

  although pylist is really extended in-place.
"""
APPEND

"""
  Extend a list by a slice of stack objects.

    Stack before:  ... pylist markobject stackslice
    Stack after:   ... pylist+stackslice

  although pylist is really extended in-place.
"""
APPENDS

"""
  Build a list out of the topmost stack slice, after markobject.

  All the stack entries following the topmost markobject are placed into
  a single Python list, which single list object replaces all of the
  stack from the topmost markobject onward.  For example,

    Stack before: ... markobject 1 2 3 'abc'
    Stack after:  ... [1, 2, 3, 'abc']
"""
LIST

"Push an empty tuple."
EMPTY_TUPLE

"""
  Build a tuple out of the topmost stack slice, after markobject.

  All the stack entries following the topmost markobject are placed into
  a single Python tuple, which single tuple object replaces all of the
  stack from the topmost markobject onward.  For example,

    Stack before: ... markobject 1 2 3 'abc'
    Stack after:  ... (1, 2, 3, 'abc')
"""
TUPLE

"""
  Build a one-tuple out of the topmost item on the stack.

  This code pops one value off the stack and pushes a tuple of
  length 1 whose one item is that value back onto it.  In other
  words:

      stack[-1] = tuple(stack[-1:])
"""
TUPLE1

"""
  Build a two-tuple out of the top two items on the stack.

  This code pops two values off the stack and pushes a tuple of
  length 2 whose items are those values back onto it.  In other
  words:

      stack[-2:] = [tuple(stack[-2:])]
"""
TUPLE2

"""
  Build a three-tuple out of the top three items on the stack.

  This code pops three values off the stack and pushes a tuple of
  length 3 whose items are those values back onto it.  In other
  words:

      stack[-3:] = [tuple(stack[-3:])]
"""
TUPLE3

"Push an empty dict."
EMPTY_DICT

"""
  Build a dict out of the topmost stack slice, after markobject.

  All the stack entries following the topmost markobject are placed into
  a single Python dict, which single dict object replaces all of the
  stack from the topmost markobject onward.  The stack slice alternates
  key, value, key, value, ....  For example,

    Stack before: ... markobject 1 2 3 'abc'
    Stack after:  ... {1: 2, 3: 'abc'}
"""
DICT

"""
  Add a key+value pair to an existing dict.

    Stack before:  ... pydict key value
    Stack after:   ... pydict

  where pydict has been modified via pydict[key] = value.
"""
SETITEM

"""
  Add an arbitrary number of key+value pairs to an existing dict.

  The slice of the stack following the topmost markobject is taken as
  an alternating sequence of keys and values, added to the dict
  immediately under the topmost markobject.  Everything at and after the
  topmost markobject is popped, leaving the mutated dict at the top
  of the stack.

    Stack before:  ... pydict markobject key_1 value_1 ... key_n value_n
    Stack after:   ... pydict

  where pydict has been modified via pydict[key_i] = value_i for i in
  1, 2, ..., n, and in that order.
"""
SETITEMS

"Push an empty set."
EMPTY_SET

"""
  Add an arbitrary number of items to an existing set.

  The slice of the stack following the topmost markobject is taken as
  a sequence of items, added to the set immediately under the topmost
  markobject.  Everything at and after the topmost markobject is popped,
  leaving the mutated set at the top of the stack.

    Stack before:  ... pyset markobject item_1 ... item_n
    Stack after:   ... pyset

  where pyset has been modified via pyset.add(item_i) = item_i for i in
  1, 2, ..., n, and in that order.
"""
ADDITEMS

"""
  Build a frozenset out of the topmost slice, after markobject.

  All the stack entries following the topmost markobject are placed into
  a single Python frozenset, which single frozenset object replaces all
  of the stack from the topmost markobject onward.  For example,

    Stack before: ... markobject 1 2 3
    Stack after:  ... frozenset({1, 2, 3})
"""
FROZENSET

"Discard the top stack item, shrinking the stack by one item."
POP

"Push the top stack item onto the stack again, duplicating it."
DUP

"""
  Push markobject onto the stack.

  markobject is a unique object, used by other opcodes to identify a
  region of the stack containing a variable number of objects for them
  to work on.  See markobject.doc for more detail.
"""
MARK

"""
  Pop all the stack objects at and above the topmost markobject.

  When an opcode using a variable number of stack objects is done,
  POP_MARK is used to remove those objects, and to remove the markobject
  that delimited their starting position on the stack.
"""
POP_MARK

"""
  Read an object from the memo and push it on the stack.

  The index of the memo object to push is given by the newline-terminated
  decimal string following.  BINGET and LONG_BINGET are space-optimized
  versions.
"""
GET

"""
  Read an object from the memo and push it on the stack.

  The index of the memo object to push is given by the 1-byte unsigned
  integer following.
"""
BINGET

"""
  Read an object from the memo and push it on the stack.

  The index of the memo object to push is given by the 4-byte unsigned
  little-endian integer following.
"""
LONG_BINGET

"""
  Store the stack top into the memo.  The stack is not popped.

  The index of the memo location to write into is given by the newline-
  terminated decimal string following.  BINPUT and LONG_BINPUT are
  space-optimized versions.
"""
PUT

"""
  Store the stack top into the memo.  The stack is not popped.

  The index of the memo location to write into is given by the 1-byte
  unsigned integer following.
"""
BINPUT

"""
  Store the stack top into the memo.  The stack is not popped.

  The index of the memo location to write into is given by the 4-byte
  unsigned little-endian integer following.
"""
LONG_BINPUT

"""
  Store the stack top into the memo.  The stack is not popped.

  The index of the memo location to write is the number of
  elements currently present in the memo.
"""
MEMOIZE

"""
  Extension code.

  This code and the similar EXT2 and EXT4 allow using a registry
  of popular objects that are pickled by name, typically classes.
  It is envisioned that through a global negotiation and
  registration process, third parties can set up a mapping between
  ints and object names.

  In order to guarantee pickle interchangeability, the extension
  code registry ought to be global, although a range of codes may
  be reserved for private use.

  EXT1 has a 1-byte integer argument.  This is used to index into the
  extension registry, and the object at that index is pushed on the stack.
"""
EXT1

"""
  Extension code.

  See [`EXT1`](@ref).  EXT2 has a two-byte integer argument.
"""
EXT2

"""
  Extension code.

  See [`EXT1`](@ref).  EXT4 has a four-byte integer argument.
"""
EXT4

"""
  Push a global object (module.attr) on the stack.

  Two newline-terminated strings follow the GLOBAL opcode.  The first is
  taken as a module name, and the second as a class name.  The class
  object module.class is pushed on the stack.  More accurately, the
  object returned by self.find_class(module, class) is pushed on the
  stack, so unpickling subclasses can override this form of lookup.
"""
GLOBAL

"Push a global object (module.attr) on the stack."
STACK_GLOBAL

"""
  Push an object built from a callable and an argument tuple.

  The opcode is named to remind of the __reduce__() method.

    Stack before: ... callable pytuple
    Stack after:  ... callable(*pytuple)

  The callable and the argument tuple are the first two items returned
  by a __reduce__ method.  Applying the callable to the argtuple is
  supposed to reproduce the original object, or at least get it started.
  If the __reduce__ method returns a 3-tuple, the last component is an
  argument to be passed to the object's __setstate__, and then the REDUCE
  opcode is followed by code to create setstate's argument, and then a
  BUILD opcode to apply  __setstate__ to that argument.

  If not isinstance(callable, type), REDUCE complains unless the
  callable has been registered with the copyreg module's
  safe_constructors dict, or the callable has a magic
  '__safe_for_unpickling__' attribute with a true value.  I'm not sure
  why it does this, but I've sure seen this complaint often enough when
  I didn't want to <wink>.
"""
REDUCE

"""
  Finish building an object, via __setstate__ or dict update.

    Stack before: ... anyobject argument
    Stack after:  ... anyobject

  where anyobject may have been mutated, as follows:

  If the object has a __setstate__ method,

      anyobject.__setstate__(argument)

  is called.

  Else the argument must be a dict, the object must have a __dict__, and
  the object is updated via

      anyobject.__dict__.update(argument)
"""
BUILD

"""
  Build a class instance.

  This is the protocol 0 version of protocol 1's OBJ opcode.
  INST is followed by two newline-terminated strings, giving a
  module and class name, just as for the GLOBAL opcode (and see
  GLOBAL for more details about that).  self.find_class(module, name)
  is used to get a class object.

  In addition, all the objects on the stack following the topmost
  markobject are gathered into a tuple and popped (along with the
  topmost markobject), just as for the TUPLE opcode.

  Now it gets complicated.  If all of these are true:

    + The argtuple is empty (markobject was at the top of the stack
      at the start).
    + The class object does not have a __getinitargs__ attribute.

  then we want to create an old-style class instance without invoking
  its __init__() method (pickle has waffled on this over the years; not
  calling __init__() is current wisdom).  In this case, an instance of
  an old-style dummy class is created, and then we try to rebind its
  __class__ attribute to the desired class object.  If this succeeds,
  the new instance object is pushed on the stack, and we're done.

  Else (the argtuple is not empty, it's not an old-style class object,
  or the class object does have a __getinitargs__ attribute), the code
  first insists that the class object have a __safe_for_unpickling__
  attribute.  Unlike as for the __safe_for_unpickling__ check in REDUCE,
  it doesn't matter whether this attribute has a true or false value, it
  only matters whether it exists (XXX this is a bug).  If
  __safe_for_unpickling__ doesn't exist, UnpicklingError is raised.

  Else (the class object does have a __safe_for_unpickling__ attr),
  the class object obtained from INST's arguments is applied to the
  argtuple obtained from the stack, and the resulting instance object
  is pushed on the stack.

  NOTE:  checks for __safe_for_unpickling__ went away in Python 2.3.
  NOTE:  the distinction between old-style and new-style classes does
         not make sense in Python 3.
"""
INST

"""
  Build a class instance.

  This is the protocol 1 version of protocol 0's INST opcode, and is
  very much like it.  The major difference is that the class object
  is taken off the stack, allowing it to be retrieved from the memo
  repeatedly if several instances of the same class are created.  This
  can be much more efficient (in both time and space) than repeatedly
  embedding the module and class names in INST opcodes.

  Unlike INST, OBJ takes no arguments from the opcode stream.  Instead
  the class object is taken off the stack, immediately above the
  topmost markobject:

    Stack before: ... markobject classobject stackslice
    Stack after:  ... new_instance_object

  As for INST, the remainder of the stack above the markobject is
  gathered into an argument tuple, and then the logic seems identical,
  except that no __safe_for_unpickling__ check is done (XXX this is
  a bug).  See INST for the gory details.

  NOTE:  In Python 2.3, INST and OBJ are identical except for how they
  get the class object.  That was always the intent; the implementations
  had diverged for accidental reasons.
"""
OBJ

"""
  Build an object instance.

  The stack before should be thought of as containing a class
  object followed by an argument tuple (the tuple being the stack
  top).  Call these cls and args.  They are popped off the stack,
  and the value returned by cls.__new__(cls, *args) is pushed back
  onto the stack.
"""
NEWOBJ

"""
  Build an object instance.

  The stack before should be thought of as containing a class
  object followed by an argument tuple and by a keyword argument dict
  (the dict being the stack top).  Call these cls and args.  They are
  popped off the stack, and the value returned by
  cls.__new__(cls, *args, *kwargs) is  pushed back  onto the stack.
"""
NEWOBJ_EX

"""
  Protocol version indicator.

  For protocol 2 and above, a pickle must start with this opcode.
  The argument is the protocol version, an int in range(2, 256).
"""
PROTO

"""
  Stop the unpickling machine.

  Every pickle ends with this opcode.  The object at the top of the stack
  is popped, and that's the result of unpickling.  The stack should be
  empty then.
"""
STOP

"""
  Indicate the beginning of a new frame.

  The unpickler may use this opcode to safely prefetch data from its
  underlying stream.
"""
FRAME

"""
  Push an object identified by a persistent ID.

  The pickle module doesn't define what a persistent ID means.  PERSID's
  argument is a newline-terminated str-style (no embedded escapes, no
  bracketing quote characters) string, which *is* "the persistent ID".
  The unpickler passes this string to self.persistent_load().  Whatever
  object that returns is pushed on the stack.  There is no implementation
  of persistent_load() in Python's unpickler:  it must be supplied by an
  unpickler subclass.
"""
PERSID

"""
  Push an object identified by a persistent ID.

  Like PERSID, except the persistent ID is popped off the stack (instead
  of being a string embedded in the opcode bytestream).  The persistent
  ID is passed to self.persistent_load(), and whatever object that
  returns is pushed on the stack.  See PERSID for more detail.
"""
BINPERSID


text:
"""
This module is a submodule of `GridVisualize`.

It manages a layoutscene with interactive layout and blocking functionality.

Thanks to Julius Krumbiegel for providing  a basic implementation of focus switching.

`GridVisualize` avoids creating dependencies on plotting backends.
So we provide a way to emulate "import Makie" by allowing
to set it as a global variable in the [`setmakie!`](@ref).
As a consequence, we can't use Makie types at compile time.
"""
module FlippableLayout
using DocStringExtensions

Makie=nothing

"""
$(TYPEDEF)

Struct describing flippable layout data.
We don't type annotate with Makie types as they are
unknown at start time.

$(FIELDS)
"""
mutable struct FLayout
    """
    Visible GridLayout
    """
    visible #::GridLayout
    """
    Hidden GridLayot
    """
    offscreen #::GridLayout

    """
    Has the layout been blocked by the block key ?
    """
    blocked::Bool

    """
    Layoutables attached to layout
    """
    layoutables::Dict{Tuple{Int64,Int64},Any} # Union{Makie.MakieLayout.Layoutable,Makie.GridLayout}

    """
    Condition variable working together with the blocked field.
    """
    condition::Condition

    FLayout(visible;blocked=false)=new(visible,
                                       Makie.GridLayout(bbox = Makie.BBox(-500, -400, -500, -400)),
                                       blocked,
                                       Dict{Tuple{Int64,Int64},Any}(),
                                       Condition())
end

function Base.setindex!(flayout::FLayout,layoutable,i,j)
    if isnothing(layoutable)
        flayout.offscreen[1, 1] = flayout.layoutables[(i,j)]
        delete!(flayout.layoutables,(i,j))
    elseif !isa(layoutable,Union{Makie.MakieLayout.Layoutable,Makie.GridLayout})
        error("can only set layoutables")
    else
        flayout.layoutables[(i,j)]=layoutable
        _showall(flayout)
        yield()
    end
end

Base.getindex(flayout::FLayout,i,j)=flayout.layoutables[(i,j)]

function _showall(flayout::FLayout)
    for (pos, layoutable) in flayout.layoutables
        flayout.visible[pos...]=layoutable
    end
    Makie.trim!(flayout.visible)
end

#
# Check if mouse position is  within pixel area of scene
#
function _inarea(area,pos)
    pos[1]>area.origin[1] &&
        pos[1] < area.origin[1]+area.widths[1] &&
        pos[2]>area.origin[2] &&
        pos[2] < area.origin[2]+area.widths[2]
end

function _inscene(l,pos)
    if isa(l,Makie.MakieLayout.Layoutable)
        _inarea(l.scene.px_area[],pos)
    elseif isa(l,Makie.GridLayout)
        _inarea(l.layoutobservables.computedbbox[],pos)
    end
end

"""
    flayoutscene(;blocked=false, kwargs...)

Layoutscene with interactive layout and blocking functionality.

The `,` key switches between focused view showing only one subscene
and "gallery view" showing all layoutables at once.

The space key toggles blocking of the execution of the main therad
when `yield` is replaced by `yieldwait`. Initial blocking state is
set by the `blocked` kwarg.

The `kwargs...` are the same as of `AbstractPlotting.layoutscene`.

The idea is that this can work in some cases as a drop-in replacement
of `layoutscene`.
"""
function flayoutscene(;blocked=false,
                      focuskey=Makie.Keyboard.comma,
                      blockingkey=Makie.Keyboard.space,
                      kwargs...)

    parent = Makie.Scene(; camera = Makie.campixel!, kwargs...)
    layout = Makie.GridLayout(parent,alignmode = Makie.Outside(5))

    Makie.rowgap!(layout,Makie.Relative(0.0))
    Makie.colgap!(layout,Makie.Relative(0.0))

    flayout=FLayout(layout,blocked=blocked)

    gallery_view=Makie.Node(true)

    # Watch mouse position
    mouseposition=Makie.Node((0.0,0.0))

    Makie.on(parent.events.mouseposition) do m
        mouseposition[]=m
        false
    end

    # Switch focus to subscene  at pos
    function _focus(focus)
        for (key,layoutable) in flayout.layoutables
            if key==focus
                flayout.visible[1, 1] = layoutable
            else
                flayout.offscreen[1, 1] = layoutable
            end
        end
        Makie.trim!(flayout.visible)
    end

    # Figure out to which subscene the mouse position
    # corresponds
    function _subscene(mouseposition)
        for (key,layoutable) in flayout.layoutables
            if _inscene(layoutable,mouseposition)
                return key
            end
        end
        return nothing
    end

    function _toggle_block(flayout::FLayout)
        if flayout.blocked
            flayout.blocked=false
            notify(flayout.condition)
        else
            flayout.blocked=true
        end
    end

    # Handle global key events for `,` (focus/gallery view)
    # and space (toggle blocking)
    Makie.on(parent.events.keyboardbutton) do buttons
        if Makie.ispressed(parent, focuskey)
            if gallery_view[]
                pos=_subscene(mouseposition[])
                if !isnothing(pos)
                    _focus(pos)
                end
                gallery_view[]=false
            else
                _showall(flayout)
                gallery_view[]=true
            end
            return true
        end
        if Makie.ispressed(parent, blockingkey)
            _toggle_block(flayout)
            return true
        end
        return false
    end
    parent,flayout
end

"""
     yieldwait(fliplayoutscene)

Yield and wait in case of scene being blocked via space key toggle
"""
function yieldwait(flayout::FLayout)
    yield()
    if flayout.blocked
        wait(flayout.condition)
    end
end

"""
    setmakie!(MyMakie)

Set the Makie module.
This Makie can be GLMakie,WGLMakie,CairoMakie
"""
setmakie!(MyMakie) = global Makie=MyMakie

end

ext: jl


text:
# This file was generated, do not modify it.

using MLJ, RDatasets, PrettyPrinting, Random

Random.seed!(3203)
X = randn(20, 2)
y = vcat(-ones(10), ones(10))

using PyPlot
figure(figsize=(8,6))

ym1 = y .== -1
ym2 = .!ym1
plot(X[ym1, 1], X[ym1, 2], ls="none", marker="o")
plot(X[ym2, 1], X[ym2, 2], ls="none", marker="x")

savefig("assets/literate/ISL-lab-9-g1.svg") # hide

X = MLJ.table(X)
y = categorical(y);

@load SVC pkg=LIBSVM

svc_mdl = SVC()
svc = machine(svc_mdl, X, y)

fit!(svc);

ypred = predict(svc, X)
misclassification_rate(ypred, y)

rc = range(svc_mdl, :cost, lower=0.1, upper=5)
tm = TunedModel(model=svc_mdl, ranges=[rc], tuning=Grid(resolution=10),
                resampling=CV(nfolds=3, rng=33), measure=misclassification_rate)
mtm = machine(tm, X, y)

fit!(mtm)

ypred = predict(mtm, X)
misclassification_rate(ypred, y)


text:
## boundary

∂{A<:Interval,B<:Interval}(d::ProductDomain{Tuple{A,B}})=PiecewiseInterval(d[1].a+im*d[2].a,d[1].b+im*d[2].a,d[1].b+im*d[2].b,d[1].a+im*d[2].b,d[1].a+im*d[2].a)
∂{A<:Interval,B<:PeriodicInterval}(d::ProductDomain{Tuple{A,B}})=UnionDomain((d[1].b+im*d[2],d[1].a+im*reverse(d[2])))
∂{A<:PeriodicInterval,B<:Interval}(d::ProductDomain{Tuple{A,B}})=UnionDomain((d[1]+im*d[2].a,reverse(d[1])+im*d[2].b))
∂{A<:PeriodicInterval,B<:PeriodicInterval}(d::ProductDomain{Tuple{A,B}})=EmptyDomain()

## Union

function Base.join{IT<:Interval}(p1::AbstractVector{IT},p2::AbstractVector{IT})
    for k=length(p1):-1:1,j=length(p2):-1:1
        if p1[k]==reverse(p2[j])
            deleteat!(p1,k)
            deleteat!(p2,j)
            break
        end
    end
    [p1;p2]
end

function ∂{PD1<:ProductDomain,PD2<:ProductDomain}(d::UnionDomain{Tuple{PD1,PD2}})
    bnd=map(d->pieces(∂(d)),d.domains)
    PiecewiseInterval(reduce(join,bnd))
end


text:
#using Plots, LinearAlgebra
#pyplot(framestyle = :box, legend = false, color = :black)

#A node structure has an x,y position in the global space
struct node
    x::Float64 #global x coordinate of node
    y::Float64 #global y coordinate of node
    dof_x::Bool #unconstrained DOF = True; fixed = False
    dof_y::Bool #unconstrained DOF = True; fixed = False
end

#An element structure has two nodes that define the start and end positions of the system
struct element
    a::node #node at start
    b::node #node at end
    A::Float64 #Cross sectional area
    E::Float64 #Young's Modulus
end

#load structure indicates the node that is being loaded + global XY components of the load
struct load
    n::node #node at load point
    x_val::Float64 #X component of force
    y_val::Float64 #Y component of force
end

#Get length of element
function elem_length(elem::element)
    return sqrt((elem.b.x - elem.a.x)^2 + (elem.b.y - elem.a.y)^2)
end

#get mid point of element
function element_mid(e::element)
    x = (e.b.x - e.a.x)/2 + e.a.x
    y = (e.b.y - e.a.y)/2 + e.a.y
    return x, y
end

#Get angle of element
function θ(elem::element)

    #### Angle must be in degrees

    x_global = [1,0] #unit vector in X direciton
    vec_element = [elem.b.x - elem.a.x, elem.b.y - elem.a.y] #element as a vector

    #take dot product of unit x vector and the directional vector of the element
    cos = dot(x_global, vec_element) / (norm(x_global) * norm(vec_element))

    #get the included angle
    angle = acosd(cos)

    #This ensures that the clockwise-θ convention is met
    if vec_element[2] < 0
        angle = 360 - angle
    end

    return angle
end

# Which DOFs at a given node are active/constrained?
function activecheck(n::node)
    #returns tuple for DOF status of global x and y at the given node
    return (n.dof_x, n.dof_y)
end

# Single dimension array of 2n (n = number of nodes) values that indicate the active
# DOFs in order of definition in the array of nodes
function active_dofs(active_map::Array{Tuple{Bool,Bool},1})
    # active_map is an array of tuples indicating the active status of the DOF for each node

    #What is returned is a single 1D array of all active nodal DOFs
    #ie active_dof = [(true, true), (true, false), (false, true)]
    #gives the output [true, true, true, false, false, true]

    return vcat(collect.(active_map)...)
end

function node_element_link(n::node, e::element)
    #Returns a tuple of whether or not the start or end point of a given element is a given node
    return(e.a == n, e.b == n)
end

#Determine which nodes are at the start/end of a given element
function startnode(node_links::Array{Tuple{Bool,Bool},1})
    ans = findfirst([node == (true, false) for node in node_links])
    return ans
end

function endnode(node_links::Array{Tuple{Bool,Bool},1})
    ans = findfirst([node == (false, true) for node in node_links])
    return ans
end

function start_end_node_idx(element, nodes)
    nodelink = [node_element_link(node, element) for node in nodes]
    n_start = startnode(nodelink)
    n_end = endnode(nodelink)
    return [n_start, n_end]
end

# Create the vector of loads given the array of load structures
function loadvector(loads, nodes, n_dof)
    f = zeros(n_dof)
    for l in loads
        idx = 2*findfirst([node == l.n for node in nodes]) - 1
        f[idx:idx+1] = [l.x_val, l.y_val]
    end
    return f
end

## Local coordinate elemental stiffness
function k_localxy(elem::element,length::Float64)
    k = elem.A * elem.E / length .* [1 -1; -1 1]
    return k
end

## Transformation matrix
function Γ(ang)
    return [cosd(ang) sind(ang) 0 0;0 0 cosd(ang) sind(ang)]
end

function k_globalxy(k_local::Array{Float64, 2}, Γ; unit = :degrees)
    #Global_xy stiffness matrix
    k = transpose(Γ) * k_local * Γ

    #Note the indices of the matrix are in the form of:
    # 1: node A, global X
    # 2: node A, global Y
    # 3: node B, global X
    # 4: node B, global Y

    return k
end

## Global coordinate elemental stiffness
## This is the main function that simply outputs the global stiffness matrix
function build_K(K_dim, dof_active, nodes, elements, k_element_global)

    K = zeros(K_dim, K_dim)

    for i = 1:length(elements)

        #for element i, determine which nodes are connected to i, and at what end
        n_start, n_end = start_end_node_idx(elements[i], nodes)

        #Find the index values for the active DOF in the global XY element stiffness matrix
        idx_local_start = [1,2] .* dof_active[n_start]
        idx_local_start = idx_local_start[idx_local_start .!= 0]

        #Find the index values for the active DOf in the global XY element stiffness matrix
        idx_local_end = [3,4] .* dof_active[n_end]
        idx_local_end = idx_local_end[idx_local_end .!= 0]

        #Concatenate the activated index values for the global XY element stiffness matrix
        idx_local = vcat(idx_local_start, idx_local_end)

        #Convert the above calculated values to the index of the GLOBAL stiffness matrix
        idx_global_start =  (n_start * 2 - 2) .+ idx_local_start
        idx_global_end =  (n_end * 2 - 4) .+ idx_local_end
        idx_global = vcat(idx_global_start, idx_global_end)

        #Extracted sub matrix from the global XY element stiffness matrix
        k_extracted = k_element_global[i][idx_local, idx_local]

        K[idx_global, idx_global] .+= k_extracted
    end

    #return stored matrices
    return K
end

## displacement of nodes
function d_glob(disp_long, element, nodes)
    a, b = start_end_node_idx(element, nodes)

    a = a*2 - 1
    b = b*2 - 1
    out_idx = [a, a+1, b, b+1]

    return disp_long[out_idx]
end

## reaction forces
function reactions(F_e, dof_list, nodes, elements)
#create empty storage matrix
F_all = zeros(length(dof_list))

#for each element
for i = 1:length(elements)
    #extract the global XY elemental end forecs
    F = F_e[i]

    #extract the given element information
    e = elements[i]

    #Determine the start and end nodes for the given element
    n_start, n_end = start_end_node_idx(e, nodes)

    #Convert to the indexes in the global DOF list
    idx_start = 2 * n_start - 1 .+ [0,1]
    idx_end = 2 * n_end - 1 .+ [0,1]

    #Add the end forces of each element to the appropriate DOF:
    F_all[idx_start] .+= F[1:2]
    F_all[idx_end] .+= F[3:4]
end

#Return the force values of each DOF that are at REACTIONS
return F_all .* (dof_list .== 0)
end

## Nodes
function displaced_nodes(nodes, disp_long; scale_x = 1, scale_y = 1)
    n_disp = []
    for i = 1:length(nodes)
        n = nodes[i]
        n_new = node((n.x + disp_long[i*2 - 1]) * scale_x, (n.y + disp_long[i*2]) * scale_y, n.dof_x, n.dof_y)
        push!(n_disp, n_new)
    end
    return n_disp
end

## Elements
function displaced_elements(elements, nodes, displaced_nodes)
    elem_disp = []
    for i = 1:length(elements)
        elem = elements[i]
        a, b = start_end_node_idx(elem, nodes)
        elem_new = element(displaced_nodes[a], displaced_nodes[b], elem.A, elem.E)
        push!(elem_disp, elem_new)
    end
    return elem_disp
end

function nodecollector(nodestore::Array)
    L = length(nodestore)
    nodes = [([nodestore[i].x], [nodestore[i].y]) for i = 1:L]
    return nodes
end

function elementcollector(elementstore::Array)
    L = length(elementstore)
    elements = [([elementstore[i].a.x, elementstore[i].b.x], [elementstore[i].a.y, elementstore[i].b.y]) for i = 1:L]
    return elements
end

function loadplotter(loadstore, scale)
    if typeof(loadstore) !== Array{load,1}
        println("No loads specified/incorrect type.")
        return [([0,0], [0,0])]
    end
    #Normalize length of load arrows
    normalizer = maximum([sqrt(load.x_val^2 + load.y_val^2) for load in loadstore])

    L = length(loadstore)

    x1(i) = loadstore[i].n.x
    x2(i) = loadstore[i].n.x + loadstore[i].x_val / normalizer * scale

    y1(i) = loadstore[i].n.y
    y2(i) = loadstore[i].n.y + loadstore[i].y_val / normalizer * scale

    loads = [([x1(i), x2(i)], [y1(i), y2(i)]) for i = 1:L]

    return loads
end

function bc_plotter(nodes)
    #X-direction BCs
    x_fixed_idx = findall([node.dof_x == 0 for node in nodes])
    #Y-direction BCs
    y_fixed_idx = findall([node.dof_y == 0 for node in nodes])

    x_bcs = [([nodes[i].x], [nodes[i].y]) for i in x_fixed_idx]
    y_bcs = [([nodes[i].x], [nodes[i].y]) for i in y_fixed_idx]

    return x_bcs, y_bcs
end

function trussplotter(nodes, elements, loads;
        nodesize = 3,
        nodecolor = :black,
        elementlabels = false,
        elementcolor = :gray,
        elementstyle = :solid,
        lwscale = 2,
        loadscale = :auto,
        bcscale = 10,
        bccolor = :green)

    #initiate plot
    canvas = plot()

    #plot Boundary Conditions
    xbc, ybc = bc_plotter(nodes)
    xmarker = [(-1, sqrt(2)/2), (0,0), (-1,-sqrt(2)/2)]
    ymarker = [(-sqrt(2)/2,-1), (0,0), (sqrt(2)/2,-1)]

    scatter!(xbc, marker = (Shape(xmarker), bcscale, bccolor))
    scatter!(ybc, marker = (Shape(ymarker), bcscale, bccolor))

    #plot elements
    #Linewidth scale
    areas = [element.A for element in elements]
    widths = transpose(areas ./ maximum(areas)) * lwscale

    plot!(elementcollector(elements),
        color = elementcolor,
        linestyle = elementstyle,
        lw = widths,
        label = "")

    #plot element labels
    if elementlabels
        for i = 1:length(elements)
            annotate!((element_mid(elements[i])[1], element_mid(elements[i])[2], Plots.text(string(i), :black)))
        end
    end

    #plot loads
    #Determine scale for loads
    if loadscale == :auto
        xrange = maximum(node.x for node in nodes) - minimum(node.x for node in nodes)
        yrange = maximum(node.y for node in nodes) - minimum(node.y for node in nodes)

        scale = 0.10 * max(xrange, yrange)
    else
        scale = loadscale
    end

    plot!(loadplotter(loads, scale),
        color = :red,
        arrows = true,
        label = "")

    #plot nodes
    scatter!(nodecollector(nodes),
        markercolor = nodecolor,
        label = "")

    canvas
    return canvas
end

function fem_init(nodes, elements)
    #extract lengths of each element
    lengths = [elem_length(e) for e in elements]
    #extract CC angle (degrees) of each element from global axis
    angles = [θ(e) for e in elements]

    #Determine which DOFs are active for each node (X, Y)
    dof_active = [activecheck(node) for node in nodes]

    #Determine array of 2*n_nodes indicating which global DOFs are active
    dof_list = active_dofs(dof_active)

    #Number of total DOFs (active and inactive)
    n_dof = length(dof_list)

    return lengths, angles, dof_active, dof_list, n_dof
end

function stiffness_init(nodes, elements, lengths, angles, dof_active, n_dof)
    #Transformation Matrix
    T = [Γ(angle) for angle in angles]

    #Stiffness Matrices
    k_element_local = [k_localxy(elements[i], lengths[i]) for i = 1:length(elements)]
    k_element_global = [k_globalxy(k_element_local[i], T[i]) for i = 1:length(elements)]

    #Global stiffness matrix
    K = build_K(n_dof, dof_active, nodes, elements, k_element_global)

    return T, k_element_local, k_element_global, K
end

function displacements(K, F, dof_list, n_dof)
    #Reduce K matrix and force vector to non-zero items:
    K_reduced = K[dof_list, dof_list]
    F_reduced = F[dof_list]

    #compute displacement vector
    disp = K_reduced \ F_reduced #K^-1 F

    compliance = transpose(F_reduced) * disp #F^T d

    #Extended disp vector
    disp_long = zeros(n_dof)
    disp_long[dof_list] = disp

    #Convert global displacements to elemental global XY displacements
    disp_global = [d_glob(disp_long, elements[i], nodes) for i = 1:length(elements)]

    return disp_long, disp_global, compliance
end

function forces(nodes, elements, k_global, disp_global, T, dof_list)
    #Global element forces
    F_e = [k_global[i] * disp_global[i] for i = 1:length(elements)]

    #Convert global element forces t olocal element forces
    f_e = [T[i] * F_e[i] for i = 1:length(elements)]

    #Single axial force value
    f = [force[2] for force in f_e]

    #Stress values
    σ = f ./ [element.A for element in elements]

    #Determing Reaction Forces
    rxns = reactions(F_e, dof_list, nodes, elements)

    return F_e, f_e, f, σ, rxns
end

function eq_check(rxns, F, tolerance)
    #Check equilibrium
    if sum(rxns) + sum(F) > tolerance
        println("Equilibrium not met!!")
        return true
    else
        println("Equilibrium met.")
    end
end

function new_positions(nodes, elements, displacements)
    new_nodes = displaced_nodes(nodes, displacements)
    new_elements = displaced_elements(elements, nodes, new_nodes)

    return new_nodes, new_elements
end

function fem2d_solver(nodes, elements, loads; tol = 1e-6)

    #Initial geometric parameters + activity check of Degrees of Freeodm
    lengths, angles, dof_active, dof_list, n_dof = fem_init(nodes, elements)

    #Return global stiffness matrix (including 0 rows/columns)
    T, k_local, k_global, K = stiffness_init(nodes, elements, lengths, angles, dof_active, n_dof)

    #Create force vector
    F = loadvector(loads, nodes, n_dof)

    #Displacements in global coordinate system + elemental DOF displacements
    disp, disp_elemental, compliance = displacements(K, F, dof_list, n_dof)

    F_global, f_elemental, f_axial, stress_axial, rxns = forces(nodes, elements, k_global, disp_elemental, T, dof_list)

    #If equilibrium is not reached, exit solver
    if eq_check(rxns, F, tol) == true
        return
    end

    #Deformed nodes and elements
    new_nodes, new_elements = new_positions(nodes, elements, disp)

    return new_nodes, new_elements, disp, f_axial, stress_axial, rxns, compliance
end

function fem2d_solver_modified(nodes, elements, loads; tol = 1e-6)

    #Initial geometric parameters + activity check of Degrees of Freeodm
    lengths, angles, dof_active, dof_list, n_dof = fem_init(nodes, elements)

    #Return global stiffness matrix (including 0 rows/columns)
    T, k_local, k_global, K = stiffness_init(nodes, elements, lengths, angles, dof_active, n_dof)

    #Create force vector
    F = loadvector(loads, nodes, n_dof)

    #Displacements in global coordinate system + elemental DOF displacements
    disp, disp_elemental, compliance = displacements(K, F, dof_list, n_dof)

    return compliance
end

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


text:
module HeterogeneousAgentsContinuousTime

greet() = print("Hello World!")

using SparseArrays, LinearAlgebra, Plots

function HJBequation(r0, lambda2)
    #This function solves the HJB equation using the implicit method in Achdou et al (2017)
    #r0 is the interest rate and lambda2 the probability to pass to state 2

    #Maximum iterations
    maxiter=100
    #Tolerance
    tol=1e-6
    #Delta to calculate vprime
    delta=100

    #Price
    #r=0.03
    r=r0

    #Discount
    ρ=0.05

    #Parameter for the CRRA
    σ=2.0
    #Inverse of the marginal utility
    uprimeinv(x)=x^(-1/σ)
    #utility
    u(x)=(x^(1-σ))/(1-σ)
    #marginal utility
    uprime(x)=x^(-σ)

    z=[0.1,0.2]#States of the stochastic process
    Nz=length(z)

    #λ=[1.2,1.2]#Transition probabilities
    λ=[0.6,lambda2]#Transition probabilities

    amin=-0.15
    amax=4.0
    Na=500 #number of points in a
    a=[range(amin,stop=amax,length=Na);]
    deltaa=(amax-amin)/(Na-1) #Assumes that a is equally spaced

    #Create big vector with z
    zVec=[]
    for j=1:Nz
        zVec=[zVec;z[j]*ones(Na)]
    end

    v0=zeros(eltype(a),Na*Nz)
    for j=1:Nz
        v0[(j-1)*Na+1:j*Na]=((z[j].+max(r,0.01).*a).^(1-σ))/(1-σ)/ρ
    end

    #definition of the positions that are different from zero in the matrix A
    rows=[1:Na;Na+1:2*Na;1:2*Na;2:Na;Na+2:2*Na;1:Na-1;Na+1:2*Na-1]#Positions of \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
    cols=[Na+1:2*Na;1:Na;1:2*Na;1:Na-1;Na+1:2*Na-1;2:Na;Na+2:2*Na]#Positions of \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
    #Initialization of A
    valsA=ones(eltype(v0),length(rows))
    A=sparse(rows,cols,valsA)

    v1=zeros(eltype(a),Na*Nz)
    c=zeros(eltype(a),Na*Nz)
    s=zeros(eltype(a),Na*Nz)

    for i=1:maxiter
        vprimeF=ones(eltype(a),Na*Nz)
        vprimeB=ones(eltype(a),Na*Nz)
        for j=1:Nz
            vprimeF[(j-1)*Na+1:j*Na-1]=(v0[(j-1)*Na+2:j*Na].-v0[(j-1)*Na+1:j*Na-1])./deltaa #In the last position the forward difference is not defined
            vprimeF[j*Na]=(z[j] + r.*amax).^(-σ) #Will never be used, but just in case
            vprimeB[(j-1)*Na+2:j*Na]=(v0[(j-1)*Na+2:j*Na].-v0[(j-1)*Na+1:j*Na-1])./deltaa #In the first position the backward difference is not defined
            vprimeB[(j-1)*Na+1]=(z[j] + r.*amin).^(-σ) #state constraint boundary condition
        end
        sF=zeros(eltype(vprimeF),Na*Nz)
        sB=zeros(eltype(vprimeB),Na*Nz)

        ######CAREFUL. SOMETIMES VPRIME IS NEGATIVE, THIS ENSURES THIS DOESN'T HAPPEN
        vprimeF[vprimeF.<=0].=tol
        vprimeB[vprimeB.<=0].=tol

        for j=1:Nz
            for k=1:Na
                sF[(j-1)*Na+k]=z[j]+r*a[k]-uprimeinv(vprimeF[(j-1)*Na+k])
                sB[(j-1)*Na+k]=z[j]+r*a[k]-uprimeinv(vprimeB[(j-1)*Na+k])
            end
        end
        vprime=similar(vprimeF)
        for j=1:Nz
            for k=1:Na
                if sF[(j-1)*Na+k]>0
                    vprime[(j-1)*Na+k]=vprimeF[(j-1)*Na+k]
                elseif sB[(j-1)*Na+k]<0
                    vprime[(j-1)*Na+k]=vprimeB[(j-1)*Na+k]
                else
                    vprime[(j-1)*Na+k]=uprime(z[j]+r*a[k])
                end
            end
        end
        #Consumption
        c=uprimeinv.(vprime)

        #Calculation of savings
        s=zVec+r*repeat(a,Nz,1)-c
        b=u.(c)+v0/delta

        x=-sB/deltaa
        x[x.<=0].=0
        w=sF/deltaa #This is the z that goes in the matrix A. It has a conflict with the stochastic process
        w[w.<=0].=0
        #Making the first element of x (for each value of z) 0
        for j=1:Nz
            x[(j-1)*Na+1]=0
        end
        for j=1:Nz
            w[j*Na]=0
        end
        #Definition of y without substracting the probability lambda
        y=-x-w
        #Adjusting y by the probability lambda
        for j=1:Nz
            y[(j-1)*Na+1:j*Na]=y[(j-1)*Na+1:j*Na].-λ[j]
        end

        #Creation of matrix a
        #THIS MATRIX ASSUMES LENGTH(Z)=2
        valsA=vec([λ[1]*ones(eltype(λ[1]),Na,1);λ[2]*ones(eltype(λ[1]),Na,1);y;x[2:Na];x[Na+2:2*Na];w[1:Na-1];w[Na+1:2*Na-1]]) #vector with the values of A. \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
        A=sparse(rows,cols,valsA)

        B=(1/delta+ρ)*I-A
        v1=B\b

        if maximum(abs.(v1-v0))<tol
            break
        end
        v0=copy(v1)
    end

    #p1=plot(a,[s[1:Na] s[Na+1:2*Na]],xlabel="Wealth",title="Savings",label=["y1" "y2"])
    #p2=plot(a,[c[1:Na] c[Na+1:2*Na]],xlabel="Wealth",title="Consumption",label=["y1" "y2"],legend=:topleft)
    #display(plot(p1,p2,layout=(1,2)))

    return (a=a,z=z,c=c,s=s,A=A,v=v0)
end

function KF(a,z,A)
    Na=length(a)
    deltaa=(a[Na]-a[1])/(Na-1)
    Nz=length(z)
    g=zeros(eltype(a),Nz*Na)
    zeroVec=zeros(eltype(a),Nz*Na) #Zero vector in equation 20 of the Numerical appendix
    zeroVec[1]=0.1
    (rowsA,colsA,valsA)=findnz(A) #retrieving rows, columns and values
    m,n=size(A)
    rowsAt=[colsA;ones(n)] #the columns of A and a vector of ones to change first row
    colsAt=[rowsA;collect(1:n)]
    valsAt=[valsA;1;zeros(eltype(a),n-1)]
    At=sparse(rowsAt,colsAt,valsAt)
    gtilda=At\zeroVec
    g=gtilda/(deltaa*sum(gtilda))

    #p1=plot(a[1:Na],g[1:Na],xlabel="Wealth",title="Distribution of agents",label=["y1" "y2"])
    #p2=plot!(a,g[Na+1:2*Na])
    #display(plot(p1,p2))

    return g
end

function AssetSupply()

    #Values to iterate over r
    rmin=-0.05
    rmax=0.04
    Nr=20
    vecr=[range(rmin,stop=rmax,length=Nr);]

    #Supply function
    S=zeros(eltype(vecr),Nr)

    lambda2=0.5

    for i in eachindex(vecr)
        results=HJBequation(vecr[i],lambda2)
        g=KF(results.a,results.z,results.A)
        Na=length(results.a)
        deltaa=(results.a[Na]-results.a[1])/(Na-1)
        S[i]=deltaa*sum(g.*repeat(results.a,2,1))
    end
    plot(S,vecr,xlabel="Net demand", ylabel="r", title="Asset supply", legend=false)
    plot!([0], seriestype="vline")
    #@show vecr
end

function equilibriumR(lambda2)
    #Calculates the equilibrium interest rate

    #Maximum r to consider
    rmax=0.04
    #Minimum r to consider
    rmin=0.01
    #weight of previous r
    wr=0.5

    r0=0.03 #Initial guess

    #Maximum iterations
    itermax=100
    tol=1e-6

    #Excess demand
    S=0

    #Transition probability
    #lambda2=1.2

    for i=1:itermax
        results=HJBequation(r0,lambda2)
        g=KF(results.a,results.z,results.A)
        Na=length(results.a)
        deltaa=(results.a[Na]-results.a[1])/(Na-1)
        S=deltaa*sum(g.*repeat(results.a,2,1))#Total supply
        #@show r0 rmin rmax S
        #@show results.c[1:10]
        if S>tol/100 #&& (rmax-rmin)>tol
            rmax=copy(r0)
            r0=wr*r0+(1-wr)*rmin
        elseif S<-tol/100 #&& (rmax-rmin)>tol #Excess demand
            rmin=copy(r0)
            r0=wr*r0+(1-wr)*rmax
        else
            break
        end
    end
    return r0
    #println(r0, S)
end

function TimeHJBequation(rl, resultsT,T,lambda2)
    #Given a time path for r (rl) between t=0 and t=T, it computes the value and policy functions for each t
    #resultsT are a tuple given by HJBequation which contains the results for the steady state (t=T)
    Nt=length(rl)
    deltat=T/Nt

    #Maximum iterations
    maxiter=100
    #Tolerance
    tol=1e-6

    #Discount
    ρ=0.05

    #Parameter for the CRRA
    σ=2.0
    #Inverse of the marginal utility
    uprimeinv(x)=x^(-1/σ)
    #utility
    u(x)=(x^(1-σ))/(1-σ)
    #marginal utility
    uprime(x)=x^(-σ)

    z=[0.1,0.2]#States of the stochastic process
    Nz=length(z)

    λ=[0.6,lambda2]#Transition probabilities

    amin=-0.15
    amax=4.0
    Na=500 #number of points in a
    a=[range(amin,stop=amax,length=Na);]
    deltaa=(amax-amin)/(Na-1) #Assumes that a is equally spaced

    #Create big vector with z
    zVec=[]
    for j=1:Nz
        zVec=[zVec;z[j]*ones(Na)]
    end

    v=zeros(eltype(a),Na*Nz,Nt+1)

    #Last v is the stationary state value function

    v[:,Nt+1]=resultsT.v

    #definition of the positions that are different from zero in the matrix A
    rows=[1:2*Na;1:2*Na;2:Na;Na+2:2*Na;1:Na-1;Na+1:2*Na-1]#Positions of \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
    cols=[Na+1:2*Na;1:Na;1:2*Na;1:Na-1;Na+1:2*Na-1;2:Na;Na+2:2*Na]#Positions of \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
    #Initialization of A
    valsA=ones(eltype(v),length(rows))
    #In this version A is a vector of sparse matrices
    #A=zeros(2*Na,2*Na,Nt+1)

    A=[]
    push!(A,resultsT.A)
    #A[:,:,Nt+1]=resultsT.A

    #Initialization of Atemp
    #Atemp=sparse(rows,cols,valsA)

    c=zeros(eltype(a),Na*Nz,Nt)
    s=zeros(eltype(a),Na*Nz,Nt)

    for i=Nt:-1:1
        vprimeF=ones(eltype(a),Na*Nz)
        vprimeB=ones(eltype(a),Na*Nz)
        for j=1:Nz
            vprimeF[(j-1)*Na+1:j*Na-1]=(v[(j-1)*Na+2:j*Na,i+1].-v[(j-1)*Na+1:j*Na-1,i+1])./deltaa #In the last position the forward difference is not defined
            vprimeF[j*Na]=(z[j] + rl[i].*amax).^(-σ) #Will never be used, but just in case
            vprimeB[(j-1)*Na+2:j*Na]=(v[(j-1)*Na+2:j*Na,i+1].-v[(j-1)*Na+1:j*Na-1,i+1])./deltaa #In the first position the backward difference is not defined
            vprimeB[(j-1)*Na+1]=(z[j] + rl[i].*amin).^(-σ) #state constraint boundary condition
        end
        sF=zeros(eltype(vprimeF),Na*Nz)
        sB=zeros(eltype(vprimeB),Na*Nz)

        ######CAREFUL. SOMETIMES VPRIME IS NEGATIVE, THIS ENSURES THIS DOESN'T HAPPEN
        vprimeF[vprimeF.<=0].=tol
        vprimeB[vprimeB.<=0].=tol

        for j=1:Nz
            for k=1:Na
                sF[(j-1)*Na+k]=z[j]+rl[i]*a[k]-uprimeinv(vprimeF[(j-1)*Na+k])
                sB[(j-1)*Na+k]=z[j]+rl[i]*a[k]-uprimeinv(vprimeB[(j-1)*Na+k])
            end
        end
        vprime=similar(vprimeF)

        for j=1:Nz
            for k=1:Na
                if sF[(j-1)*Na+k]>0
                    vprime[(j-1)*Na+k]=vprimeF[(j-1)*Na+k]
                elseif sB[(j-1)*Na+k]<0
                    vprime[(j-1)*Na+k]=vprimeB[(j-1)*Na+k]
                else
                    vprime[(j-1)*Na+k]=uprime(z[j]+rl[i]*a[k])
                end
            end
        end
        #Consumption
        c[:,i]=uprimeinv.(vprime)

        #Calculation of savings
        s[:,i]=zVec+rl[i]*repeat(a,Nz,1)-c[:,i]
        b=u.(c[:,i])+v[:,i+1]/deltat

        x=-sB/deltaa
        x[x.<=0].=0
        w=sF/deltaa #This is the z that goes in the matrix A. It has a conflict with the stochastic process
        w[w.<=0].=0
        #Making the first element of x (for each value of z) 0
        for j=1:Nz
            x[(j-1)*Na+1]=0
        end
        for j=1:Nz
            w[j*Na]=0
        end
        #Definition of y without substracting the probability lambda
        y=-x-w
        #Adjusting y by the probability lambda
        for j=1:Nz
            y[(j-1)*Na+1:j*Na]=y[(j-1)*Na+1:j*Na].-λ[j]
        end

        #Creation of matrix a
        #THIS MATRIX ASSUMES LENGTH(Z)=2
        #valsA=vcat(λ[1]*ones(eltype(λ[1]),Na,1),λ[2]*ones(eltype(λ[1]),Na,1),y,x[2:Na],x[Na+2:2*Na],w[1:Na-1],w[Na+1:2*Na-1]) #vector with the values of A. \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
        valsA=vec([λ[1]*ones(eltype(λ[1]),Na,1);λ[2]*ones(eltype(λ[1]),Na,1);y;x[2:Na];x[Na+2:2*Na];w[1:Na-1];w[Na+1:2*Na-1]]) #vector with the values of A. \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
        Atemp=sparse(rows,cols,valsA)

        B=(1/deltat+ρ)*I-Atemp
        v[:,i]=B\b
        pushfirst!(A,Atemp)
    end
    #@show typeof(A[:,:,3])
    return (a=a,z=z,c=c,s=s,A=A[1:end-1],v=v[:,1:Nt])
end

function TimeKF(a,z,A,g0,rl,T)
    #Calculates distribution over time period implied by rl. Taking g0 as initial distribution and A as the matrices at each period in time
    Na=length(a)
    deltaa=(a[Na]-a[1])/(Na-1)
    Nz=length(z)
    Nt=length(rl)
    deltat=T/Nt

    g=zeros(eltype(a),Nz*Na,Nt+1)
    g[:,1]=g0

    for i=2:Nt+1
        Atemp1=A[i-1]
        (rowsA,colsA,valsA)=findnz(Atemp1) #retrieving rows, columns and values
        Atemp=sparse(colsA,rowsA,valsA)#Defining the tranpose
        g[:,i]=(I-deltat*sparse((Atemp)))\g[:,i-1]
    end
    return g
end

function shockSim()
    #Maximum iterations
    maxiter=100
    #Tolerance
    tol=0.001

    #Initial value of lambda2
    lambda2_0=0.5 #1.0
    #Final value of lambda2
    lambda2_T=0.9 #1.5

    r0=equilibriumR(lambda2_0)
    results0=HJBequation(r0, lambda2_0)
    a=results0.a
    Na=length(a)
    deltaa=(a[end]-a[1])/(Na-1)

    rT=equilibriumR(lambda2_T)
    resultsT=HJBequation(rT, lambda2_T)

    g0=KF(results0.a,results0.z,results0.A)
    gT=KF(resultsT.a,resultsT.z,resultsT.A)

    T=20
    #Number of points in the time interval
    N=100

    #Matrix to save g
    g=repeat(g0,1,N+1)

    #speed of updating the interest rate. Taken directly from  Moll's code
    xi = 20*(exp.(-0.05*(1:N)) .- exp.(-0.05*N))

    #First guess on the path of r
    rl0=rT*ones(eltype(rT),N)
    rl1=copy(rl0)
    @show r0 rT

    for i=1:maxiter
        results=TimeHJBequation(rl0, resultsT,T, lambda2_T)
        g=TimeKF(results0.a,results0.z,results.A,g0,rl0,T)
        S=deltaa.*sum(g.*repeat(a,2,N+1),dims=1)
        diffS=S[2:end]-S[1:end-1]
        #@show maximum(abs.(diffS))
        rl1=rl0-xi.*diffS
        #@show maximum(abs.(rl1-rl0))
        if maximum(abs.(diffS))<tol || maximum(abs.(rl1-rl0))<tol
            break
        end
        #println(rl1)
        rl0=copy(rl1)
    end
    #@show rl0
    p1=plot([range(-5*T/N,stop=T,length=N+5);],[r0*ones(5); rl1],xlabel="T",ylabel="r",legend=false)
    p2=plot(a[1:100], [g[1:100,1] g[Na+1:Na+100,1] g[1:100,3] g[Na+1:Na+100,3]],xlabel="Wealth",title="t=0.06",label=["y1 t=0" "y2  t=0" "y1 t=0.06" "y2  t=0.06"],linecolor=["red" "blue" "red" "blue"],linestyle=[:solid :solid :dash :dash])
    p3=plot(a[1:100], [g[1:100,1] g[Na+1:Na+100,1] g[1:100,N+1] g[Na+1:Na+100,N+1]],xlabel="Wealth",title="Long run",label=["y1 t=0" "y2  t=0" "y1 t=T" "y2  t=T"],linecolor=["red" "blue" "red" "blue"],linestyle=[:solid :solid :dash :dash])
    display(plot(p1,p2,p3,layout=(1,3)))
end

function optimalh(a, ϕ, r, p, h,f)
    #Finds the optimal h given a value of wealth a. In this case the input a is a value, while h is a vector of possible values of h.
    #From the paper it should have this shape={0,[hmin,\infty]}
    #f is the function of housing in the utility function

    #Possible values of utility for all h
    temp=f.(h)
    temp[p*h.>ϕ*a].=-Inf
    #Position of the optimal h
    hopt=argmax(temp)
    #Value at optimum
    ftilda=f(h[hopt])

    return (ftilda,hopt)

end

function HJBequationHousing(r0, lambda2)
        #This function solves the HJB equation using the implicit method in Achdou et al (2017)
        #r0 is the interest rate and lambda2 the probability to pass to state 2

        #Maximum iterations
        maxiter=120
        #Tolerance
        tol=1e-6
        #Delta to calculate vprime
        delta=1000

        #Price
        #r=0.03
        r=r0

        #Discount
        ρ=0.05

        #Parameter for the CRRA
        σ=2.0
        #Inverse of the marginal utility
        uprimeinv(x)=x^(-1/σ)
        #utility
        u(x)=(x^(1-σ))/(1-σ)
        #marginal utility
        uprime(x)=x^(-σ)

        z=[0.1,0.135]#States of the stochastic process
        Nz=length(z)

        #λ=[1.2,1.2]#Transition probabilities
        λ=[0.5,lambda2]#Transition probabilities

        amin=0.0
        amax=3.0
        Na=500 #number of points in a
        a=[range(amin,stop=amax,length=Na);]
        deltaa=(amax-amin)/(Na-1) #Assumes that a is equally spaced

        #Create big vector with z
        zVec=[]
        for j=1:Nz
            zVec=[zVec;z[j]*ones(Na)]
        end

        #Housing parameters
        hmin=2.3
        p=1
        ϕ=2
        #Parameters of f
        α=1/3
        η=0.2
        #Utility of housing
        f(x)=η*(max(x-hmin,0))^α - r*p*x

        #Vector with possible values of housing
        Nh=10000
        h=[0;range(hmin,stop=amax*ϕ/p,length=(Nh-1));]

        #Finding optimal values for h given each value of a
        hopt=zeros(Float64,Na)
        fopt=zeros(eltype(h),Na)
        hopt=min.((r*p/(α*η))^(1/(α-1))+hmin,ϕ*a./p)
        hopt[hopt.<hmin].=0
        fopt=f.(hopt)
        #for i=1:Na
        #    (fopt[i],hopt[i])=optimalh(a[i], ϕ, r, p, h,f)
        #end
        #@show hopt
        v0=zeros(eltype(a),Na*Nz)
        for j=1:Nz
            v0[(j-1)*Na+1:j*Na]=((z[j].+max(r,0.01).*a).^(1-σ))/(1-σ)/ρ
        end

        #definition of the positions that are different from zero in the matrix A
        rows=[1:Na;Na+1:2*Na;1:2*Na;2:Na;Na+2:2*Na;1:Na-1;Na+1:2*Na-1]#Positions of \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
        cols=[Na+1:2*Na;1:Na;1:2*Na;1:Na-1;Na+1:2*Na-1;2:Na;Na+2:2*Na]#Positions of \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
        #Initialization of A
        valsA=ones(eltype(v0),length(rows))
        A=sparse(rows,cols,valsA)

        v1=zeros(eltype(a),Na*Nz)
        cx=zeros(eltype(a),Na*Nz)#This is defined as c+f(h)
        c=zeros(eltype(a),Na*Nz)
        s=zeros(eltype(a),Na*Nz)

        for i=1:maxiter
            vprimeF=ones(eltype(a),Na*Nz)
            vprimeB=ones(eltype(a),Na*Nz)
            for j=1:Nz
                vprimeF[(j-1)*Na+1:j*Na-1]=(v0[(j-1)*Na+2:j*Na].-v0[(j-1)*Na+1:j*Na-1])./deltaa #In the last position the forward difference is not defined
                vprimeF[j*Na]=(z[j] + fopt[Na] +r.*amax).^(-σ) #Will never be used, but just in case
                vprimeB[(j-1)*Na+2:j*Na]=(v0[(j-1)*Na+2:j*Na].-v0[(j-1)*Na+1:j*Na-1])./deltaa #In the first position the backward difference is not defined
                vprimeB[(j-1)*Na+1]=(z[j] + fopt[1] +r.*amin).^(-σ) #state constraint boundary condition
            end
            sF=zeros(eltype(vprimeF),Na*Nz)
            sB=zeros(eltype(vprimeB),Na*Nz)

            ######CAREFUL. SOMETIMES VPRIME IS NEGATIVE, THIS ENSURES THIS DOESN'T HAPPEN
            vprimeF[vprimeF.<=0].=10^(-10)
            vprimeB[vprimeB.<=0].=10^(-10)

            for j=1:Nz
                for k=1:Na
                    #@show j k
                    sF[(j-1)*Na+k]=z[j]+fopt[k]+r*a[k]-uprimeinv(vprimeF[(j-1)*Na+k])
                    sB[(j-1)*Na+k]=z[j]+fopt[k]+r*a[k]-uprimeinv(vprimeB[(j-1)*Na+k])
                end
            end
            vprime=similar(vprimeF)
            #vprime for the first value
            #for j=1:Nz
            #    if sF[(j-1)*Na+1]>0
            #        vprime[(j-1)*Na+1]=vprimeF[(j-1)*Na+1]
            #    else
            #        vprime[(j-1)*Na+1]=uprime(z[j]+r*amin)
            #    end
            #end
            #vprime for the last value
            #for j=1:Nz
            #    if sB[j*Na]<0
            #        vprime[j*Na]=vprimeB[j*Na]
            #    else
            #        vprime[j*Na]=uprime(z[j]+r*amax)
            #    end
            #end
            #vprime for intermediate values
            indexB=zeros(Int64,Na*Nz) #1 if using backward
            indexF=zeros(Int64,Na*Nz) #1 if using forward
            for j=1:Nz
                for k=1:Na
                    #utility plus continuation value forward
                    uF=u(uprimeinv(vprimeF[(j-1)*Na+k]))+vprimeF[(j-1)*Na+k]*sF[(j-1)*Na+k]
                    uB=u(uprimeinv(vprimeB[(j-1)*Na+k]))+vprimeB[(j-1)*Na+k]*sB[(j-1)*Na+k]
                    if ((sF[(j-1)*Na+k]>0) & (sB[(j-1)*Na+k]>=0)) | ((sF[(j-1)*Na+k]>0) & (sB[(j-1)*Na+k]<0) & (uF>=uB))
                        vprime[(j-1)*Na+k]=vprimeF[(j-1)*Na+k]
                        indexF[(j-1)*Na+k]=1
                    elseif ((sB[(j-1)*Na+k]<0) & (sF[(j-1)*Na+k]<=0)) | ((sB[(j-1)*Na+k]<0) & (sF[(j-1)*Na+k]>0) & (uB>=uF))
                        vprime[(j-1)*Na+k]=vprimeB[(j-1)*Na+k]
                        indexB[(j-1)*Na+k]=1
                    else
                        vprime[(j-1)*Na+k]=uprime(z[j]+fopt[k]+r*a[k])
                    end
                end
            end
            #Consumption
            cx=uprimeinv.(vprime)

            #Calculation of savings
            s=zVec+repeat(fopt,Nz,1)+r*repeat(a,Nz,1)-cx
            b=u.(cx)+v0/delta

            x=-sB/deltaa
            x[indexB.==0].=0
            #@show x
            w=sF/deltaa #This is the z that goes in the matrix A. It has a conflict with the stochastic process
            w[indexF.==0].=0
            #Making the first element of x (for each value of z) 0
            for j=1:Nz
                x[(j-1)*Na+1]=0
            end
            for j=1:Nz
                w[j*Na]=0
            end
            #Definition of y without substracting the probability lambda
            y=-x-w
            #Adjusting y by the probability lambda
            for j=1:Nz
                y[(j-1)*Na+1:j*Na]=y[(j-1)*Na+1:j*Na].-λ[j]
            end

            #Creation of matrix a
            #THIS MATRIX ASSUMES LENGTH(Z)=2
            #valsA=vcat(λ[1]*ones(eltype(λ[1]),Na,1),λ[2]*ones(eltype(λ[1]),Na,1),y,x[2:Na],x[Na+2:2*Na],w[1:Na-1],w[Na+1:2*Na-1]) #vector with the values of A. \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
            valsA=vec([λ[1]*ones(eltype(λ[1]),Na,1);λ[2]*ones(eltype(λ[1]),Na,1);y;x[2:Na];x[Na+2:2*Na];w[1:Na-1];w[Na+1:2*Na-1]]) #vector with the values of A. \lambda1 and \lambda 2. Positions of y. Positions for x1. Positions for x2. Positions for z1. Positions for z2
            A=sparse(rows,cols,valsA)

            B=(1/delta+ρ)*I-A
            v1=B\b

            if maximum(abs.(v1-v0))<tol
                break
            end
            v0=copy(v1)
        end

        c=cx.-repeat(η*(max.(hopt.-hmin,0)).^α,Nz)

        #p1=plot(a,[s[1:Na] s[Na+1:2*Na]],xlabel="Wealth",title="Savings",label=["y1" "y2"],legend=:topright)
        #p2=plot(a,[v1[1:Na] v1[Na+1:2*Na]],legend=false)
        #p3=plot(a, hopt,xlabel="Wealth",title="Housing",legend=false)
        #p4=plot(a,[c[1:Na] c[Na+1:2*Na]],xlabel="Wealth",title="Consumption",label=["y1" "y2"],legend=:bottomright)
        #p5=plot(a,fopt,legend=false)
        #display(plot(p1,p3,p4,layout=(1,3)))
        #@show p*hmin/ϕ
        #@show h[hopt]
        return (a=a,z=z,c=c,s=s,h=hopt,A=A,v=v0)
end

function KFHousing(a,z,A)
    Na=length(a)
    deltaa=(a[Na]-a[1])/(Na-1)
    Nz=length(z)
    g=ones(eltype(a),Nz*Na)
    g=g/(deltaa*sum(g))
    gtilda=copy(g)
    N=1000#Number of times until convergence to distribution
    dt=10
    At=A'
    for i=1:N
        gtilda=(I - At*dt)\g
        g=copy(gtilda)
    end
    #g=gtilda/(deltaa*sum(gtilda))
    p1=plot(a[1:Na],g[1:Na],xlabel="Wealth",title="Distribution of agents",label=["y1" "y2"])
    p2=plot!(a,g[Na+1:2*Na])
    display(plot(p1,p2))

    return g
end

end # module


text:
module Homebrew

import Base: show
using Unicode, JSON, Libdl, InteractiveUtils

# Find homebrew installation prefix
const brew_prefix = abspath(joinpath(dirname(@__FILE__),"..","deps", "usr"))
const brew_exe = joinpath(brew_prefix,"bin","brew")

# Types and show() overrides, etc..
include("types.jl")

# Utilities
include("util.jl")

# The public API that uses Homebrew like a blackbox, only through the `brew` script
include("API.jl")

# The private API that peeks into the internals of Homebrew a bit
include("private_API.jl")
include("translation.jl")

# Include our own, personal bindeps integration stuff
include("bindeps_integration.jl")

"""
`__init__()`

Initialization function.  Calls `install_brew()` to ensure that everything we
need is downloaded/installed, then calls `update_env()` to set the environment
properly so that packages being installed can find their binaries.
"""
function __init__()
    if Sys.isapple()
        # Let's see if Homebrew is installed.  If not, let's do that first!
        (isdir(brew_prefix) && isdir(tappath)) || install_brew()

        # Update environment variables such as PATH, DL_LOAD_PATH, etc...
        update_env()
    else
        # change this to an error in future
        @warn("""
Homebrew.jl can only be used on Apple macOS. Suggested usage is

    @static if Sys.isapple()
        using Homebrew
        # Homebrew specific code goes here
    end
""")
    end
end

end # module


text:
function qrsolve!(A::AbstractMatrix{T},b::AbstractVector{T}) where {T}
    # Fast least-squares solver for full column rank Hessenberg-like matrices
    m, n = size(A)
    m < n && error("Column dimension exceeds row dimension")
    _, τ = LinearAlgebra.LAPACK.geqrf!(A)
    T <: Complex ? tran = 'C' : tran = 'T'
    LinearAlgebra.LAPACK.ormqr!('L',tran,A,τ,view(b,:,1:1))
    return UpperTriangular(triu(A[1:n,:]))\b[1:n]
end

function poldeg1(v::Vector{T}) where {T <: Number}
    # Degree of a polynomial plus one.
    n1 = findlast(!iszero,v) # degree(v) + 1
    n1 === nothing && (n1 = 0)
    return n1
end
poldeg(v::Vector{T}) where {T <: Number} = poldeg1(v)-1  # Degree of a polynomial
function poldiv(a::Vector{T}, b::Vector{S}) where {T <: Number, S <: Number }
    # Quotient of the exact division of two polynomials a and b.
    # q = poldiv(a,b) returns the quotient q = a/b, which is computed by solving the
    # linear least-square problem ||C(b)*q - a|| = min, where C(b) is the convolution
    # matrix of appropriate size.

    na1 = poldeg1(a) # degree(a) + 1
    nb1 = poldeg1(b) # degree(b) + 1
    nb1 == 0 && throw(DivideError())
    na1 < nb1 && error("Degree of second polynomial exceeds the degree of first polynomial")
    a1, b1  = promote(float(a[1:na1]), float(b[1:nb1]))
    #return Matrix(Toeplitz([den;zeros(n1-length(den))],[den[1];zeros(n1-m1)]))\num[1:n1]
    #return convmtx(b[1:nb1],na1-nb1+1)\a[1:na1]
    return qrsolve!(convmtx(b1,na1-nb1+1),a1)
end
function convmtx(v::Vector{T}, n::Int) where {T <: Number}

    #   Convolution matrix.
    #   C = convmtx(v,n) returns the convolution matrix C for a vector v.
    #   If q is a column vector of length n, then C*q is the same as conv(v,q).

    #   Form C as the Toeplitz matrix
    #   C = Toeplitz([v; zeros(n-1)],[v[1]; zeros(n-1));  put Toeplitz code inline

    nv = length(v)-1
    C = zeros(T, n+nv, n)
    @inbounds for j = 1:n
        C[j:j+nv,j] = v
    end
    return C
end
function poldivrem(a::Vector{T}, b::Vector{S}) where {T <: Number, S <: Number }
    # Quotient and remainder of the division of two polynomials a and b.
    # q, r = poldivrem(a,b) returns the quotient q and remainder r such that a = q*b + r.
    na1 = poldeg1(a) # degree(a) + 1
    nb1 = poldeg1(b) # degree(b) + 1
    nb1 == 0 && throw(DivideError())

    R = eltype(one(T)/one(S))

    if na1 < nb1
       return zeros(R, 1), a
    end
    m = nb1-1
    q = zeros(R, na1-m)
    r = R[ a[i] for i in 1:na1 ]
    #r = copy_oftype(num[1:n1],R)

    @inbounds for i in na1:-1:nb1
        s = r[i] / b[nb1]
        q[i-m] = s
        @inbounds for j in 1:m
            r[i-nb1+j] -= b[j] * s
        end
    end

    return q, m == 0 ? zeros(R,1) : r[1:m]
end
function polgcdvw(a::Vector{T}, b::Vector{S}; atol::Real = 0, rtol::Real = Base.rtoldefault(float(real(T))), maxnit::Int = 0) where {T <: Number, S <: Number }
    # Greatest common divisor of two polynomials.
    #     polgcdvw(a, b; atol = 0, rtol, maxnit = 0) -> (d, v, w, δ)
    # Compute the greatest common divisor `d` of two polynomials `a` and `b`, and
    # the polynomials v and w such that a = d * v * ka and b = d * w *kb, where ka and kb
    # are scalar scaling factors (ka and kb are not explicitly provided).
    # The accuracy estimate δ is computed as δ = norm([a/ka-d*v; b/kb-d*w]).
    # A SVD-based method adapted from [1] is employed and an iterative accuracy refinement
    # of δ, with a maximum number of maxnit iterations, is additionally performed.
    #  [1] Z. Zeng, The numerical greatest common divisor of univariate polynomials,
    #      in: L. Gurvits, P. Pébay, J.M. Rojas, D. Thompson (Eds.),
    #      Randomization, Relaxation,and Complexity in Polynomial Equation Solving,
    #      Contemporary Mathematics,vol.556,AMS,2011,pp.187–217.

    na1 = poldeg1(a) # degree(a) + 1
    nb1 = poldeg1(b) # degree(b) + 1

    R = eltype(one(T)/one(S))

    a1, b1  = promote(float(a[1:na1]), float(b[1:nb1]))

    na1 <= 0 && (return ones(R,1), zeros(R,1), b1, zero(R))
    nb1 <= 0 && (return ones(R,1), a1, zeros(R,1), zero(R))
    a1 = a1/norm(a1)
    b1 = b1/norm(b1)
    switch = (na1 < nb1)
    switch && ((a1, b1, na1, nb1) = (b1, a1, nb1, na1))

    atol == 0 ? tol = rtol : tol = atol
    # determine the degree of GCD as the nullity of the Sylvester matrix
    nd = na1 + nb1 - 2 - rank([convmtx(a1,nb1-1) convmtx(b1,na1-1)], atol = tol)
    nd == 0 && (switch ? (return [one(R)], b1, a1, zero(R)) : (return [one(R)], a1, b1, zero(R)))
    # determine [w; -v] from an orthogonal/unitary nullspace basis of dimension 1
    # of a reduced Sylvester matrix using the last row of Vt from its (partial) SVD
    _, sv, Vt = LAPACK.gesvd!('N', 'A', [convmtx(a1,nb1-nd) convmtx(b1,na1-nd)])
    k = na1 + nb1 -2*nd - count(x -> x > tol, sv) # expected nullity k = 1
    k == 1 || error("GCD computation failure")
    wv = Vt[end,:]
    # wv = LAPACK.gesvd!('N', 'A', [convmtx(a1,nb1-nd) convmtx(b1,na1-nd)])[3][end,:]
    eltype(wv) <: Complex && (wv = conj(wv))
    v = -wv[nb1-nd+1:end]
    w = wv[1:nb1-nd]
    # determine the GCD d with high accuracy as the solution of a well-conditioned
    # linear least-squares problem
    # d = convmtx(w,nd+1)\b1
    # d = [convmtx(v,nd+1) ; convmtx(w,nd+1)] \ [a1; b1]
    d = qrsolve!([convmtx(v,nd+1) ; convmtx(w,nd+1)],[a1; b1])
    if maxnit > 0
       d, v, w, δ = gcdvwupd(a1, b1, d, v, w, maxnit = maxnit)
    else
       δ = norm( [a1; b1] - [ conv(v,d); conv(w,d) ])
    end

    switch ? (return d, w, v, δ) : (return d, v, w, δ)
end
function gcdvwupd(a::Vector{T}, b::Vector{T}, d::Vector{T}, v::Vector{T}, w::Vector{T}; maxnit::Int = 10) where {T <: Number}
    # Iterative refinement of the accuracy of the greatest common divisor of two polynomials.
    #     gcdvwupd(a, b, d, v, w; maxnit = 0) -> (dupd, vupd, wupd, δ)
    # Given the polynomials a and b, and a triple of polynomials (d, v, w) such that
    # `d` is an approximation of a greatest common divisor of `a` and `b`, and
    # v and w are approximate quotients of a/d and b/d, respectively, compute the updated triple
    # (dupd, vupd, wupd) which minimizes the accuracy estimate δ = norm([a-dupd*vupd; b-dupd*wupd]).
    # A maximum number of maxnit iterations are performed using the Gauss-Newton iteration method,
    # in the form proposed in [1].
    #
    #  [1] Z. Zeng, The numerical greatest common divisor of univariate polynomials,
    #      in: L. Gurvits, P. Pébay, J.M. Rojas, D. Thompson (Eds.),
    #      Randomization, Relaxation,and Complexity in Polynomial Equation Solving,
    #      Contemporary Mathematics,vol.556,AMS,2011,pp.187–217.

    na1 = length(a)
    nb1 = length(b)
    nd1 = length(d)
    nv1 = length(v)
    nw1 = length(w)
    na1 == nd1+nv1-1 || error("Incompatible dimensions between a, d and v")
    nb1 == nd1+nw1-1 || error("Incompatible dimensions between b, d and w")

    h = d/norm(d)^2
    f = [a;b;fill(one(T),1)]
    fh = [conv(v,d); conv(w,d); fill(dot(h,d),1)]
    Δf = f-fh;
    δ = norm(Δf)
    d0 = copy(d); v0 = copy(v); w0 = copy(w);
    nvw1 = nv1+nw1
    ndvw1 = nd1+nvw1
    nab1 = na1+nb1
    Jh = [zeros(nab1,ndvw1); zeros(1,nvw1) h' ]

    j1 = 1:nv1
    j2 = nv1+1:nvw1
    j3 = nvw1+1:ndvw1
    i1 = 1:nv1+nd1-1
    i2 = nv1+nd1:nvw1+2*nd1-2
    Jh11 = view(Jh,i1,j1)
    Jh22 = view(Jh,i2,j2)
    Jh13 = view(Jh,i1,j3)
    Jh23 = view(Jh,i2,j3)

    for it = 1:maxnit
        # Jh = [ convmtx(d0,nv1)          0     convmtx(v0,nd1) ]
        #      [  0             convmtx(u0,nw1) convmtx(w0,nd1) ]
        #      [  0                       0            h'       ]
        Jh11[:,:] = convmtx(d0,nv1)
        Jh22[:,:] = convmtx(d0,nw1)
        Jh13[:,:] = convmtx(v0,nd1)
        Jh23[:,:] = convmtx(w0,nd1)
        # Δz = Jh\Δf;   # this does not exploit the zero structure
        Δz = qrsolve!(copy(Jh),copy(Δf));  # this is probably the fastest
        # Δz = qr(Jh)\Δf;  # this is usually slower, but faster for full matrices
        v1 = v0 + Δz[j1];
        w1 = w0 + Δz[j2];
        d1 = d0 + Δz[j3];
        fh1 = [conv(v1,d1); conv(w1,d1); fill(dot(h,d1),1)];
        Δf1 = f-fh1;
        δ1 = norm(Δf1)
        if δ1 < δ
           d0 = copy(d1); v0 = copy(v1); w0 = copy(w1); δ = copy(δ1); Δf = copy(Δf1);
        else
           break
        end
    end
    return d0, v0, w0, δ
end
function pollcm(a::Vector{T}, b::Vector{S}; atol::Real = 0, rtol::Real = Base.rtoldefault(float(real(T))), maxnit::Int = 10) where {T <: Number, S <: Number }
    # Least common multiple of two polynomials.
    #     m = pollcm(a,b; atol = 0, rtol, maxnit = 0)
    # Compute the least common multiple `m` of two polynomials `a` and `b`. Employ iterative refinement of
    # the accuracy of GCD(a,b) if maxnit > 0.

    na1 = poldeg1(a) # degree(a) + 1
    nb1 = poldeg1(b) # degree(b) + 1

    na1 == 0 && (return zeros(T,1))
    nb1 == 0 && (return zeros(S,1))
    na1 == 1 && (return b)
    nb1 == 1 && (return a)

    return conv(a, polgcdvw(a, b, atol = atol, rtol = rtol, maxnit = maxnit)[3])
end
function conv(a::Vector{T}, b::Vector{S}) where {T <: Number, S <: Number }
    # Convolution of two vectors (or product of two polynomials).
    # c = conv(a,b) returns the convolution of vectors a and b.

    na1 = length(a)
    nb1 = length(b)

    R = promote_type(T, S)
    c = zeros(R, na1 + nb1 - 1)
    for i in 1:na1, j in 1:nb1
        @inbounds c[i + j - 1] += a[i] * b[j]
    end
    return c
end
function polcoeffval(r::AbstractVector{T},val::Number) where {T}
    # Compute the coefficients of a polynomial from its roots and evaluate the polynomial
    # for a given value of its argument. Both are equal to one for an empty vector r.
    T1 = promote_type(T,eltype(val))
    n = length(r)
    c = zeros(T1, n+1)
    ONE = one(T1)
    c[1] = ONE
    pval = ONE
    for j = 1:n
        pval = pval*(val-r[j])
        for i = j:-1:1
            c[i+1] = c[i+1]-r[j]*c[i]
        end
    end
    return reverse(c), pval
 end


text:
struct Font
    surface::Cairo.CairoSurface
    charDict::Dict{Char, Tuple{Int, Int, Int, Int}}
    charSpacing::Int
    lineSpacing::Int

    Font(filename::String, fontString::String, charSpacing::Int=1, lineSpacing::Int=1) = new(open(Cairo.read_from_png, filename), generateFontDict(fontString), charSpacing, lineSpacing)
end

const stringLineAlignOffsets = Dict{Symbol, Function}(
    :left => (widest, width) -> 0,
    :center => (widest, width) -> round(Int, (widest - width) / 2),
    :right => (widest, width) -> widest - width,
)

function getQuadTuple(row::Int, col::Int, width::Int=3, height::Int=5)
    return (8 * row - 8, 8 * col - 8, width, height)
end

# Use first occurance of a symbol as the target quad
function generateFontDict(fontString::String)
    fontMap = Dict{Char, Tuple{Int, Int, Int, Int}}()

    for (i, row) in enumerate(split(fontString, "\n"))
        for (j, c) in enumerate(row)
            get!(fontMap, c) do
                getQuadTuple(j, i, 3, 5)
            end
        end
    end

    return fontMap
end

# Multiline string
function lineWidth(font::Font, s::AbstractString, replaceMissing::Char=' ')
    return isempty(s) ? 0 : sum([get(font.charDict, c, font.charDict[replaceMissing])[3] for c in s]) + (length(s) - 1) * font.charSpacing
end

function lineHeight(font::Font, s::AbstractString, replaceMissing::Char=' ')
    return isempty(s) ? 0 : maximum([get(font.charDict, c, font.charDict[replaceMissing])[4] for c in s])
end

function Base.size(font::Font, s::AbstractString, replaceMissing::Char=' ')
    lines = split(s, "\n")

    width = maximum(lineWidth.(Ref(font), lines))
    height = sum(lineHeight.(Ref(font), lines)) + (length(lines) - 1) * font.lineSpacing

    return width, height
end

const stringSurfaceCache = Dict{String, Cairo.CairoSurface}()

function getStringSurface(font::Font, s; replaceMissing::Char=' ', align::Symbol=:left)
    return get!(stringSurfaceCache, s) do
        accX = 0
        accY = 0

        lineAllignFunc = get(stringLineAlignOffsets, align, :left)
        widest, height = size(font, s, replaceMissing)

        surface = CairoARGBSurface(widest, height)
        ctx = getSurfaceContext(surface)

        for line in split(s, "\n")
            width = lineWidth(font, line, replaceMissing)
            accX = lineAllignFunc(widest, width)

            for c in line
                quadX, quadY, width, height = get(font.charDict, c, font.charDict[replaceMissing])
                drawImage(ctx, font.surface, accX, accY, quadX, quadY, width, height)

                accX += width + font.charSpacing
            end

            accY += lineHeight(font, line, replaceMissing) + font.lineSpacing
        end

        return surface
    end
end

function drawString(ctx::Cairo.CairoContext, font::Font, s, x=0, y=0; scale=1.0, replaceMissing::Char=' ', align::Symbol=:left, tint=colors.canvas_font_color)
    shouldSave = x != 0 || y != 0 || scale != 1

    if shouldSave
        Cairo.save(ctx)
        Cairo.translate(ctx, x, y)
        Cairo.scale(ctx, scale, scale)
    end

    surface = getStringSurface(font, s, replaceMissing=replaceMissing, align=align)
    drawImage(ctx, surface, 0, 0, tint=tint)

    if shouldSave
        Cairo.restore(ctx)
    end
end

const pico8FontString = raw"""

 !"#$% ´()*+,-./
0123456789:;<=>?
@abcdefghijklmno
pqrstuvwxyz[\]^_
`ABCDEFGHIJKLMNO
PQRSTUVWXYZ{|}~
"""

const pico8Font = Font(abs"../assets/pico8_font.png", pico8FontString)

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
"""
    writewdir(fun::Function, gl::GitLink; tout = 60.0))

Allow to access the working directory of the GitLink just before a `push event`.
The function `fun(working_dir)` will be executed and it should access
the files into `working_dir`.
Any changes to those files will be commited (please do not mess up the `.git` folder).
It is recommended that `fun` not to be an expensive function.
This method will update all events tokens `pull/stage/push`.
This method will sleep till (or timeout `tout`) the GitLink lock is free (which must by must of the time, but...).
Returns `true` if the action was successful.
"""
function writewdir(fun::Function, gl::GitLink; verbose = false, tout = 60.0)

    wdir = repo_dir(gl)
    mkpath(wdir)

    before_push = () -> fun(wdir)
    return sync_link(gl; before_push, tout, verbose, force = true)
end


text:
using Canopy
using Documenter

makedocs(
    format=Documenter.HTML(prettyurls=false),  # NOTE: prettyurls does not work
    sitename="Canopy.jl documentation",
    pages=[
        "Home" => "index.md",
        "Library" => Any[
            "lib.md",
            "lib/constants.md",
            "lib/water.md",
            "lib/air.md",
            "lib/physchem.md",
            "lib/radtrans.md",
            "lib/transfer.md",
            "lib/leaf.md",
            "lib/watercycle.md",
        ]
    ]
)