text:
export allAtributesToLayout, AtributesToLayout, CopyDict, DomToLayout

# boxes = [ "div", "aside", "audio", "button", "canvas", "form", "header", "img",
#           "input", "ul", "ol", "li", "menu", "menuitem", "progress", "textarea",
#           "table", "tbody", "td", "th", "thead", "tr", "tfoot", "video"]
# textboxes = [  "p", "a", "em", "h1", "h2", "h3", "h4", "h5", "h6",
#           "i", "s", "span", "small", "strong", "title", "abbr", "address",
#           "b", "blockquote", "q", "title", "sub", "sup"]
# ======================================================================================
# Test for string in array (this probably already exists somewhere).
# ======================================================================================
function isAny(key::String, strings::Array{String})
     for str in strings
       if key == str
         return str
       end
     end
     return nothing
end
# ======================================================================================
#
# ======================================================================================
function allAtributesToLayout(document::Page, parent)
    isa(parent.shape, TextLine) && return # TODO: may be able to make this part of multiple dispatch later

    parent.rows  = []  # Clean out old data (?)
        for child in parent.children
                child.parent = parent # attach to parent
                AtributesToLayout(document, child) # Create Child's DOM (DomToLayout.jl)
                allAtributesToLayout(document, child) # Create child's children
        end
end
# ======================================================================================
# Copy an entire nested dictionary to another giving preference to values of "primary"
# ======================================================================================
function CopyDict(primary::Dict, secondary::Dict)
   allKeys = union([key for key in keys(primary)],[key for key in keys(secondary)])
    for key in allKeys
          if haskey(primary, key) && haskey(secondary, key)
                if isa(primary[key], Dict) # its a sub-Dict
                    CopyDict(primary[key], secondary[key])
                end
          elseif haskey(secondary, key) # secondary must have key. Just copy it!
                primary[key] = secondary[key]
          end
    end
end
# ======================================================================================
# Generate a layout Tree from the data stored in the DOM. This does not position the nodes
# in the page but only parses the DOM attributes to the Layout tree.
# ======================================================================================
function CreateShape(form::String, node::Element, h, w)

    if form == "circle"
        node.shape = Circle()
    # elseif form == "path"
    elseif form == "arc"
        node.shape = Arc()
    elseif form == "page" # alternatively I could use the familiar "body" tag for this!
        node.shape = NBox()
          node.shape.color = [1,1,1]
          node.shape.border = Border(0,0, 0,0, 0,0, "solid",[.8,.8,.8,1],[0,0,0,0])
          node.shape.width   = w
          node.shape.height  = h
          node.shape.flags[IsVScroll] = true # IsHScroll, IsVScroll
          node.shape.flags[FixedHeight] = true # Because it's the window!
          node.shape.flags[Clip] = true
    # elseif isAny(tag, boxes) !== nothing
    #             node.shape = NBox()
            # end
    else
        node.shape = NBox()
        if haskey(node.DOM, "text")
            # println(node.DOM["text"])
                node.font = Font()
                textnode = TextElement()
                textnode.shape = TextLine(node, node.DOM["text"])
                push!(node.children, textnode)
                # TODO: set text width as a default
                # something like this:
                    # ctx = CairoContext( CairoRGBSurface(0,0) )
                    # setTextContext(ctx, node)
                    # text_extents(ctx, node.shape.text )[3]
        end
    end
    # if  isAny(form, text) !== nothing # Any from the list
    #     #node.shape = NText()
    #     node.font = Font() # The actual text nodes are created later.
    # end
end
# ======================================================================================
# Parse something like this:   "M-2,2 L8,2.8"   to   [('M', [-2,2]),('L', [8,2.8])]
# WAS: ([A-Za-z])+([\s\d\.,-]*)
# ======================================================================================
function parsePathData(str)
           steps = matchall(r"([a-zA-Z] *(?:(?:-?\d*(?:\.\d+)?(?: *|,?)))*)", str)
           s = []
    for step in steps
        values = matchall(r"[\d\.-]+", step)
        fvals = map(x->parse(Float64,x), values)
        l = step[1]
        push!(s, (l, fvals))
    end
    return s
end
# ======================================================================================
#
# ======================================================================================
function DomToLayout(document::Page, node)

          DOM = node.DOM
          PageStyles = document.styles
          h, w = document.height ,document.width

      if haskey(DOM, "float")
          if DOM["float"] == "right"
              node.shape.flags[FloatRight] = true
          elseif DOM["float"] == "left"
              node.shape.flags[FloatLeft] = true
          end
      end
# "path":"M2,2 L8,2 L2,5 L8,5 L2,8 L8,8"
# b = "M2,2 L-8,2 L2,5 L8,5 L2,8 L8,8"

        # display:inline-block; height:21; padding:3;
        # width:100; color:0.6 0.6 0.6; margin:0 1 0 0;
        # border-radius:7 0 0 3; border-width:0 2 0 0; border-style:solid;
        # border-color:0.2 0.6 0.99;

      if haskey(DOM, "path")
                    node.shape.flags[HasPath] = true
                    properties
      end

      if haskey(DOM, "image")
                    node.shape.flags[HasImage] = true
      end
      if haskey(DOM, "radius")
          r =  DOM["radius"]
                    node.shape.radius = r
                    node.shape.width  = r*2
                    node.shape.height = r*2
      end

      if haskey(DOM, "opacity")
                    node.shape.opacity = DOM["opacity"]
                    node.shape.flags[HasOpacity] = true
      end
      if haskey(DOM, "color")
                    node.shape.color = GetTheColor(node.shape, DOM["color"])
      end

      if haskey(DOM, "linear-gradient")
        grad = DOM["linear-gradient"]
      #  "linear-gradient":{ "start":"left" "stops":[30,60,90], "colors":["red","green","blue"] },
        colorstops = []

              if haskey(grad, "start")
                str = grad["start"]
              end
              if haskey(grad, "colors")
                    steps = length(grad["colors"])
                    stops = []
                    if haskey(grad, "stops")
                        stops = fill(1.0, steps)
                        for j in 1:steps
                           stops[j] = grad["stops"][j]
                        end
                    else # automatically generate stop points
                        X = 1/(steps+1) # y = X/2
                        stops = collect(X:X:1)
                    end

                for i in 1:steps
                  color = GetTheColor(node.shape, grad["colors"][i])
                  if length(color) == 3
                    push!(color, 1.0)
                  end
                 push!(colorstops, [stops[i], color...]   )

                end
              end
              node.shape.gradient = colorstops
                node.shape.flags[LinearGrad] = true
      end
      if haskey(DOM, "radial-gradient")
                node.shape.flags[RadialGrad] = true
      end
      # This is for testing purposes!
      if haskey(DOM, "marked")
              node.shape.flags[Marked] = true
      end
      if haskey(DOM, "top")
              node.shape.top = DOM["top"]
      end
      if haskey(DOM, "left")
              node.shape.left = DOM["left"]
      end
      if haskey(DOM, "bottom")
              node.shape.top = DOM["bottom"]
              node.shape.flags[Bottom] = true
      end
      if haskey(DOM, "right")
              node.shape.left = DOM["right"]
              node.shape.flags[Right] = true
      end
      if haskey(DOM, "width")
              node.shape.width = DOM["width"]
      end
      if haskey(DOM, "height")
              node.shape.height = DOM["height"]
              node.shape.flags[FixedHeight] = true
      end

      if haskey(DOM, "vertical-align")
          if DOM["vertical-align"] == "bottom"
              node.shape.flags[AlignBase] = true
          elseif DOM["vertical-align"] == "middle"
              node.shape.flags[AlignMiddle] = true
          #else node.shape.flags[TextJustify] = true
          end
      end

      if haskey(DOM, "display")
          # like a <span> ...height/width are ignored
         if DOM["display"] == "inline"
              node.shape.width = 0
              node.shape.height = 0
              node.shape.flags[FixedHeight] = false
              node.shape.flags[DisplayInline] = true
         end
         # Like inline but with height/width
         if DOM["display"] == "inline-block"
             node.shape.flags[DisplayInlineBlock] = true
             node.shape.flags[LineBreakAfter] = true
         end
          # LineBreakBefore, LineBreakAfter,
          # set width has it's own row and set width/max-width ...default width is 100%
          if DOM["display"] == "block"
              node.shape.flags[DisplayBlock] = true
              node.shape.flags[LineBreakBefore] = true
              node.shape.flags[LineBreakAfter] = true
          end

          if DOM["display"] == "none"
              node.shape.flags[DisplayNone] = true
          end
          if DOM["display"] == "table"
              node.shape.flags[DisplayTable] = true
          end
         # DisplayBlock, DisplayInlineBlock, DisplayNone, DisplayTable, DisplayFlex,

      end
      #.........................................................................
      if haskey(DOM, "position")
          # like a <span> ...height/width are ignored
         if DOM["position"] == "fixed"
              node.shape.flags[Fixed] = true
         end
         if DOM["position"] == "absolute"
              node.shape.flags[Absolute] = true
              node.parent.shape.flags[HasAbsolute] = true
         end
         if DOM["position"] == "relative"
              node.shape.flags[Relative] = true
         end
     end

       #.........................................................................
#        b = false
#        bcolor = [0,0,0]
#        bradius = [0,0,0,0]
#        bstyle = "solid"
#        bwidth = [0,0,0,0, 0,0]
# if haskey(DOM, "border-color")
#     bcolor = GetTheColor(node.shape, DOM["border-color"]) #DOM["border-color"]
#     b=true
# end
# if haskey(DOM, "border-radius")
#     r = DOM["border-radius"]
#    if isa(r, Array)
#      bradius = r
#    else
#      bradius = [r,r,r,r]
#    end
#    node.shape.flags[IsRoundBox] = true # IsRoundBox
#     b=true
# end
#
# if haskey(DOM, "border-width")
#     w = DOM["border-width"]
#    if isa(w, Array) && length(w) == 4
#      bwidth = [ w[1],w[2],w[3],w[4],  w[1]+w[3], w[2]+w[4] ]
#    elseif isa(w, Number)
#      bwidth = [ w,w,w,w,  w*2, w*2 ]
#    else
#      if w == "thin"
#        bwidth = [1,1,1,1,  2,2]
#      end
#      if w == "medium"
#        bwidth = [3,3,3,3,  6,6]
#      end
#      if w == "thick"
#        bwidth = [4,4,4,4,  8,8]
#      end
#    end
#
#     b=true
# end
# if haskey(DOM, "border-style")
#     bstyle = DOM["border-style"]
#     b=true
# end
# if b == true
#     node.shape.border  = Border(bwidth... , bstyle, bcolor, bradius)
# end

if haskey(DOM, "border")
          border = DOM["border"]

          if haskey(border, "color")
            color = GetTheColor(node.shape, border["color"]) # only works for: [.5,.8,.8]
          end
          if haskey(border, "radius")
               r = border["radius"]
              if isa(r, Array)
                radius = r
            else
                radius = [r,r,r,r]
              end
              node.shape.flags[IsRoundBox] = true # IsRoundBox
          else
              radius = [0,0,0,0]
          end

          if haskey(border, "style")
              style = border["style"]
          else
              style = border["solid"] # default
          end

          if haskey(border, "width")
                  w = border["width"]
                 if isa(w, Array) && length(w) == 4
                   width = [ w[1],w[2],w[3],w[4],  w[1]+w[3], w[2]+w[4] ]
                 elseif isa(w, Number)
                   width = [ w,w,w,w,  w*2, w*2 ]
                 else
                   if w == "thin"
                     width = [1,1,1,1,  2,2]
                   end
                   if w == "medium"
                     width = [3,3,3,3,  6,6]
                   end
                   if w == "thick"
                     width = [4,4,4,4,  8,8]
                   end
                 end
          end
          node.shape.border  = Border(width... , style, color, radius)
      end
      #.........................................................................
      if haskey(DOM, "padding")
        p = DOM["padding"]
        padding = []
               if isa(p, Array)
                     if length(p) == 4
                       padding = [ p[1],p[2],p[3],p[4],  p[1]+p[3], p[2]+p[4] ]
                     end
                     if length(p) == 2
                       padding = [ p[1],p[2],p[1],p[2],  p[1]*2, p[2]*2 ]
                     end
               else
                 padding = [ p,p,p,p,  p*2, p*2 ]
               end
        node.shape.padding = BoxOutline(padding...)
      end
      #.........................................................................
      if haskey(DOM, "margin")
        p = DOM["margin"]
        margin = []
               if isa(p, Array)
                     if length(p) == 4
                       margin = [ p[1],p[2],p[3],p[4],  p[1]+p[3], p[2]+p[4] ]
                     end
                     if length(p) == 2
                       margin = [ p[1],p[2],p[1],p[2],  p[1]*2, p[2]*2 ]
                     end
               else
                 margin = [ p,p,p,p,  p*2, p*2 ]
               end
        node.shape.margin = BoxOutline(margin...)
      end
      #.........................................................................
# font-color:black; font-size:15; align:left; lineHeight:1.4; font-family:sans;

if haskey(DOM, "font-color") || haskey(DOM, "font-size") || haskey(DOM, "font-family")
    node.font = Font()

    if haskey(DOM, "font-color")
        node.font.color = GetTheColor(node.font, DOM["font-color"])
    end
    if haskey(DOM, "font-size")
        node.font.size = DOM["font-size"]
    end
    if haskey(DOM, "font-family")
        node.font.family = DOM["font-family"]
    end
    if haskey(DOM, "lineHeight")
        node.font.lineHeight = DOM["lineHeight"]
    end
    if haskey(DOM, "align")
        if DOM["align"] == "center"
            node.font.flags[TextCenter] = true
        elseif DOM["align"] == "right"
            node.font.flags[TextRight] = true
        end
    end

end
# "font-weight"  "font-size"  "font-family"  "font-style"  "font-color"  "font-gradient"
# "font-opacity"  "font-padding"  "font-border"  "font-margin"  "font-offset"
# "font-lineHeight"  "font-lineWidth"  "font-fill"  "vertical-align"  "align"
      if haskey(DOM, "font") # && isa(node.shape, F o n t)

          if !isdefined(node, :font)
              node.font = Font()
          end
        font = DOM["font"]
# "fill" "linewidth" "color" "align" "vertical-align" "style" "size" "lineHeight" "weight" "family"
        if haskey(font, "fill")
                node.font.fill = GetTheColor(node.shape, font["fill"])
                node.font.flags[TextPath] = true
        end
        if haskey(font, "linewidth")
                node.font.lineWidth = font["linewidth"]
                node.font.flags[TextPath] = true
        end
        if haskey(font, "color")
                node.font.color = GetTheColor(node.font, font["color"])
        end
        if haskey(font, "align")
            if font["align"] == "center"
                node.font.flags[TextCenter] = true
            elseif font["align"] == "right"
                node.font.flags[TextRight] = true
            #else node.shape.flags[TextJustify] = true
            end
        end

        if haskey(font, "vertical-align")
            if font["vertical-align"] == "bottom"
                node.font.flags[AlignBase] = true
            elseif font["vertical-align"] == "middle"
                node.font.flags[AlignMiddle] = true
            #else node.shape.flags[TextJustify] = true
            end
        end
        if haskey(font, "style")
                if font["style"] == "italic"
                    node.font.flags[TextItalic] = true
                elseif font["style"] == "oblique"
                    node.font.flags[TextOblique] = true
                end
        end
        if haskey(font, "size")
                node.font.size = font["size"]
        end
        if haskey(font, "lineHeight")
                node.font.lineHeight = font["lineHeight"]
                #node.font.height = node.shape.size * node.shape.lineHeight
        end
        if haskey(font, "weight")
          if font["weight"] == "bold"
              node.font.flags[TextBold] = true
          end
        end
        if haskey(font, "family")
                node.font.family = font["family"]
        end
      end

#...............Event Registration..................
      if haskey(DOM, "hover")
          #  = DOM["hover"]
      end
end # function

# ======================================================================================
#
# ======================================================================================
function AtributesToLayout(document::Page, node, generative::Bool=true)

    DOM = node.DOM # PageStyles = document.styles
    h, w = document.height, document.width

    if haskey(DOM, ">")
        tag = DOM[">"]

        # Attach events:
        # TODO: expound and move to separate function.
        if haskey(DOM, "mousedown")
            push!(document.eventsList.mousedown, node)
        end
        if haskey(DOM, "click")
            push!(document.eventsList.click, node)
        end
        if haskey(DOM, "hover")
            push!(document.eventsList.mouseover, node)
            push!(document.eventsList.mouseout, node)
        end
        # Create actual node.
        if generative
            if isa(tag, Array)
                println(DOM)
            else
                CreateShape(tag, node, h, w)
            end
        end

        # Add user defined styles to node's DOM.
        # if haskey(DOM, "style")
        #     style = DOM["style"]
        #     CopyDict(DOM, PageStyles[style])
        # end

        # Add default values generally associated with this tag
        CopyDict(DOM, Tags_Default[tag])
        DomToLayout(document, node)
    end
end
# ======================================================================================
#
# ======================================================================================

# ======================================================================================
# Dict{Any,Any}("nodes"=>Any[],">"=>Any["h", 1.0],"font"=>Dict{Any,Any}("nodes"=>Any[],"color"=>"white","align"=>"center"),"text"=>"Why reinvent the wheel?")


text:
"""
A complete statistical model for one river, with three-stage modelling: yearly average, seasonal averages, daily averages.
"""
immutable Adam2015RiverStatisticalModel
  year::GeneralizedExtremeValue # Annual average; Weibull when xi = -1.
  season_wet::GeneralizedExtremeValue # Difference between seasonal and annual averages, normalised wrt annual average.
  season_dry::GeneralizedExtremeValue
  daily_wet::GeneralizedExtremeValue
  daily_dry::GeneralizedExtremeValue

  seasons_separator::SeasonsSeparator
end

getYearlyAverage(m::Adam2015RiverStatisticalModel) = m.year
function getSeasonAverageDelta(m::Adam2015RiverStatisticalModel, season::Symbol)
  if season == :Wet
    return m.season_wet
  elseif season == :Dry
    return m.season_dry
  else
    error("Unknown season: " * string(season))
  end
end
function getSeasonDailyDelta(m::Adam2015RiverStatisticalModel, season::Symbol)
  if season == :Wet
    return m.daily_wet
  elseif season == :Dry
    return m.daily_dry
  else
    error("Unknown season: " * string(season))
  end
end
getSeasonSeparator(m::Adam2015RiverStatisticalModel) = m.seasons_separator

"""
Based on the data contained in `TS` (one time series per year!), fit a statistical model.

  * `TS`: array of time series (one series per year) to fit the statistical model on.
  * `verbosity`: level of output.
    * `0`: no output.
    * `1`: summary output
    * `2`: detailed output of this algorithm's details
    * `3`: detailed output of this algorithm's and its sub-algorithms' details
    * `4`: including solver output

TODO: implement flooding.
"""
function fit_Adam2015(TS::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1};
                      seasonsSeparator::SeasonsSeparator=HydrologicalWesternEurope,
                      verbosity::Int=1)
  ## Separate every year in two seasons: dry and wet.
  sep_ts = [separate(ts, seasonsSeparator, guaranteeTimeStep=true) for ts in TS]
  wet_TS = [ts[:Wet] for ts in sep_ts]
  dry_TS = [ts[:Dry] for ts in sep_ts]

  # Compute a few useful averages within those seasons.
  ann_avg = [mean(values(ts)) for ts in TS]
  wet_avg = [mean(values(ts)) for ts in wet_TS]
  dry_avg = [mean(values(ts)) for ts in dry_TS]

  if verbosity >= 1
    println(" == ")
    println(" ==== Statistical model fitting on ", length(TS), " years")

    if verbosity >= 2
      println(" ====== Length of wet season: ", length(wet_TS[1]))
      println(" ====== Length of dry season: ", length(dry_TS[1]))
    end
  end

  ## First deal with annual average.
  if verbosity >= 1
    println(" ==== Fitting a generalised extreme values model on the annual average.")
  end

  d_year_params = fit_mle_optim(ExtremeValueDistributions.GeneralizedExtremeValue, ann_avg, [0.5, 0.5, -1], verbose=(verbosity >= 3))
  d_year = Distributions.GeneralizedExtremeValue(d_year_params[1], exp(d_year_params[2]), d_year_params[3]) # ~ <Q>_Y

  if verbosity >= 2
    println(" ====== Annual average: generalised extreme values.")
    println(" ======    Parameters: μ = ", d_year.μ, "; σ = ", d_year.σ, "; ξ = ", d_year.ξ, ".")
  end

  ## Second with the seasonal variations with respect to the annual averages.
  if verbosity >= 1
    println(" ==== Fitting generalised extreme values models on the seasonal averages.")
  end

  wet_deltas = (wet_avg .- ann_avg) ./ ann_avg
  d_season_wet_params = fit_mle_optim(ExtremeValueDistributions.GeneralizedExtremeValue, wet_deltas, [0.5, 0.5, 0.5], verbose=(verbosity >= 3))
  d_season_wet = Distributions.GeneralizedExtremeValue(d_season_wet_params[1], exp(d_season_wet_params[2]), d_season_wet_params[3]) # ~ \Delta Q_S^W

  if verbosity >= 2
    println(" ====== Wet season average: generalised extreme values.")
    println(" ======    Parameters: μ = ", d_season_wet.μ, "; σ = ", d_season_wet.σ, "; ξ = ", d_season_wet.ξ, ".")
  end

  dry_deltas = (dry_avg .- ann_avg) ./ ann_avg
  d_season_dry_params = fit_mle_optim(ExtremeValueDistributions.GeneralizedExtremeValue, dry_deltas, [0.5, 0.5, 0.5], verbose=(verbosity >= 3))
  d_season_dry = Distributions.GeneralizedExtremeValue(d_season_dry_params[1], exp(d_season_dry_params[2]), d_season_dry_params[3]) # ~ \Delta Q_S^D

  if verbosity >= 2
    println(" ====== Dry season average: generalised extreme values.")
    println(" ======    Parameters: μ = ", d_season_dry.μ, "; σ = ", d_season_dry.σ, "; ξ = ", d_season_dry.ξ, ".")
  end

  ## Finally with the daily variations within each season.
  # Here, working with matrices of values, as each day is considered (and not a single point per year).
  if verbosity >= 1
    println(" ==== Fitting generalised extreme values models on the daily inflow.")
  end

  wet_daily_deltas = vcat([(values(wet_TS[i]) - wet_avg[i]) / wet_avg[i] for i in 1:length(wet_TS)]...)
  d_daily_wet_params = fit_mle_optim(ExtremeValueDistributions.GeneralizedExtremeValue, wet_daily_deltas, [0.5, 0.5, 0.5], verbose=(verbosity >= 3))
  d_daily_wet = Distributions.GeneralizedExtremeValue(d_daily_wet_params[1], exp(d_daily_wet_params[2]), d_daily_wet_params[3]) # ~ \Delta Q_d^W

  if verbosity >= 2
    println(" ====== Daily variations during the wet season: generalised extreme values.")
    println(" ======    Parameters: μ = ", d_daily_wet.μ, "; σ = ", d_daily_wet.σ, "; ξ = ", d_daily_wet.ξ, ".")
  end

  dry_daily_deltas = vcat([(values(dry_TS[i]) - dry_avg[i]) / dry_avg[i] for i in 1:length(dry_TS)]...)
  d_daily_dry_params = fit_mle_optim(ExtremeValueDistributions.GeneralizedExtremeValue, dry_daily_deltas, [0.5, 0.5, 0.5], verbose=(verbosity >= 3))
  d_daily_dry = Distributions.GeneralizedExtremeValue(d_daily_dry_params[1], exp(d_daily_dry_params[2]), d_daily_dry_params[3]) # ~ \Delta Q_d^D

  if verbosity >= 2
    println(" ====== Daily variations during the wet season: generalised extreme values.")
    println(" ======    Parameters: μ = ", d_daily_dry.μ, "; σ = ", d_daily_dry.σ, "; ξ = ", d_daily_dry.ξ, ".")
  end

  ## Done. Prepare the output data structure.
  if verbosity >= 1
    println(" ==== Model fitted. ")
    println(" == ")
  end

  return Adam2015RiverStatisticalModel(d_year, d_season_wet, d_season_dry, d_daily_wet, d_daily_dry, seasonsSeparator)
end

function fit_Adam2015(r::River; kwargs...)
  if haskey(r.ext, "adam2015")
    return r.ext["adam2015"]
  else
    m = fit_Adam2015(r.scenarios; kwargs...)
    r.ext["adam2015"] = m
    return m
  end
end

"""
Prints a report analysing the statistical model with observations made on the raw data.
"""
function printReport(r::River, stat::Adam2015RiverStatisticalModel)
  sc = r.scenarios
  wet = separate(sc, stat.seasons_separator, :Wet)
  dry = separate(sc, stat.seasons_separator, :Dry)
  wet_avg = Float64[mean(values(ts)) for ts in wet]
  dry_avg = Float64[mean(values(ts)) for ts in dry]

  ymean_sc = Float64[mean(values(ts)) for ts in sc]
  wet_smean_sc = (wet_avg .- ymean_sc) ./ ymean_sc
  dry_smean_sc = (dry_avg .- ymean_sc) ./ ymean_sc
  wet_dmean_sc = vcat([(values(wet[i]) .- wet_avg[i]) ./ wet_avg[i] for i in 1:length(wet)]...)
  dry_dmean_sc = vcat([(values(dry[i]) .- dry_avg[i]) ./ dry_avg[i] for i in 1:length(dry)]...)

  # Printing helpers.
  n2s = n -> @sprintf("%7.3f", n)
  function printReport_distr(title::AbstractString, distr::ContinuousUnivariateDistribution, data::Array{Float64, 1})
    # Expected output:
	# 	 ++++ Yearly averages: mean 0.511 vs 0.509
	# 	 ++++              variance 0.014 vs 0.014
	# 	 ++++              kurtosis 0.192 vs 1.689
	#                      ^^^^ (4 characters)
	#         ^^^^^^^^^^^^^^^   length(title), followed by " :"

    title_blank = " " ^ (length(title) - 4)

    println(" ++++ " * title       *     ": mean " * n2s(mean(distr))      * " vs " * n2s(mean(data)))
    println(" ++++ " * title_blank * "  variance " * n2s(var(distr))       * " vs " * n2s(var(data)))
    println(" ++++ " * title_blank * "  skewness " * n2s(skewness(distr))  * " vs " * n2s(skewness(data)))
    println(" ++++ " * title_blank * "  kurtosis " * n2s(kurtosis(distr))  * " vs " * n2s(kurtosis(data)))
  end

  # Validity report.
  println(" ++ Validity report for the statistical model (model vs reality): ")
  printReport_distr("Yearly averages",   stat.year,       ymean_sc)
  printReport_distr("Wet season deltas", stat.season_wet, wet_smean_sc)
  printReport_distr("Dry season deltas", stat.season_dry, dry_smean_sc)
  printReport_distr("Wet daily deltas",  stat.daily_wet,  wet_dmean_sc)
  printReport_distr("Dry daily deltas",  stat.daily_dry,  dry_dmean_sc)

  # Covariance report.
  corr_adj  = c -> (abs(c) < 0.3) ? "weak" : (abs(c) < 0.8) ? "moderate" : "strong"
  corr_sign = c -> (c <= 0.0) ? "negative" : "positive"

  pearson_ymean_wetsmean = cor(ymean_sc, wet_smean_sc)
  pearson_ymean_drysmean = cor(ymean_sc, dry_smean_sc)

  println(" ++ ")
  println(" ++ Covariance report for the statistical model: ")
  println(" ++++ Pearson's coefficient between yearly average and wet season average: " * n2s(pearson_ymean_wetsmean))
  println(" ++++      " * corr_adj(pearson_ymean_wetsmean) * " " * corr_sign(pearson_ymean_wetsmean) * " correlation")
  println(" ++++ Pearson's coefficient between yearly average and dry season average: " * n2s(pearson_ymean_drysmean))
  println(" ++++      " * corr_adj(pearson_ymean_drysmean) * " " * corr_sign(pearson_ymean_drysmean) * " correlation")
end

# Base functions: from the statistical model to actual scenarios.
# Those functions allow scalar *and* vector arguments! (Only useful for Qd, however.)
# 213 days in wet season:     Date(2016, April, 30) - Date(2015, October, 1) + Day(1)
# 152 left in dry season.
adam2015_wet_season_value(QY, QS) = QY * (1 + QS)
adam2015_dry_season_value(QY, QS) = QY * (1 - QS * 213 / 152)
adam2015_wet_day_value(QY, QS, Qd) = QY * (1 + QS) * (1 + Qd)
adam2015_dry_day_value(QY, QS, Qd) = QY * (1 - QS * 213 / 152) * (1 + Qd)

"""
Based on the given statistical `model`, generate one yearly scenarios, with a value each `period`.
This scenario starts at the beginning of the wet season (as per the seasons in `model`), and occur in the year
`currentYear`.
"""
function generateScenarios_Adam2015(model::Adam2015RiverStatisticalModel, period::Period, currentYear::Int; kwargs...)
  # Pick a year average and a season average at random. The wet average has a high likelihood of being negative,
  # so generate until something positive is generated.
  year = rand(getYearlyAverage(model))

  # dry_season = rand(getSeasonAverageDelta(model, :Dry)) # Not fitted, using a fraction of the other season.
  wet_season = rand(getSeasonAverageDelta(model, :Wet))
  while adam2015_wet_season_value(year, wet_season) <= 0.0
    wet_season = rand(getSeasonAverageDelta(model, :Wet))
  end

  # Then work on days.
  bs, es = seasonBoundariesGuarantee(getSeasonSeparator(model), period, currentYear)
  n_wet = length(bs[:Wet] : period : es[:Wet])
  n_dry = length(bs[:Dry] : period : es[:Dry])

  gen_wet = adam2015_wet_day_value(year, wet_season, rand(getSeasonDailyDelta(model, :Wet), n_wet))
  gen_dry = adam2015_dry_day_value(year, wet_season, rand(getSeasonDailyDelta(model, :Dry), n_dry))

  gen = vcat(gen_wet, gen_dry)
  gen[gen .<= 0] = 0.
  return (collect(bs[:Wet] : period : es[:Dry]), gen)
end

generateScenarios_Adam2015(r::River, period::Period, currentYear::Int; kwargs...) =
  generateScenarios_Adam2015(fit_Adam2015(r; kwargs...), period, currentYear)

"""
Computes confidence intervals based on the statistical model.

Assumption: all three stages of probability distributions are not correlated.
This assumption seems mostly valid, according to a few correlation tests.
"""
function computeCI_Adam2015(r::River, model::Adam2015RiverStatisticalModel, level::Float64=.95)
  TS = r.scenarios

  ## Start by computing the actual values to this confidence level using quantiles.
  q_low = (1 - level) / 2
  q_high = 1 - q_low

  # First the annual confidence interval.
  year_low = quantile(model.year, q_low)
  year_mean = mean(model.year)
  year_high = quantile(model.year, q_high)

  # Then the same for the seasonal average normalised deltas.
  wet_season_low = quantile(model.season_wet, q_low)
  wet_season_mean = mean(model.season_wet)
  wet_season_high = quantile(model.season_wet, q_high)

  dry_season_low = quantile(model.season_dry, q_low)
  dry_season_mean = mean(model.season_dry)
  dry_season_high = quantile(model.season_dry, q_high)

  # And finally for the seasonal daily normalised deltas.
  wet_daily_low = quantile(model.daily_wet, q_low)
  wet_daily_mean = mean(model.daily_wet)
  wet_daily_high = quantile(model.daily_wet, q_high)

  dry_daily_low = quantile(model.daily_dry, q_low)
  dry_daily_mean = mean(model.daily_dry)
  dry_daily_high = quantile(model.daily_dry, q_high)

  ## Reconstruct the complete confidence interval as the sum of the previous random variables.
  # Here is used the hypothesis of absence of correlation.
  # Technique:
  #   - start from the annual average bound
  #   - transform one season's distribution based on this to get the seasonal average's bound
  #   - sum those bounds to get the total bound for seasonal variations
  # On average, the daily variations are zero: they are not taken into account when computing the mean.
  # However, these have a larger impact on the bounds of the confidence interval.
  dwet_low  = adam2015_wet_day_value(year_low,  wet_season_low,  wet_daily_low)
  dwet_high = adam2015_wet_day_value(year_high, wet_season_high, wet_daily_high)
  ddry_low  = adam2015_dry_day_value(year_low,  dry_season_low,  dry_daily_low)
  ddry_high = adam2015_dry_day_value(year_high, dry_season_high, dry_daily_high)

  # NOTE: dry_season seems to be replaced by wet_season in Adam's code!?
  dwet_mean = adam2015_wet_season_value(year_mean, wet_season_mean)
  ddry_mean = adam2015_dry_season_value(year_mean, wet_season_mean)

  ## Prepare the outputs, respecting the two seasons.
  wet_begin = begins(model.seasons_separator, :Wet, 1)
  wet_end = ends(model.seasons_separator, :Dry, 1)
  dry_begin = begins(model.seasons_separator, :Dry, 1)
  dry_end = ends(model.seasons_separator, :Dry, 1)
  period = getPeriod(r)

  # The length of the reconstructed series must match that of the time series!
  while length(vcat(wet_begin:period:wet_end, dry_begin:period:dry_end)) > length(timestamp(TS[1]))
    wet_begin += period
  end
  if length(vcat(wet_begin:period:wet_end, dry_begin:period:dry_end)) < length(timestamp(TS[1]))
    error("Reconstructed period too small (" * length(vcat(wet_begin:period:wet_end, dry_begin:period:dry_end)) * ") " *
          "with respect to existing time series (" * length(timestamp(TS[1])) * ").")
  end

  means   = vcat([dwet_mean for d in wet_begin:period:wet_end], [ddry_mean for d in dry_begin:period:dry_end])
  ci_low  = vcat([dwet_low  for d in wet_begin:period:wet_end], [ddry_low  for d in dry_begin:period:dry_end])
  ci_high = vcat([dwet_high for d in wet_begin:period:wet_end], [ddry_high for d in dry_begin:period:dry_end])

  means[means .<= 0.0] = 0.0
  ci_low[ci_low .<= 0.0] = 0.0
  ci_high[ci_high .<= 0.0] = 0.0

  ts_means   = TimeArray(timestamp(TS[1]), means,   colnames(TS[1]), meta(TS[1]))
  ts_ci_low  = TimeArray(timestamp(TS[1]), ci_low,  colnames(TS[1]), meta(TS[1]))
  ts_ci_high = TimeArray(timestamp(TS[1]), ci_high, colnames(TS[1]), meta(TS[1]))

  return ts_means, ts_ci_low, ts_ci_high
end

computeCI_Adam2015(r::River, level::Float64=.95; kwargs...) =
  computeCI_Adam2015(r, fit_Adam2015(r; kwargs...), level)

ext: jl


text:
#=
   The Exponential Functions `exp(x)` and `exp(z)`.
   The associaed functions `expm1(x)` and `expm1(z)`.
=#

const Ln2_D64 = Double64((0.6931471805599453, 2.3190468138462996e-17))

function Base.exp(x::Double64)
   exp_specialvalues(Hi(x))
   # 4.3108e-78 <= absx < 709.9
   absx = abs(x)

   # Solve for k satisfying x = 2^k * log(2) r
   #   with |r| <= 1/2 log(2)
   k = ceil(Int64, (absx / Ln2) - 0.5)
   # p = 2^k
   p = Int64(1) << k
   # Determine r using the k we computed
   r = absx - (k * Ln2)

end

function exp_specialvalues(x::Float64)
   isnan(x) && return x
   absx = abs(x)
   absx > 709.0 && return signbit(x) ? -InfD64 : InfD64
   abx < 2.0^-257 && return one(Double64)
   return nothing
end

#=
   source: Efficient implementation of elementary functions in the medium-precision range
           by Fredrik Johansson
           section 3: Argument reduction

   result = exp(x)

   m = fld(x, log(2))
   t = x - m * log(2)

  julia> FF4(log(BigFloat(2)))
  (0.6931471805599453, 2.3190468138462996e-17, 5.707708438416212e-34, -3.5824322106018114e-50)
  julia> FF4(exp(BigFloat(1)))
  (2.718281828459045, 1.4456468917292502e-16, -2.1277171080381768e-33, 1.5156301598412191e-49)
  julia> FF4(exp(BigFloat(-1)))
  (0.36787944117144233, -1.2428753672788363e-17, -5.830044851072742e-34, -2.8267977849017436e-50)
  log2(exp(1))
  (1.4426950408889634, 2.0355273740931033e-17, -1.0614659956117258e-33, -1.3836716780181402e-50)

   exp(x) = exp(t) * (2.0^m)
            where t  in [0, log(2)) # log(2) ~0.6931471805599453

   exp(t) = (exp(t/(2.0^r)))^(2.0^r)
            where t/(2.0^r) in [0, 2.0^(-r))   # at the expense of r squarings

julia> 2.0^(-8) = 0.00390625
julia> 2.0^(-9) = 0.001953125

   if we precompute exp(i/(2^r)) for i = 0..(2^r - 1)
      we can write exp(x) = exp(x - i/2^r) * exp(i/2^r)  wher i = floor(x * 2^r)
      this achieves r halvings of argument reduction for the cost of a multiplication

sum([x^k / factorial(k) for k=0:n]) # perhaps 11
ex2(x,n) = sum([(x^(2*k) * (1 + 2*k +x)) / factorial(1+2*k) for k=0:n])  # 5 perhaps 6

=#
#=
ir(r) = collect(0:(2^r - 1));
ir2(r1,r2) = collect((2^r1):(2^r2 - 1))
exps(r) = exp.(ir(r) / big"2"^r)
exps2(r1,r2) = exp.(ir2(r1,r2) / big"2"^r2)

expi1 = exps(8);
expi2 = exps2(8,9);
expi1s = expi1[1:128]; # length 128
expi2s = expi2; # leingth 256

expi1s[1],expi1s[end]
(1.0, 1.642293515618909941146175177943829065898476010090130767464052857406618001282157)
expi2s[1],expi2s[end]
(1.648721270700128146848650787814163571653776100710148011575079311640661021194214, 2.712977865600149830047947983547967696224441828481010944922219814560357243651112)

# 256 values
expi = [
big"1.0",
big"1.0039138893383475734436096039034602818988136",
big"1.0078430972064479776934535597601235791933921",
big"1.0117876835593314915141138433716067245783427",
big"1.0157477085866857474585350720823517489067163",
big"1.0197232327137741547831816318953623989389768",
big"1.0237143166023579169688505321657689580745566",
big"1.0277210211516216579160287393882763310060287",
big"1.0317434074991026709387478152815071441944983",
big"1.0357815370216238047361682727501377703668922",
big"1.0398354713362300005766220098751917925969306",
big"1.0439052723011284949845543336459217728408202",
big"1.0479910020166327022767382918196849669390532",
big"1.0520927228261097913502841465226716758711456",
big"1.0562104973169319711813367037681297731537515",
big"1.0603443883214314995499437191070563220034387",
big"1.0644944589178594295633905946428896731007254",
big"1.0686607724313481086073309335198834903610129",
big"1.072843392434877444411300095010813766928624",
big"1.077042382750244952972680570165861926965078",
big"1.0812578074490396031408946517428286548338917",
big"1.0854897308536194727215323833540410098046093",
big"1.0897382175380932310182820254185122044095245",
big"1.0940033323293054627889171570399958544544903",
big"1.0982851403078258486502099342567867188098263",
big"1.1025837068089422170254848406286398617008396",
big"1.106899097423657482787602393741778963630765",
big"1.1112313779996904878094686123406995149077797",
big"1.1155806146424807586937045083794438798390782",
big"1.1199468737161971970128813548111838167408411",
big"1.1243302218447507174517329054449565537714046",
big"1.1287307259128108493029960233035576751121704",
big"1.1331484530668263168290072278117938725655031",
big"1.1375834707160496140618954216220641739226107",
big"1.1420358465335655896761614313488004988947112",
big"1.1465056484573240576286239229647318858420707",
big"1.1509929446911764493221396657803508410236923",
big"1.1554978037059165231111749561077458362413083",
big"1.1600202942403251470292152138940958101001022",
big"1.1645604853022191706801522774745177736791979",
big"1.1691184461695044022981846915118914445916047",
big"1.1736942463912327070434062632478574421242414",
big"1.1782879557886632426631433081958559090232306",
big"1.1828996444563278487122322779162386958636926",
big"1.1875293827631006055888078228331305661056264",
big"1.1921772413532715797057977592610942010914603",
big"1.1968432911476247711821968527700683464587539",
big"1.201527603344520280502316774392647998778235",
big"1.2062302494209807106555860104464335480403936",
big"1.2109513011337818213341018932134198490219412",
big"1.215690830520547451830018254540396832863575",
big"1.2204489099008487293399874776109514959456667",
big"1.2252256118773075794492659275775547905529261",
big"1.2300210093367045556337378772080063588467264",
big"1.2348351754510910046840161108760027424078237",
big"1.2396681836789055850219383947130699551294591",
big"1.2445201077660951549461989520766544581233531",
big"1.2493910217472400479105339951661133401110873",
big"1.2542809999466837520048212530906856313930507",
big"1.2591901169796670108766563253994731033558361",
big"1.2641184477534663633984346034590176916127065",
big"1.2690660674685371394526974695586839474881616",
big"1.2740330516196609292764965387163400150598443",
big"1.2790194759970975438737908883703017739782832",
big"1.2840254166877414840734205680624364583362809",
big"1.2890509500762829358789962405026343223583583",
big"1.2940961528463733098261106268100709645926577",
big"1.299161101981795342131613566113939048057404",
big"1.3042458747676377754893000110524664638914931",
big"1.3093505487914746374362402281617419935314083",
big"1.314475201944549134284134921979872958797392",
big"1.3196199124229621786805060243653563910202619",
big"1.3247847587288655689352375606889644453814644",
big"1.3299698196716598383189614010596868898353478",
big"1.335175174369196792611040910663613496051334",
big"1.3404009022489867542464426159453927247594736",
big"1.345647083049410531482603094316017723141301",
big"1.3509137968209361310794964701714017243589766",
big"1.3562011239273402330584882593857570013204578",
big"1.3615091450469344461782249525402249833846504",
big"1.3668379411737963628387567727212086721727333",
big"1.3721875936190054321983245998746874054136052",
big"1.3775581840118836703607622378892548520433121",
big"1.3829497943012412265652731346321226949971137",
big"1.3883625067566268243844374753000036641797537",
big"1.3937964039695830970106923862996621838753668",
big"1.3992515688549068357862059454951895801349576",
big"1.4047280846519141712060359345473907335921464",
big"1.4102260349257107056997279341815912975975336",
big"1.4157455035684666175720656019875394367736763",
big"1.4212865748006967555595399377075767858324404",
big"1.4268493331725457435352551903480186121366873",
big"1.4324338635650781149714379568249381105086462",
big"1.4380402511915734968454641297483462771412234",
big"1.4436685815988268627523668434427752933667452",
big"1.4493189406684538750641386180656840287525084",
big"1.454991414618201336053793691987518508346842",
big"1.4606860900032627679801132473277319831728503",
big"1.4664030537175991422072580622261772341235882",
big"1.4721423929952647775120012602512819589278704",
big"1.4779041954117384278102094711614577636044574",
big"1.4836885488852595796133850715926159360653458",
big"1.4894955416781699796055764474013653746366884",
big"1.4953252623982604128107686243383379457966073",
big"1.5011778000001227519009843681945126839847195",
big"1.5070532437865072982757571820568627066188943",
big"1.5129516834096854356243837541143330540396894",
big"1.518873208872817616763425566660741826592259",
big"1.5248179105313267046233088024596115443256307",
big"1.5307858790942766883395696202202495812292224",
big"1.5367772056257567954863095631977568976482128",
big"1.5427919815462710215717645653643101046597946",
big"1.5488302986341330979985519845954923375583037",
big"1.5548922490268669197741445831814759100397908",
big"1.5609779252226124543404296589529359617213441",
big"1.5670874200815371529748468766060916540358334",
big"1.5732208268272528862995610345959599306077298",
big"1.5793782390482384255194173094123070388118751",
big"1.5855597506992674910940477323437952296789346",
big"1.5917654561028423906344500651913320024805135",
big"1.5979954499506332678996451470842104572159357",
big"1.6042498273049229848546374859283226869322636",
big"1.6105286836000576588368576715138990015137373",
big"1.6168321146439028769645554044834985475711929",
big"1.6231602166193056100072398828952478310491124",
big"1.6295130860855618480252312878508519071585423",
big"1.6358908199798899801726944886054803066336858",
big"1.642293515618909941146175177943829065898476",
big"1.6487212707001281468486507878141635716537761",
big"1.6551741833034282419274450659151529161544963",
big"1.661652351892567681933037464039017553408944",
big"1.6681558753166801729348278516747020845493453",
big"1.6746848528117839915192948821484471915616397",
big"1.6812393840022962081857139672252104808506591",
big"1.6878195689025528372446796287495748265689688",
big"1.694425507918334936415108367064433784118984",
big"1.7010573018484006794061834954848285349306956",
big"1.7077150518860234248618440232984793229206938",
big"1.7143988596205358051369170173112717453784532",
big"1.7211088270388798584658483299957258431122154",
big"1.7278450565271632281772015524381343960105588",
big"1.7346076508722214526996709396061731128182903",
big"1.7413967132631863701982922765711232180120641",
big"1.7482123472930606617728376253439213210762717",
big"1.7550546569602985572440470365989676887382375",
big"1.7619237466703927276473840580130012741833298",
big"1.768819721237467388648403656081329732171722",
big"1.7757426858858776391885924315490433000564346",
big"1.7826927462518150597656831961070971624870767",
big"1.7896700083849195948479605415087339387426964",
big"1.7966745787498977440179624291990663920674301",
big"1.8037065642281470865372465220624332490665196",
big"1.8107660721193871641205304388476395701804874",
big"1.8178532101432967468045338108486506072754985",
big"1.8249680864411575068942484398874703584971162",
big"1.8321108095775041260671424827506305635923199",
big"1.8392814885417808608139669119715291919683603",
big"1.8464802327500045914933790240679012662617646",
big"1.8537071520464343803765299877394712853248167",
big"1.8609623567052475641570828556931428635577371",
big"1.8682459574322224065018356201881044531149723",
big"1.8755580653664273363172222955086314557144833",
big"1.8828987920819167975074551877605054582064688",
big"1.890268249589433736100954997008217074841387",
big"1.8976665503381187507229937281862390730851715",
big"1.9050938072172259324941501108153118386642345",
big"1.9125501335578454205362498947690069141803877",
big"1.9200356431346326993699355576741048656485595",
big"1.9275504501675446645908831925937307660930033",
big"1.9350946693235824833149602121576936617744414",
big"1.9426684157185412759862975831367569870468297",
big"1.9502718049187666462463361756410743039045826",
big"1.9579049529429180856664000618841160115832721",
big"1.965567976263739280251251825234587055462156",
big"1.9732609918098353457263977417630639167722963",
big"1.980984116967457018727635680625228052334163",
big"1.9887374695822918311174773496469253668482552",
big"1.9965211679612622947596307079391523309846108",
big"2.0043353308743311241896996041635187235500792",
big"2.0121800775563135247276476094443781622950744",
big"2.0200555277086965736853832365247240976992445",
big"2.0279618015014657224310559166446580204412786",
big"2.0358990195749384471803078942778508475105447",
big"2.0438673030416050764938082552508055357536932",
big"2.0518667734879768235699032903388202787086325",
big"2.0598975529764410515311539851427396253726037",
big"2.0679597640471238000138982954032159493019471",
big"2.076053529719759601480774698971190974886495",
big"2.0841789734955686157873760020828031213826762",
big"2.0923362193591411116458702156591367747823646",
big"2.1005253917803293237405302109843070944234469",
big"2.1087466157161467143626575238964185488599182",
big"2.1170000166126746685453698198370956101344916",
big"2.1252857204069766517921478827499607656345219",
big"2.1336038535290198596069082796475876104173024",
big"2.141954542903604388147683817205012365763677",
big"2.1503379159522999554407572913112531245709817",
big"2.1587541005953902027073065862709664584426364",
big"2.1672032252538246054702826653026868576426094",
big"2.1756854188511780242253581730150091331184588",
big"2.1842008108156179245763550154521229783421722",
big"2.1927495310818792968525861728144356991424522",
big"2.2013317100932473053430319197695160332670769",
big"2.2099474788035476974002153710531479045782357",
big"2.2185969686791450027850486355161375002671637",
big"2.2272803117009485537427906566397217584734165",
big"2.2359976403664263564195928556324133682743723",
big"2.2447490876916268443489107954443761003418463",
big"2.2535347872132085448573310784292862413086861",
big"2.2623548729904776893601044121401684567152138",
big"2.2712094796074337986378900692397035812348409",
big"2.2800987421748232743079056783226487661726783",
big"2.2890227963322010278248412694130178222761383",
big"2.2979817782500001784695396251036787167203531",
big"2.3069758246316098519065681271390235943873022",
big"2.3160050727154611110154123174147474653913428",
big"2.32506966027712105082411019794800030841958",
big"2.3341697256313950894987207698035707748063283",
big"2.3433054076344374874670823571053781999975876",
big"2.352476845685870126880867787415201591593135",
big"2.3616841797309095837459864196889344674000366",
big"2.3709275502625025251779192489894964270765329",
big"2.3802070983234694643656048039370752132460841",
big"2.3895229655086569059550222268607428894033582",
big"2.3988752939670979146916457337162967482854659",
big"2.4082642264041811402894735446052064189541318",
big"2.4176899060838283316233663173665756270227349",
big"2.4271524768306803734709670750917334949139414",
big"2.4366520830322918791605185701090931800750108",
big"2.4461888696413343726114469563210060388105308",
big"2.4557629821778080933856445408266406336725324",
big"2.4653745667312624584989612544143600966570358",
big"2.4750237699630252148745063264598957734822349",
big"2.4847107391084403164519704885869968967885934",
big"2.4944356219791145601003068865880476028893769",
big"2.5041985669651730146147577829358268760753641",
big"2.5139997230375232772133861220671538815112637",
big"2.5238392397501285920829681547893051234396304",
big"2.533717267242289865659327631968429446833967",
big"2.5436339562409366134629456444676263089122303",
big"2.5535894580629268734469650775485663612146274",
big"2.5635839246173561209515269431795734112490795",
big"2.5736175084078752204957296405891146914479675",
big"2.5836903625350174497763935698217882499996026",
big"2.5938026406985346313812445890522234238021264",
big"2.6039544971997424078631026762839444413754604",
big"2.6141460869438746959611789503655914445178259",
big"2.6243775654424473558956470538559895767179813",
big"2.6346490888156311118022659383545698498240101",
big"2.6449608137946337595149924670414734811886784",
big"2.6553128977240916980462361132966847731583179",
big"2.6657054985644708212566765507744459129605862",
big"2.676138774894476806349390270046054291415995",
big"2.6866128859134748359664166992354320756243384",
big"2.6971279914439187908098398388318050409222763",
big"2.7076842519337899498539703405097468565320357",
];

julia> function accurateC1C2(C,p=53)
         R3 = FF3(1/C)
         R = sum(BigFloat.(R3))
         C1a = FF3(1/R)
         C1b = map(x->round(x, sigdigits=p-2, base=2), C1a)
         C1 = FF3(sum(BigFloat.(C1b))); C13 = sum(BigFloat.(C1))
         C2a = ((C-C13) / (8*eps(eps(C13)))) * (8*eps(eps(C13)))
         C2 = FF3(C2a)
         return C1, C2
       end
accurateC1C2 (generic function with 2 methods)

julia> accurateC1C2(g)
((1.618033988749895, -5.4321152036825055e-17, 2.6543252083815655e-33), (-6.162975822039155e-33, -3.304991997502108e-50, 6.608118949502031e-67))

julia> function accurateargreduction(x, z, C1, C2)
           u = x - z*C1[1]
           v1 = u - z*C2[1]
           p1,p2 = two_prod(z, C2[1])
           t1,t2 = two_hilo_sum(u, -p1)
           v2 = (((t1-v1)+t2)-p2)
           return v2
       end

=#

#=
   function exp(x)
      signbit(x) && return inv(exp(abs(x))
=#

#=
 function v = exp( v )
            if ~isreal( v )
                [ sin_imag_v, cos_imag_v ] = sincos( imag( v ) );
                v = exp( real( v ) ) .* ( cos_imag_v + 1i .* sin_imag_v );
                return
            end

            % Strategy:  We first reduce the size of x by noting that
            % exp(kr + m * log(2)) = 2^m * exp(r)^k
            % where m and k are integers.  By choosing m appropriately
            % we can make |kr| <= log(2) / 2 = 0.347.  Then exp(r) is
            % evaluated using the familiar Taylor series.  Reducing the
            % argument substantially speeds up the convergence.
            k = 512.0;
            inv_k = 1.0 / k;
            Threshhold = inv_k .* DoubleDouble.eps.v1;

            m = floor( v.v1 ./ DoubleDouble.log_2.v1 + 0.5 );
            r = TimesPowerOf2( v - DoubleDouble.log_2 .* m, inv_k );

            p = r .* r;
            s = r + TimesPowerOf2( p, 0.5 );
            p = p .* r;
            t = p .* DoubleDouble.Make( DoubleDouble.InverseFactorial( 1, 1 ), DoubleDouble.InverseFactorial( 1, 2 ) );
            for i = 2 : DoubleDouble.NInverseFactorial
                s = s + t;
                p = p .* r;
                t = p .* DoubleDouble.Make( DoubleDouble.InverseFactorial( i, 1 ), DoubleDouble.InverseFactorial( i, 2 ) );
                if all( abs( t.v1(:) ) <= Threshhold )
                    break
                end
            end

            s = s + t;

            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = TimesPowerOf2( s, 2.0 ) + s .* s;
            s = s + 1.0;

            v = DoubleDouble.Make( pow2( s.v1, m ), pow2( s.v2, m ) );
        end

https://github.com/tholden/DoubleDouble/blob/master/DoubleDouble.m
=#
#=
https://www.pseudorandom.com/implementing-exp

def reduced_taylor_exp(x):
    """
    Evaluates f(x) = e^x for any x in the interval [-709, 709].
    If x < -709 or x > 709, raises an assertion error because of
    underflow/overflow limits. Performs a range reduction step
    by applying the identity e^x = e^r * 2^k for |r| <= 1/2 log(2)
    and positive integer k. e^r is approximated using a Taylor series
    truncated at 14 terms, which is enough because e^r is in the
    interval [0, 1/2 log(2)]. The result is left shifted by k
    to obtain e^r * 2^k = e^x.
    Performance: In the worst case we have 51 operations:
    - 16 multiplications
    - 16 divisions
    - 14 additions
    - 2 subtractions
    - 1 rounding
    - 1 left shift
    - 1 absolute value
    Accuracy: Over a sample of 10,000 equi-spaced points in
    [-709, 709] we have the following statistics:
    - Max relative error = 7.98411243625574e-14
    - Min relative error = 0.0
    - Avg relative error = 1.7165594254816366e-14
    - Med relative error = 1.140642685235478e-14
    - Var relative error = 2.964698541882666e-28
    - 6.29 percent of the values have less than 14 digits of precision
    :param x: (int) or (float) power of e to evaluate
    :return: (float) approximation of e^x
    """
    # If x is not a valid value, exit early
    assert(-709 <= x <= 709)
    # If x = 0, we know e^x = 1 so exit early
    if x == 0:
        return 1
    # Normalize x to a positive value since we can use reciprocal
    # symmetry to only work on positive values. Keep track of
    # the original sign for return value
    x0 = np.abs(x)
    # Hard code the value of natural log(2) in double precision
    l = 0.6931471805599453
    # Solve for an integer k satisfying x = k log(2) + r
    k = math.ceil((x0 / l) - 0.5)
    # p = 2^k
    p = 1 << k
    # r is a value between 0 and 0.5 log(2), inclusive
    r = x0 - (k * l)
    # Setup the Taylor series to evaluate e^r after
    # range reduction x -> r. The Taylor series
    # only approximates on the interval [0, log(2)/2]
    Tn = 1
    # We need at most 14 terms to achieve 16 digits
    # of precision anywhere on the interval [0, log(2)/2]
    for i in range(14, 0, -1):
        Tn = Tn * (r / i) + 1
    # e^x = e^r * 2^k, so multiply p by Tn. This loses us
    # about two digits of precision because we compound the
    # relative error by 2^k.
    p *= Tn
    # If the original sign is negative, return reciprocal
    if x < 0:
        p = 1 / p
    return p

double reduced_taylor_exp(double x) {
    /*
     * Evaluates f(x) = e^x for any x in the interval [-709, 709].
     * If x < -709 or x > 709, raises an assertion error. Performs a
     * range reduction step by applying the identity e^x = e^r * 2^k
     * for |r| <= 1/2 log(2) and positive integer k. e^r is evaluated
     * using Taylor series truncated at 14 terms, which is sufficient
     * to achieve 16 digits of precision for r so close to 0. The
     * result is left shifted by k to obtain e^r * 2^k = e^x.
     * Performance: In the worst case we have 51 operations:
     * - 16 multiplications
     * - 16 divisions
     * - 14 additions
     * - 2 subtractions
     * - 1 rounding
     * - 1 left shift
     * - 1 absolute value
     * Accuracy: Over a sample of 10,000 linearly spaced points in
     * [-709, 709], we have the following error statistics:
     * - Max relative error = 7.99528e-14
     * - Min relative error = 0
     * - Avg relative error = 0
     * - Med relative error = 2.27878e-14
     * - Var relative error = 0
     * - 6.4 percent of the values have less than 14 digits of precision.
     * Args:
     *      - x (double float) power of e to evaluate
     * Returns:
     *      - (double float) approximation of e^x
     */
    // Make sure x is a valid input
    assert(-709 <= x && x <= 709);
    // When x is 0, we know e^x is 1
    if (x == 0) {
        return 1;
    }
    // Normalize x to a positive value to take advantage of
    // reciprocal symmetry, but keep track of the original value
    double x0 = abs(x);
    // Solve for k satisfying x = 2^k * log(2) r, with |r| <= 1/2 log(2)
    int k = ceil((x0 / M_LN2) - 0.5);
    // Determine r using the k we computed
    double r = x0 - (k * M_LN2);
    // Setup the Taylor series to evaluate e^r after range reduction
    // x -> r. This only approximates over the interval [0, 1/2 log(2)]
    double tk = 1.0;
    double tn = 1.0;
    // We only need 14 terms to achieve 16 digits of precision for e^r
    for (int i = 1; i < 14; i++) {
        tk *= r / i;
        tn += tk;
    };
    tn = tn * pow2(k);
    if (x < 0) {
        return 1 / tn;
    }
    return tn;
}

=#

#=
kv-master

```

julia> x
Double64((19.73721077166514, -1.6936500472509631e-15))

julia> x/Ln2
Double64((28.474776101261526, -8.19680688222626e-17))

julia> xoLn2 = ans
Double64((28.474776101261526, -8.19680688222626e-17))

julia> x_i = floor(x)
Double64((19.0, 0.0))

julia> x_f = x - x_i
Double64((0.7372107716651383, -2.831551031322832e-17))

julia> x_f = x_f - 1
Double64((-0.2627892283348618, 2.7195640918029507e-17))

julia> x_i = x_i + 1
Double64((20.0, 0.0))

julia> r = y = Double64(1);

julia> for i=1:25
         y = y * x_f
         y = y / i
                  r = r + y
       end

julia> x_i
Double64((20.0, 0.0))

julia> e1 = Double64(exp(BigFloat(1)))
Double64((2.718281828459045, 1.4456468917292502e-16))

julia> e2 = Double64(BigFloat(e1)^BigFloat(x_i))
Double64((4.851651954097903e8, 4.880277289790481e-10))

julia> r2 = r * e2
Double64((3.730454319498937e8, -6.494042474083988e-10))

julia> expa = Double64(exp(BigFloat(a)))
Double64((3.730454319498937e8, -6.49404247408404e-10))
```
	friend dd exp(const dd& x) {
		dd x_i, x_f, tmp;
		dd r, y;
		int i;

		if (x == dd(std::numeric_limits<double>::infinity(), 0.)) {
			return dd(std::numeric_limits<double>::infinity(), 0.);
		}
		if (x == -dd(std::numeric_limits<double>::infinity(), 0.)) {
			return (dd)0.;
		}

		if (x >= 0.) {
			x_i = floor(x);
			x_f = x - x_i;
			if (x_f >= 0.5) {
				x_f -= 1.;
				x_i += 1.;
			}
		} else {
			x_i = -floor(-x);
			x_f = x - x_i;
			if (x_f <= -0.5) {
				x_f += 1.;
				x_i -= 1.;
			}
		}

		r = 1.;
		y = 1.;
		for (i=1;  i<=25 ; i++) {
			y *= x_f;
			y /= i;
			r += y;
		}

		if (x_i >= 0.) {
			// r *= pow(constants<dd>::e(), (int)x_i);
			r *= ipower(e(), (double)x_i);
		} else {
			// r /= pow(constants<dd>::e(), -(int)x_i);
			r /= ipower(e(), -(double)x_i);
		}

		return r;
	}

	static dd expm1_origin(const dd& x) {
		dd r, y;
		int i;

		r = 0.;
		y = 1.;
		for (i=1; i<=25 ; i++) {
			y *= x;
			y /= i;
			r += y;
		}

		return r;
	}

	friend dd expm1(const dd& x) {
		if (x >= -0.5 && x <= 0.5) {
			return expm1_origin(x);
		} else {
			return exp(x) - 1.;
		}
	}

=#

#=
   function expsmall6(v)
    t1 = v ^ 2
    t5 = t1 ^ 2
    t13 = t5 ^ 2
    return(1 + v + (3 + v) * t1 / 6 + (5 + v) * t5 / 120 + (7 + v) * t1 * t5 / 5040 + (9 + v) * t13 / 362880 + (11 + v) * t1 * t13 / 39916800 + (13 + v) * t5 * t13 / 6227020800)
   end

   function expsmall8horner(v)
      1 + (1 + (1//2 + (1//6 + (1//24 + (1//120 + (1//720 + (1//5040 + (1//40320 + (1//362880 + (1//3628800 + (1//39916800 + (1//479001600 + (1//6227020800 +
          (1//87178291200 + v / 1307674368000) * v) * v) * v) * v) * v) * v) * v) * v) * v) * v) * v) * v) * v) * v
    end
=#


text:
function GetAttribLocation(prog::Program, name::String)
	ret = ccall( (:glGetAttribLocation, lib), Attribute,
		(Program, Ptr{GLchar}), prog, bytestring(name))
	if ret < 0
		GetError()
		error("attribute not found")
	else
		return ret
	end
end

function EnableVertexAttribArray(attr::Attribute)
	ccall( (:glEnableVertexAttribArray, lib), Void, (Attribute,), attr)
end

function VertexAttribPointer(attr::Attribute, size::Integer, type_::Integer,
	normalize::Bool, stride::Integer, ptr::Integer)

	ccall( (:glVertexAttribPointer, lib), Void,
		(Attribute, GLint, GLenum, GLboolean, GLsizei, Ptr{GLvoid}),
		attr, size, type_, normalize, stride, ptr)
	GetError()
end


text:
using StenoGraphs
using Test, SafeTestsets, Documenter

@testset "StenoGraphs.jl" begin
    DocMeta.setdocmeta!(StenoGraphs, :DocTestSetup, :(using StenoGraphs); recursive=true)
    doctest(StenoGraphs)

    @safetestset "Graph Types" begin
        using StenoGraphs
        include("graph_types.jl")
    end
    @safetestset "Syntax" begin
        using StenoGraphs
        include("syntax.jl")
    end
    @safetestset "Modify Edges" begin
        using StenoGraphs
        include("modification.jl")
    end
    @safetestset "Communative" begin
        using StenoGraphs
       include("comunative.jl")
    end
    @safetestset "Promotion and Conversion" begin
        using StenoGraphs
        include("conversion.jl")
    end
    @safetestset "Comparison" begin
        using StenoGraphs
        include("compare.jl")
    end
    @safetestset "Unmeta" begin
        using StenoGraphs
        include("keep.jl")
    end
    @safetestset "Arrow" begin
        using StenoGraphs
        include("arrow.jl")
    end
    @safetestset "Deprecate" begin
        using StenoGraphs
        include("deprecate.jl")
    end
    @safetestset "Ids" begin
        using StenoGraphs
        include("id.jl")
    end
    @safetestset "Meld" begin
        using StenoGraphs
        include("meld.jl")
    end
    @safetestset "Swap" begin
        using StenoGraphs
        include("swap.jl")
    end
end


text:
# source: https://gist.github.com/denizyuret/bf74f737ba8a7076a7f43c5080555a9a

using AutoGrad, CuArrays, Knet, Test
using AutoGrad: @gcheck
using CuArrays: @cufunc

const GConstant01 = sqrt(2/pi)
const GConstant02 = 0.044715 * sqrt(2/pi)
const GConstant03 = GConstant01 / 2

# Main definition, broadcasted version works on Arrays

gelu(x::T) where T = (x/2)*(1 + tanh(T(GConstant02)*x^3 + T(GConstant01)*x))
geluback(x::T,dy::T) where T = dy*(T(0.5)*tanh(T(GConstant02)*x^3 + T(GConstant01)*x) + (T(0.0535161)*x^3 + T(GConstant03)*x)*(1/cosh(T(GConstant02)*x^3 + T(GConstant01)*x))^2 + T(0.5))

# This defines gelu for AutoGrad

@primitive  gelu(x),dy  geluback.(x,dy)

# This makes broadcasted gelu work for CuArray, directly compiling a CUDA kernel from Julia code.
# @cufunc necessary to use a GPU compatible version of tanh etc.
# https://github.com/JuliaGPU/CuArrays.jl/issues/253

@cufunc gelu(x::T) where T = (x/2)*(1 + tanh(T(GConstant02)*x^3 + T(GConstant01)*x))
@cufunc geluback(x::T,dy::T) where T = dy*(T(0.5)*tanh(T(GConstant02)*x^3 + T(GConstant01)*x) + (T(0.0535161)*x^3 + T(GConstant03)*x)*(1/cosh(T(GConstant02)*x^3 + T(GConstant01)*x))^2 + T(0.5))

# This defines gelu for KnetArray
import Base.Broadcast: broadcasted
import Knet: KnetArray

function KnetArray(x::CuArray{T,N}) where {T,N}
    p = Base.bitcast(Knet.Cptr, x.ptr)
    k = Knet.KnetPtr(p, sizeof(x), gpu(), x)
    KnetArray{T,N}(k, size(x))
end

broadcasted(::typeof(gelu),x::KnetArray) = KnetArray(gelu.(CuArray(x)))
broadcasted(::typeof(geluback),x::KnetArray,dy::KnetArray) = KnetArray(geluback.(CuArray(x),CuArray(dy)))

# This tests all of the above
# julia -L gelu.jl -e 'testgelu()'
function testgelu()
    @testset "gelu" begin
        x = rand(10)
        a = Param(x)
        c = Param(CuArray(x))
        k = Param(KnetArray(x))
        @test @gcheck gelu.(a)
        @test @gcheck gelu.(c)
        @test @gcheck gelu.(k)
    end
end

"""
    gelu(x)

    Gaussian Error Linear Unit Activation Function
    https://arxiv.org/abs/1606.08415

    returns:
        0.5x * (1 + tanh(√(2/π) * (x + 0.044715x^3)))
"""
gelu


text:
#=

Create sysimage using

    > using PackageCompiler
    > create_sysimage([:Images, :ImageView, :Random, :SparseArrays, :LinearAlgebra, :Statistics], sysimage_path="sys_img.so")

and start julia using
    $ julia -Jsys_img.so main.jl PathToNormImage PathToActualImage =#

using ImageView, Images
using SparseArrays
using Random
using LinearAlgebra
using Statistics


text:
include("utils/misc.jl")
include("utils/paths.jl")
include("utils/fingerprint.jl")
include("utils/parse_includes.jl")
include("utils/poll.jl")
include("utils/hot_restart.jl")


text:
## ----------------------------------------------------------------------------------
function collect_lp_DAT!(DAT, Ed, iJR)

    proj = string(parentmodule(iJR))
    src = string(nameof(Ed))

    rxns_map = iJR.load_rxns_map()

    objider = iJR.BIOMASS_IDER
    costider = iJR.COST_IDER

    LP_DAT = ChE.load_LP_DAT(src)

    for method in DAT[:LP_METHODS]

        for exp in DAT[:EXPS]

            @info("Doing",
                proj, exp, method
            ); println()

            model = LP_DAT[method, :model, exp]
            fbaout = LP_DAT[method, :fbaout, exp]

            # Biomass
            fba_biom = ChU.av(model, fbaout, objider)
            DAT[method, :flx, "D", exp] = fba_biom
            DAT[method, :err, "D", exp] = 0.0

            fba_cost = ChU.av(model, fbaout, costider)
            DAT[method, :flx, "cost", exp] = fba_cost
            DAT[method, :err, "cost", exp] = 0.0

            for Ed_ider in DAT[:EXCH_FLX_IDERS]
                model_ider = rxns_map[Ed_ider]

                fba_flx = ChU.av(model, fbaout, model_ider)
                DAT[method, :flx, Ed_ider, exp] = fba_flx
                DAT[method, :err, Ed_ider, exp] = 0.0
            end

            idermap = merge(iJR.load_inners_idermap(), iJR.load_inners_idermap())
            for (exglcider, model_iders) in idermap
                # flxs
                fba_flx = ChU.av(model, fbaout, model_iders[1])
                if length(model_iders) == 2 # reversible
                    # r = r+ - r-
                    fba_flx -= ChU.av(model, fbaout, model_iders[2])
                end

                DAT[method, :flx, exglcider, exp] = fba_flx
                DAT[method, :err, exglcider, exp] = 0.0
            end

        end

    end # for method

    return DAT
end


text:
groups_text = split(read("day6/input.txt", String), "\n\n")

groups = split.(groups_text, "\n")

value = 0
for group in groups
    group_questions = Set{Char}()
    for person in group
        push!(group_questions, collect(person)...)
    end
    global value += length(group_questions)
end

println(value)


text:
__precompile__(true)

"""
"""
module PredictMD # begin module PredictMD

import Distributed
import Random

include("base/api.jl")

include("base/backends.jl")
include("base/concrete-types.jl")
include("base/fallback.jl")

include("registry_url_list.jl")
include("package_directory.jl")
include("version.jl")
include("package_list.jl")
include("import_all.jl")
include("welcome.jl")
include("init.jl")

include("toplevel/always-loaded/code_loading/requires.jl")
include("toplevel/always-loaded/code_loading/require_versions.jl")

include("toplevel/always-loaded/classimbalance/smote.jl")

include("toplevel/always-loaded/datasets/csv.jl")
include("toplevel/always-loaded/datasets/datadeps.jl")
include("toplevel/always-loaded/datasets/gzip.jl")
include("toplevel/always-loaded/datasets/juliadb.jl")
include("toplevel/always-loaded/datasets/mldatasets.jl")
include("toplevel/always-loaded/datasets/mnist.jl")
include("toplevel/always-loaded/datasets/queryverse.jl")
include("toplevel/always-loaded/datasets/rdatasets.jl")

include("toplevel/always-loaded/docs_and_examples/cache.jl")
include("toplevel/always-loaded/docs_and_examples/generate_examples.jl")

include("toplevel/always-loaded/ide/atom.jl")

include("toplevel/always-loaded/io/saveload.jl")

include("toplevel/always-loaded/linearmodel/glm.jl")
include("toplevel/always-loaded/linearmodel/ordinary_least_squares_regression.jl")

include("toplevel/always-loaded/metrics/auprc.jl")
include("toplevel/always-loaded/metrics/aurocc.jl")
include("toplevel/always-loaded/metrics/averageprecisionscore.jl")
include("toplevel/always-loaded/metrics/brier_score.jl")
include("toplevel/always-loaded/metrics/coefficientofdetermination.jl")
include("toplevel/always-loaded/metrics/cohenkappa.jl")
include("toplevel/always-loaded/metrics/getbinarythresholds.jl")
include("toplevel/always-loaded/metrics/mean_square_error.jl")
include("toplevel/always-loaded/metrics/prcurve.jl")
include("toplevel/always-loaded/metrics/risk_score_cutoff_values.jl")
include("toplevel/always-loaded/metrics/roccurve.jl")
include("toplevel/always-loaded/metrics/rocnumsmetrics.jl")
include("toplevel/always-loaded/metrics/singlelabelbinaryclassificationmetrics.jl")
include("toplevel/always-loaded/metrics/singlelabelregressionmetrics.jl")

include("toplevel/always-loaded/modelselection/crossvalidation.jl")
include("toplevel/always-loaded/modelselection/split_data.jl")

include("toplevel/always-loaded/neuralnetwork/flux.jl")
include("toplevel/always-loaded/neuralnetwork/knet.jl")
include("toplevel/always-loaded/neuralnetwork/merlin.jl")

include("toplevel/always-loaded/online/onlinestats.jl")

include("toplevel/always-loaded/ontologies/ccs.jl")

include("toplevel/always-loaded/pipeline/simplelinearpipeline.jl")

include("toplevel/always-loaded/plotting/catch_plotting_errors.jl")
include("toplevel/always-loaded/plotting/defaultapplication.jl")
include("toplevel/always-loaded/plotting/pgfplots.jl")
include("toplevel/always-loaded/plotting/pgfplotsx.jl")
include("toplevel/always-loaded/plotting/plotlearningcurve.jl")
include("toplevel/always-loaded/plotting/plotprcurve.jl")
include("toplevel/always-loaded/plotting/plotroccurve.jl")
include("toplevel/always-loaded/plotting/plotsinglelabelregressiontruevspredicted.jl")
include("toplevel/always-loaded/plotting/plotsinglelabelbinaryclassifierhistograms.jl")
include("toplevel/always-loaded/plotting/probability_calibration_plots.jl")
include("toplevel/always-loaded/plotting/unicodeplots.jl")

include("toplevel/always-loaded/postprocessing/packagemultilabelpred.jl")
include("toplevel/always-loaded/postprocessing/packagesinglelabelpred.jl")
include("toplevel/always-loaded/postprocessing/packagesinglelabelproba.jl")
include("toplevel/always-loaded/postprocessing/predictoutput.jl")
include("toplevel/always-loaded/postprocessing/predictprobaoutput.jl")

include("toplevel/always-loaded/preprocessing/dataframecontrasts.jl")
include("toplevel/always-loaded/preprocessing/dataframetodecisiontree.jl")
include("toplevel/always-loaded/preprocessing/dataframetoglm.jl")
include("toplevel/always-loaded/preprocessing/dataframetoknet.jl")
include("toplevel/always-loaded/preprocessing/dataframetosvm.jl")

include("toplevel/always-loaded/svm/libsvm.jl")

include("toplevel/always-loaded/time_series/timeseries.jl")

include("toplevel/always-loaded/tree/decisiontree.jl")

include("toplevel/always-loaded/utils/constant_columns.jl")
include("toplevel/always-loaded/utils/dataframe_column_types.jl")
include("toplevel/always-loaded/utils/filename_extension.jl")
include("toplevel/always-loaded/utils/find.jl")
include("toplevel/always-loaded/utils/fix_type.jl")
include("toplevel/always-loaded/utils/formulas.jl")
include("toplevel/always-loaded/utils/inverse-dictionary.jl")
include("toplevel/always-loaded/utils/is_debug.jl")
include("toplevel/always-loaded/utils/labelstringintmaps.jl")
include("toplevel/always-loaded/utils/linearly_dependent_columns.jl")
include("toplevel/always-loaded/utils/make_directory.jl")
include("toplevel/always-loaded/utils/maketemp.jl")
include("toplevel/always-loaded/utils/missings.jl")
include("toplevel/always-loaded/utils/nothings.jl")
include("toplevel/always-loaded/utils/openbrowserwindow.jl")
include("toplevel/always-loaded/utils/openplotsduringtestsenv.jl")
include("toplevel/always-loaded/utils/predictionsassoctodataframe.jl")
include("toplevel/always-loaded/utils/probabilitiestopredictions.jl")
include("toplevel/always-loaded/utils/runtestsenv.jl")
include("toplevel/always-loaded/utils/shufflerows.jl")
include("toplevel/always-loaded/utils/simplemovingaverage.jl")
include("toplevel/always-loaded/utils/tikzpictures.jl")
include("toplevel/always-loaded/utils/transform_columns.jl")
include("toplevel/always-loaded/utils/trapz.jl")
include("toplevel/always-loaded/utils/traviscienv.jl")
include("toplevel/always-loaded/utils/tuplify.jl")

include("submodules/Cleaning/Cleaning.jl")
include("submodules/Compilation/Compilation.jl")
include("submodules/GPU/GPU.jl")
include("submodules/Server/Server.jl")

end # end module PredictMD


text:
# Code from Hans W. Borschers
function modLangerman(x)
    #if length(x) != 10; error("ModLangerman needs a 10-elements vector"); end
    sqr(z) = z^2
    x1 = x[1]; x2 = x[2]; x3 = x[3]; x4 = x[4]; x5 = x[5]
    x6 = x[6]; x7 = x[7]; x8 = x[8]; x9 = x[9]; x10 = x[10]

    y = 0.806*exp(-0.31830988618377*(sqr((-9.681) + x1) + sqr((-0.667) + x2) +
     sqr((-4.783) + x3) + sqr((-9.095) + x4) + sqr((-3.517) + x5) + sqr((-9.325
     ) + x6) + sqr((-6.544) + x7) + sqr((-0.211) + x8) + sqr((-5.122) + x9) +
     sqr((-2.02) + x10)))*cos(3.14159265359*(sqr((-9.681) + x1) + sqr((-0.667)
      + x2) + sqr((-4.783) + x3) + sqr((-9.095) + x4) + sqr((-3.517) + x5) +
     sqr((-9.325) + x6) + sqr((-6.544) + x7) + sqr((-0.211) + x8) + sqr((-5.122
     ) + x9) + sqr((-2.02) + x10))) + 0.517*exp(-0.31830988618377*(sqr((-9.4)
      + x1) + sqr((-2.041) + x2) + sqr((-3.788) + x3) + sqr((-7.931) + x4) +
     sqr((-2.882) + x5) + sqr((-2.672) + x6) + sqr((-3.568) + x7) + sqr((-1.284
     ) + x8) + sqr((-7.033) + x9) + sqr((-7.374) + x10)))*cos(3.14159265359*(
     sqr((-9.4) + x1) + sqr((-2.041) + x2) + sqr((-3.788) + x3) + sqr((-7.931)
      + x4) + sqr((-2.882) + x5) + sqr((-2.672) + x6) + sqr((-3.568) + x7) +
     sqr((-1.284) + x8) + sqr((-7.033) + x9) + sqr((-7.374) + x10))) + 0.1*exp(
     -0.31830988618377*(sqr((-8.025) + x1) + sqr((-9.152) + x2) + sqr((-5.114)
      + x3) + sqr((-7.621) + x4) + sqr((-4.564) + x5) + sqr((-4.711) + x6) +
     sqr((-2.996) + x7) + sqr((-6.126) + x8) + sqr((-0.734) + x9) + sqr((-4.982
     ) + x10)))*cos(3.14159265359*(sqr((-8.025) + x1) + sqr((-9.152) + x2) +
     sqr((-5.114) + x3) + sqr((-7.621) + x4) + sqr((-4.564) + x5) + sqr((-4.711
     ) + x6) + sqr((-2.996) + x7) + sqr((-6.126) + x8) + sqr((-0.734) + x9) +
     sqr((-4.982) + x10))) + 0.908*exp(-0.31830988618377*(sqr((-2.196) + x1) +
     sqr((-0.415) + x2) + sqr((-5.649) + x3) + sqr((-6.979) + x4) + sqr((-9.51)
      + x5) + sqr((-9.166) + x6) + sqr((-6.304) + x7) + sqr((-6.054) + x8) +
     sqr((-9.377) + x9) + sqr((-1.426) + x10)))*cos(3.14159265359*(sqr((-2.196)
      + x1) + sqr((-0.415) + x2) + sqr((-5.649) + x3) + sqr((-6.979) + x4) +
     sqr((-9.51) + x5) + sqr((-9.166) + x6) + sqr((-6.304) + x7) + sqr((-6.054)
      + x8) + sqr((-9.377) + x9) + sqr((-1.426) + x10))) + 0.965*exp(-
     0.31830988618377*(sqr((-8.074) + x1) + sqr((-8.777) + x2) + sqr((-3.467)
      + x3) + sqr((-1.867) + x4) + sqr((-6.708) + x5) + sqr((-6.349) + x6) +
     sqr((-4.534) + x7) + sqr((-0.276) + x8) + sqr((-7.633) + x9) + sqr((-1.567
     ) + x10)))*cos(3.14159265359*(sqr((-8.074) + x1) + sqr((-8.777) + x2) +
     sqr((-3.467) + x3) + sqr((-1.867) + x4) + sqr((-6.708) + x5) + sqr((-6.349
     ) + x6) + sqr((-4.534) + x7) + sqr((-0.276) + x8) + sqr((-7.633) + x9) +
     sqr((-1.567) + x10)))

    return -y
end

using BlackBoxOptim

p = BlackBoxOptim.fixeddim_problem(modLangerman; range = (0.0, 10.0), dims = 10,
  name = "Modified Langerman", fmins = [-0.9650])

BlackBoxOptim.repeated_bboptimize(30, p, 10, [
  :generating_set_search,
  :adaptive_de_rand_1_bin_radiuslimited],
  10.0, 1e-5, {:MinDeltaFitnessTolerance => 1e-50})


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

using MLJ, RDatasets
MLJ.color_off() # hide
auto = dataset("ISLR", "Auto")
y, X = unpack(auto, ==(:MPG), col->true)
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=444);

@load LinearRegressor pkg=MLJLinearModels

using PyPlot

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-5-g1.svg")) # hide

lm = LinearRegressor()
mlm = machine(lm, select(X, :Horsepower), y)
fit!(mlm, rows=train)
rms(predict(mlm, rows=test), y[test])^2

xx = (Horsepower=range(50, 225, length=100) |> collect, )
yy = predict(mlm, xx)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy, lw=3)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-5-g2.svg")) # hide

hp = X.Horsepower
Xhp = DataFrame(hp1=hp, hp2=hp.^2, hp3=hp.^3);

@pipeline LinMod(fs = FeatureSelector(features=[:hp1]),
                 lr = LinearRegressor());

lrm = LinMod()
lr1 = machine(lrm, Xhp, y) # poly of degree 1 (line)
fit!(lr1, rows=train)

lrm.fs.features = [:hp1, :hp2] # poly of degree 2
lr2 = machine(lrm, Xhp, y)
fit!(lr2, rows=train)

lrm.fs.features = [:hp1, :hp2, :hp3] # poly of degree 3
lr3 = machine(lrm, Xhp, y)
fit!(lr3, rows=train)

get_mse(lr) = rms(predict(lr, rows=test), y[test])^2

@show get_mse(lr1)
@show get_mse(lr2)
@show get_mse(lr3)

hpn  = xx.Horsepower
Xnew = DataFrame(hp1=hpn, hp2=hpn.^2, hp3=hpn.^3)

yy1 = predict(lr1, Xnew)
yy2 = predict(lr2, Xnew)
yy3 = predict(lr3, Xnew)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy1, lw=3, label="Order 1")
plot(xx.Horsepower, yy2, lw=3, label="Order 2")
plot(xx.Horsepower, yy3, lw=3, label="Order 3")

legend(fontsize=14)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-5-g3.svg")) # hide

Xhp = DataFrame([hp.^i for i in 1:10])

cases = [[Symbol("x$j") for j in 1:i] for i in 1:10]
r = range(lrm, :(fs.features), values=cases)

tm = TunedModel(model=lrm, ranges=r, resampling=CV(nfolds=10), measure=rms)

mtm = machine(tm, Xhp, y)
fit!(mtm)
rep = report(mtm)

res = rep.plotting

@show round.(res.measurements.^2, digits=2)
@show argmin(res.measurements)

Xnew = DataFrame([hpn.^i for i in 1:10])
yy5 = predict(mtm, Xnew)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy5, lw=3)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

savefig(joinpath(@OUTPUT, "ISL-lab-5-g4.svg")) # hide

PyPlot.close_figs() # hide

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
module TestBlockArrays

include("preamble.jl")
using BlockArrays

enumndim(size...) = collect(reshape(1:prod(size), size))

@testset "BlockVector" begin
    B = BlockArray(1:4, [1, 3])
    A = collect(B)
    @test collect(Map(identity), B) == A
end

@testset "BlockMatrix" begin
    B = BlockArray(enumndim(2, 3), [1, 1], [2, 1])
    A = collect(B)
    @test collect(Map(identity), B) == A[:]
end

@testset "BlockArray{T, 3}" begin
    B = BlockArray(enumndim(2, 3, 4), [1, 1], [2, 1], [1, 1, 2])
    A = collect(B)
    @test collect(Map(identity), B) == A[:]
end

end  # module


text:
#P6
#1920 1080
#255
#<frame RGB>
#
# ASCII dump:
# 0x50 0x36   0x0a
# 0x31 0x39 0x32 0x30   0x20   0x31 0x30 0x38 0x30   0x0a
# 0x32 0x35 0x35   0x0a
#
# lixun 2017/08/09
p6 = [0x50, 0x36, 0x0a]
horizontal = "480"
vertical = "640"
fps = 30
duration = 3600.7
start = 1
depth = [0x32, 0x35, 0x35, 0x0a]

noise_floor = "alexa-f0.ppm"
digest = "alexa.digest"

function delta(a, b)
    a > b && (return a-b)
    b-a
end

h = tryparse(Int64, horizontal)
h = isnull(h) ? 0 : get(h)
v = tryparse(Int64, vertical)
v = isnull(v) ? 0 : get(v)

header = zeros(UInt8, length(p6) + length(horizontal) + 1 + length(vertical) + 1 + length(depth))
f0 = zeros(UInt8, h*v*3)
fx = zeros(UInt8, h*v*3)
rgb = zeros(Float64, convert(Int64, ceil(duration * fps)), 3)

# capture the inactive frame as basis
fn = convert(Int64, ceil(start * fps))
for i = 1:fn
    read!(stdin, header)
    read!(stdin, f0)
end
open(noise_floor, "w") do f
    write(f, header)
    write(f, f0)
end

# dealing with the rest frames
i = 1
while !eof(stdin)
    read!(stdin, header)
    read!(stdin, fx)
    #version 1: signal[i] = sum(frame) - md
    for k = 1:3
        rgb[i,k] = sum(delta.(fx[k:3:end], f0[k:3:end])) / (h*v)
    end
    #version 3: signal[i] = delta(sum(frame), md)
    i += 1
end
open(digest, "w") do f
    for j = 1:i
        write(f, "$(rgb[j,1]), $(rgb[j,2]), $(rgb[j,3])\n")
    end
end


text:
using Continuation
using Test

@testset "Continuation.jl" begin
    # Write your tests here.
end


text:
using OpenCL, CLBLAS, CLBlast, BenchmarkTools
@static if VERSION < v"0.7-"
    LA = Base.LinAlg
else
    using LinearAlgebra, Random, Printf
    LA = LinearAlgebra
end

CLBLAS.setup()

function run_scal(n::Integer, T=Float32)
    srand(12345)
    x = rand(T, n)
    α = rand(T)

    @printf("\nn = %d, eltype = %s\n", n, T)

    println("BLAS:")
    x_true = copy(x)
    LinAlg.BLAS.scal!(n, α, x_true, 1)
    xx = copy(x)
    display(@benchmark LinAlg.BLAS.scal!($n, $α, $xx, $1))
    println()

    for device in cl.devices()
        println("-"^70)
        @printf("Platform name   : %s\n", device[:platform][:name])
        @printf("Platform version: %s\n", device[:platform][:version])
        @printf("Device name     : %s\n", device[:name])
        @printf("Device type     : %s\n", device[:device_type])
        println()

        ctx = cl.Context(device)
        queue = cl.CmdQueue(ctx)

        println("CLBLAS:")
        x_cl = cl.CLArray(queue, x)
        CLBLAS.scal!(n, α, x_cl, 1)
        @assert cl.to_host(x_cl) ≈ x_true
        display(@benchmark CLBLAS.scal!($n, $α, $x_cl, $1))
        println()

        println("CLBlast:")
        x_cl = cl.CLArray(queue, x)
        CLBlast.scal!(n, α, x_cl, 1)
        # or CLBlast.scal!(α, x_cl)
        @assert cl.to_host(x_cl) ≈ x_true
        display(@benchmark CLBlast.scal!($n, $α, $x_cl, $1))
        println()
    end
end

run_scal(2^10, Float32)