text:
using VegaLite
using DataFrames
using LinearAlgebra
using Combinatorics

# These are the same plots as in the web app, data is added after the fact, so stubbed values are used at first in the specs
# syntax highlighting not so great...

function barplot(stats::Stats)
  spec = vl"""
  {
    "title": {
      "text": "Top $(length(stats.mesh_names)) Co-Occuring Terms by Count"
    },
    "data": {
      "values": []
    },
    "mark": "bar",
    "encoding": {
      "x": {
        "field": "x",
        "type": "nominal",
        "sort": { "field": "y" }
      },
      "y": {
        "field": "y",
        "type": "quantitative"
      },
      "color": {
        "value": "\u00233B5B8C"
      },
      "tooltip" : [
        {
          "field": "x",
          "type": "nominal",
          "title": "MeSH"
        },
        {
          "field": "y",
          "type": "quantitative",
          "title": "Count"
        }
      ]
    }
  }
  """
  df = DataFrame(x=stats.mesh_names, y=stats.mesh_counts)
  spec(df)

  return nothing
end

function stats_matrix(stats::Stats, measure::Symbol)

  maxExtent = maximum(abs.(vals))
  minExtent = -1*maxExtent

  data = []
  for i = 1:length(stats.mesh_names), j = 1:length(stats.mesh_names)
    if i != j
      push!(data, (x=i,y=j,val=stats.measure[i,j]))
    end
  end

  table = data |> columntable

  spec = vl"""
  {
    "title": {
      "text": $string(measure)
    },
    "width": 700,
    "height": 700,
    "data": {
      "values": []
    },
    "transform": [
      {
        "filter": "datum.x != datum.y"
      }
    ],
    "config": { "axis": { "domain": false }},
    "mark": "rect",
    "encoding": {
      "x": {
        "field": "x",
        "type": "nominal",
        "axis": {
          "orient": "top",
          "ticks": false,
          "domain": false,
          "labelPadding": 4,
          "minExtent": 180
        }
       },
      "y": {
        "field": "y",
        "type": "nominal",
        "axis": {
          "ticks": false,
          "domain": false,
          "labelPadding": 4,
          "minExtent": 180
        }
      },
      "color": {
        "field": "val",
        "type": "quantitative",
        "legend": {
            "title": null
        },
        "scale": {
          "domain": [$minExtent, 0, $maxExtent],
          "range": ["\u00234d4d4d", "\u0023bababa", "\u0023ffffff", "\u0023ffffff", "\u002367a9cf","\u00233B5B8C"]
        }
      },
      "tooltip": [
        {
          "field": "x",
          "type": "nominal",
          "title": "MeSH 1"
        },
        {
          "field": "y",
          "type": "nominal",
          "title": "MeSH 2"
        },
        {
          "field": "val",
          "type": "quantitative",
          "title": "Value"
        }
      ]
    }
  """
  spec(table)
end

function chord_diagram(stats::Stats)

  chord_rules = []
  ruleID = 0
  for row in eachrow(stats.mh_rules_df)
    lhs = row[:lhs][2:end-1]
    lhsArr = split(lhs, " | ")

    rhs = row[:rhs]
    shortName = length(rhs) > 25 ? rhs[1:25] * "..." : rhs

    perms = permutations(lhsArr) |> collect

    for perm in perms
      newLhs = join(perm, " | ")
      fset = rhs * " | " * newLhs
      ruleID += 1
      push!(chordRules, (id=ruleID, lhs=newLhs, rhs=rhs, sname=shortName, fset=fset, fsetSize=length(perm), conf=row[:conf], lift=row[:lift], supp=row[:supp]))
    end
  end

  sort!(chordRules, by= x -> x.supp, rev=true)

  chordIDs = []

  for i in 1:length(stats.topn_mesh_labels) # topn_codes
    shortName = length(stats.topn_mesh_labels[i]) > 25 ? stats.topn_mesh_labels[i][1:25] * "..." : stats.topn_mesh_labels[i]

    push!(chordIDs, (id=i, name=stats.topn_mesh_labels[i], sname=shortName, supp=(stats.topn_mesh_counts[i] / stat.mh_rules_pmid_count)))

    ruleID += 1
    push!(chordRules, (id=ruleID, lhs=missing, rhs=stats.topn_mesh_labels[i], sname=shortName, fset=stats.topn_mesh_labels[i], fsetSize=1, conf=0, lift=0, supp=(stats.topn_mesh_counts[i] / stat.mh_rules_pmid_count)))
  end

  chordRels=[]
  for i = 1:length(stats.topn_mesh_labels), j = 1:length(stats.topn_mesh_labels)
    if i < j && stats.top_coo_sp[i][j] > 0
      push!(chordRels, (source=i, target=j, size=stats.top_coo_sp[i][j]))
  end

  ruleStr = json(chordRules)
  idStr = json(chordIDs)
  relStr = json(chordRels)

  dataSize = length(chordIDs)
  maxEdgeSize = maximum(x->x.size, chordIDs)
  maxFreqItems = maximum(x->x.fsetSize, chordRules)

  spec = """
  {
    "\$schema": "https://vega.github.io/schema/vega/v4.json",
    "autosize": "pad",
    "padding": 5,
    "height": 700,

    "signals": [
      { "name": "chordWidth", "value": 600},
      { "name": "originX", "update": "chordWidth / 2" },
      { "name": "originY", "update": "height / 2" },
      { "name": "radius", "update": "chordWidth / 2 - 175" },
      { "name": "dataSize", "update": $dataSize },
      { "name": "maxEdgeSize", "update": $maxEdgeSize },
      { "name": "maxFreqItems", "value": $maxFreqItems - 1 },
      { "name": "textOffset", "value": 5 },
      { "name": "textSize", "value": 11 },
      { "name": "treeWidth", "value": ($maxFreqItems-1)*175 + 320},
      {
        "name": "clicked",
        "value": {},
        "on": [
          { "events": "@cell:click", "update": "clicked && clicked.id === datum.id ? {} : datum" }
        ]
      },
      {
        "name": "active",
        "value": {},
        "update": "clicked",
        "on": [
          { "events": "@cell:mouseover", "update": "clicked.name ? clicked : datum" },
          { "events": "@cell:mouseout", "update": "clicked.name ? clicked : {}" }
        ]
      },
      {
        "name": "freqItems",
        "value": { "fset": "", "supp": 0, "url": "" },
        "update": "active ? {'fset': active.name, 'supp': active.supp, 'url': active.url } : (clickedRule && {'fset': clickedRule.fset, 'supp': clickedRule.supp, 'url': clickedRule.url})",
        "on": [
          { "events": "@ruleNodes:click", "update": "{'fset': datum.rule.fset, 'supp': datum.rule.supp, 'url': datum.rule.url}" },
          { "events": "@ruleLabels:click", "update": "{'fset': datum.rule.fset, 'supp': datum.rule.supp, 'url': datum.rule.url}" }
        ]
      },
      {
        "name": "clickedRule",
        "value": {},
        "update": "clicked && {}",
        "on": [
          { "events": "@ruleNodes:click", "update": "clickedRule.fset && clickedRule.fset === datum.fset ? {} : datum" },
          { "events": "@ruleLabels:click", "update": "clickedRule.fset && clickedRule.fset === datum.fset ? {} : datum" }
        ]
      },
      {
        "name": "activeRule",
        "value": {},
        "update": "clickedRule",
        "on": [
          { "events": "@ruleNodes:mouseover", "update": "(clickedRule.fset && datum.depth <= clickedRule.depth) ? clickedRule : datum" },
          { "events": "@ruleNodes:mouseout", "update": "clickedRule.fset ? clickedRule : {}" },
          { "events": "@ruleLabels:mouseover", "update": "(clickedRule.fset && datum.depth <= clickedRule.depth) ? clickedRule : datum" },
          { "events": "@ruleLabels:mouseout", "update": "clickedRule.fset ? clickedRule : {}" }
        ]
      }
    ],

    "data": [
      {
        "name": "edges",
        "values" : $relStr
        "transform": [
          {
            "type": "formula", "as": "strokeSize",
            "expr": "max(min(datum.size / maxEdgeSize * 10, 7),.3)"
          }
        ]
      },
      {
        "name": "sourceDegree",
        "source": "edges",
        "transform": [
          {"type": "aggregate", "groupby": ["source"]}
        ]
      },
      {
        "name": "targetDegree",
        "source": "edges",
        "transform": [
          {"type": "aggregate", "groupby": ["target"]}
        ]
      },
      {
        "name": "nodes",
        "values": $idStr
        "transform": [
          { "type": "window", "ops": ["rank"], "as": ["order"] },
          {
            "type": "lookup", "from": "sourceDegree", "key": "source",
            "fields": ["id"], "as": ["sourceDegree"],
            "default": {"count": 0}
          },
          {
            "type": "lookup", "from": "targetDegree", "key": "target",
            "fields": ["id"], "as": ["targetDegree"],
            "default": {"count": 0}
          },
          {
            "type": "formula", "as": "degree",
            "expr": "datum.sourceDegree.count + datum.targetDegree.count"
          },
          {
            "type": "formula", "as": "angle",
            "expr": "( 360 * datum.id / dataSize + 270) % 360"
          },
          {
            "type": "formula",
            "expr": "inrange(datum.angle, [90, 270])",
            "as":   "leftside"
          },
          {
            "type": "formula", "as": "x",
            "expr": "originX + radius * cos(PI * datum.angle / 180)"
          },
          {
            "type": "formula", "as": "y",
            "expr": "originX + radius * sin(PI * datum.angle / 180)"

          },
          {
            "type": "voronoi",
            "x": "x",
            "y": "y",
            "size": [{"signal": "height"}, {"signal": "height"}]
          }
        ]
      },
      {
        "name": "selected",
        "source": "edges",
        "transform": [
          {
            "type": "filter",
            "expr": "datum.source === active.id || datum.target === active.id"
          }
        ]
      },
      {
        "name": "activeNode",
        "source": "nodes",
        "transform": [
          {
            "type": "filter",
            "expr": "datum.id === active.id"
          }
        ]
      },
      {
        "name": "rules",
        "values": $ruleStr
      },
      {
        "name": "rulesFilter",
        "source": "rules",
        "transform": [
          {
            "type": "filter",
            "expr": "(datum.lhs === null && datum.rhs === active.name) || datum.lhs === active.name ||  (freqItems.fset && datum.lhs === freqItems.fset) || (freqItems.fset && (datum.lhs + ' |') === slice(freqItems.fset, 0, length(datum.lhs + ' |')))"
          },
          {
            "type": "formula",
            "as": "lhsArr",
            "expr": "datum.lhs ? split(datum.lhs, ' | ') : null"
          },
          {
            "type": "window",
            "groupby": ["lhs"],
            "sort": {"field": "supp", "order": "descending"},
            "ops": ["row_number"],
            "as": ["rank"]
          },
          {
            "type": "formula",
            "as": "hasChildren",
            "expr": "indata('rules', 'lhs', datum.fset)"
          }
        ]
      },
      {
        "name": "tree",
        "source": "rulesFilter",
        "transform": [
          {
            "type": "stratify",
            "key": "fset",
            "parentKey": "lhs"
          },
          {
            "type": "tree",
            "method": "tidy",
            "size": [{"signal": "height-70"}, {"signal": "treeWidth"}],
            "as": ["y", "origX", "depth", "children"]
          },
          {
            "type": "formula",
            "as": "x",
            "expr": "(datum.depth / maxFreqItems) * (treeWidth - 320)"
          },
          {
            "type": "lookup",
            "from": "rulesFilter",
            "key": "fset",
            "fields": ["fset"],
            "as": ["rule"]
          }
        ]
      },
      {
        "name": "links",
        "source": "tree",
        "transform": [
          { "type": "treelinks" },
          {
            "type": "linkpath",
            "shape": "diagonal",
            "orient": "horizontal"
          }
        ]
      }
    ],

    "layout": {
      "padding": 10,
      "offset": 20,
      "bounds": "full",
      "align": "none",
      "center": {"row": true}
    },

    "marks": [

      // BEGIN CHORD DIAGRAM
      {
        "type": "group",
        "name": "chord",
        "title": {
          "text": "Co-Occurences of MeSH Terms in Publications",
          "frame": "group",
          "fontSize": 22
        },

        "encode": {
          "update": {
            "width": { "signal": "chordWidth" },
            "height": { "signal": "height" },
          }
        },

        "marks": [
          {
            "type": "symbol",
            "name": "layout",
            "interactive": false,
            "from": {"data": "nodes"},
            "encode": {
              "enter": {
                "opacity": {"value": 0}
              },
              "update": {
                "x": {"field": "x"},
                "y": {"field": "y"}
              }
            }
          },
          {
            "type": "path",
            "from": {"data": "edges"},
            "encode": {
              "update": {
                "stroke": [
                  {"value": "#3B5B8C"}
                ],
                "strokeOpacity": [
                  {"test": "datum.source === active.id || datum.target === active.id", "value": 1 },
                  {"test": "active.name", "value": 0},
                  {"value": 0.07}
                ],
                "strokeWidth": {"field": "strokeSize"},
                "strokeCap": "round"
              }
            },
            "transform": [
              {
                "type": "lookup", "from": "layout", "key": "datum.id",
                "fields": ["datum.source", "datum.target"],
                "as": ["sourceNode", "targetNode"]
              },
              {
                "type": "linkpath",
                "sourceX": {"expr": "datum.sourceNode.x"},
                "targetX": {"expr": "datum.targetNode.x"},
                "sourceY": {"expr": "datum.sourceNode.y"},
                "targetY": {"expr": "datum.targetNode.y"},
                "shape": "diagonal"
              }
            ]
          },
          {
            "type": "text",
            "name": "textNodes",
            "from": {"data": "nodes"},
            "encode": {
              "enter": {
                "text": {"field": "sname"},
                "baseline": {"value": "middle"}
              },
              "update": {
                "x": {"field": "x"},
                "y": {"field": "y"},
                "dx": {"signal": "textOffset * (datum.leftside ? -1 : 1)"},
                "angle": {"signal": "datum.leftside ? datum.angle - 180 : datum.angle"},
                "align": {"signal": "datum.leftside ? 'right' : 'left'"},
                "fontSize": [
                  {"test": "datum.id === active.id", "signal": "textSize*1.2" },
                  {"signal": "textSize"}
                ],
                "fontWeight": [
                  {"test": "indata('selected', 'source', datum.id)", "value": "bold"},
                  {"test": "indata('selected', 'target', datum.id)", "value": "bold"},
                  {"value": null}
                ],
                "fill": [
                  {"test": "indata('selected', 'source', datum.id)", "value": "black"},
                  {"test": "indata('selected', 'target', datum.id)", "value": "black"},
                  { "value": "LightSlateGray" }
                ]
              }
            }
          },
          {
            "type": "path",
            "name": "cell",
            "from": {"data": "nodes"},
            "encode": {
              "enter": {
                "fill": {"value": "transparent"},
                "strokeWidth": {"value": 0.35}
              },
              "update": {
                "path": {"field": "path"}
              }
            }
          }
        ]
      },

      // BEGIN TREE CODE
      {
        "type": "group",
        "name": "activeTree",
        "title": {
          "text": "Frequent Itemset Tree of Selected Code",
          "frame": "group",
          "fontSize": 22
        },

        "encode": {
          "enter": {
            "width": { "signal": "treeWidth" },
            "height": { "signal": "height" },
          }
        },

        "marks": [
          {
            "type": "path",
            "from": {"data": "links"},
            "encode": {
              "enter": {
                "stroke": {"value": "#ccc"},
                "x": {"value": 160},
                "y": {"value": 30}
              },
              "update": {
                "path": {"field": "path"}
              }
            }
          },
          {
            "type": "symbol",
            "name": "ruleNodes",
            "from": {"data": "tree"},
            "encode": {
              "enter": {
                "size": {"value": 200},
                "fill": {"value": "#3B5B8C"},
                "stroke": {"value": "#3B5B8C"}
              },
              "update": {
                "x": {"signal": "datum.x + 160"},
                "y": {"signal": "datum.y + 30"},
                "fillOpacity": {"signal": "datum.hasChildren ? 1 : 0"}
              },
              "hover": {
                "fillOpacity": {"value": .5}
              }
            }
          },
          {
            "type": "text",
            "name": "ruleLabels",
            "from": {"data": "tree"},
            "encode": {
              "enter": {
                "text": {"field": "sname"},
                "fontSize": {"value": 10},
                "baseline": {"value": "middle"},
                "opacity": {"value": 1}
              },
              "update": {
                "x": {"signal": "datum.x + 160"},
                "y": {"signal": "datum.y + 30"},
                "dx": {"signal": "(datum.depth === 0 ? -1 : 1) * 12"},
                "align": {"signal": "datum.depth === 0 ? 'right' : 'left'"},
                "fontWeight": {"signal": "(datum.children || datum.fset === activeRule.fset) ? 'bold': ''"}
              }
            }
          },

          // Top center frequent item set and support text
          {
            "type": "text",
            "encode": {
              "enter": {
                "fontSize": {"value": 16},
                "baseline": {"value": "middle"},
                "x": {"signal": "treeWidth / 2"},
                "y": {"signal": "5"},
                "align": {"value": "center"},
                "opacity": {"value": 1}
              },
              "update": {
                "text": {"signal": "'{ ' + (activeRule.fset ? activeRule.fset : freqItems.fset) + ' }'"}
              }
            }
          },
          {
            "type": "text",
            "encode": {
              "enter": {
                "fontSize": {"value": 14},
                "baseline": {"value": "middle"},
                "x": {"signal": "treeWidth / 2"},
                "y": {"signal": "20"},
                "align": {"value": "center"},
                "opacity": {"value": 1}
              },
              "update": {
                "text": {"signal": "'Support: ' + (activeRule.supp ? activeRule.supp : freqItems.supp)"}
              }
            }
          },

          // Legend for filled/un-filled tree nodes
          {
            "type": "symbol",
            "encode": {
              "enter": {
                "size": {"value": 200},
                "stroke": {"value": "#3B5B8C"},
                "xc": {"value": 20},
                "yc": {"signal": "height - 30"},
                "fill": {"value": "#3B5B8C"},
                "fillOpacity": {"value": 1}
              }
            }
          },
          {
            "type": "symbol",
            "encode": {
              "enter": {
                "size": {"value": 200},
                "stroke": {"value": "#3B5B8C"},
                "xc": {"value": 20},
                "yc": {"signal": "height - 10"},
                "fillOpacity": {"value": 0}
              }
            }
          },
          {
            "type": "text",
            "encode": {
              "enter": {
                "text": {"value": "Expandable"},
                "fontSize": {"value": 10},
                "x": {"value": 30},
                "y": {"signal": "height - 30"},
                "baseline": {"value": "middle"}
              }
            }
          },
          {
            "type": "text",
            "encode": {
              "enter": {
                "text": {"value": "Not Expandable"},
                "fontSize": {"value": 10},
                "x": {"value": 30},
                "y": {"signal": "height - 10"},
                "baseline": {"value": "middle"}
              }
            }
          },

          // Text to display if frequent item set tree is blank
          {
            "type": "text",
            "encode": {
              "enter": {
                "text": {"value": "Explore the Chord Diagram on the Left to see the Frequent Itemset Tree"},
                "fontSize": {"value": 14},
                "x": {"signal": "treeWidth / 2"},
                "y": {"signal": "height / 2"},
                "baseline": {"value": "middle"},
                "align": {"value": "center"},
                "fill": {"value": "DarkSlateGray"}
              },
              "update": {
                "text": {"signal": "(activeRule.fset || freqItems.fset) ? '' : 'Explore the Chord Diagram on the Left to see the Frequent Item Set Tree'"}
              }
            }
          }
        ]
      }
    ]
  }
  """

  return VegaLite.VGSpec(JSON.parse(spec))

end


text:
module SigmaBenchmarking

using Sigma
using Lens
using Compat
# using DataFrames
# using Gadfly

# Split Functions
import Sigma: weighted_mid_split, rand_partial_split,
              weighted_partial_split
# Solvers
import Sigma: Solver, DRealSolver, DRealSolverBinary, SigmaSolver

using DynamicAnalysis
import DynamicAnalysis: benchmark

benchdir = pwd()

include("algorithms.jl")
include("church.jl")

benchmarks = ["simplex"
              "polyfrompixels"]

println("Running Benchmarks:")

for t in benchmarks
  benchmark_fn = joinpath(t,"$t.jl")
  println(" * $benchmark_fn")
  include(benchmark_fn)
end

"""
Evaluate the performance of Sigma on a bench marks.

- each benchmark is contained within its own folder (e.g. ./simplex)
- each folder contains a julia file of the same name (e.g. simplex/simplex.jl)

"""
SigmaBenchmarking

end

ext: jl


text:
function f!(resid, x)
    for i in eachindex(x)
        resid[i] = sin(x[i]) + x[i]^3
    end
end
x = ones(5)
@test Sundials.kinsol(f!, x, linear_solver = :Band, jac_upper = 0, jac_lower = 0) ==
      Sundials.kinsol(f!, x)


text:
# Length and Step concept from Andrés Riedemann and Mason Protter

export Length, Step

struct Length{T <: Union{Nothing,Integer} }
    length::T
end
Length(value::T) where T = Length{T}(value)
value(length::Length{T}) where T = length.length
Base.length(length::Length{T}) where T = length.length

struct Step{T}
    step::T
    Step(value::T) where T = isone(value) ? new{Nothing}(nothing) : Step{T}(value)
    Step{T}(value::T) where T = new{T}(value)
end
value(step::Step{T}) where T = step.step
Base.step(step::Step{T}) where T = step.step

AnyNothing = (Any, Nothing)

for (A,B,C) in Iterators.product( AnyNothing, AnyNothing, AnyNothing )
    Base._range(start::A, step::Step, stop::B, length::C) = Base._range(start, step.step, stop, length)
    if A != Nothing && C != Nothing
        Base._range(start::A, step::B, stop::C, length::Length) = Base._range(start, step, stop, length.length == stop - start + 1 ? nothing : length.length)
    else
        Base._range(start::A, step::B, stop::C, length::Length) = Base._range(start, step, stop, length.length)
    end
end

for (A,B) in Iterators.product( AnyNothing, AnyNothing )
    if A != Nothing && B != Nothing
        Base._range(start::A, step::Step, stop::B, length::Length) = Base._range(start, step.step, stop, length.length == stop - start + 1 ? nothing : length.length)
    else
        Base._range(start::A, step::Step, stop::B, length::Length) = Base._range(start, step.step, stop, length.length)
    end
end

(::Base.Colon)(start, stop, length::Length) = Base._range(start, nothing, stop, length)
(::Base.Colon)(start::T, length::Length, stop::T) where T<:Real = Base._range(start, nothing, stop, length)
(::Base.Colon)(length::Length, start, stop) = Base._range(start, nothing, stop, length)

(::Base.Colon)(start, stop, ::Length{Nothing}) = Base._range(start, nothing, stop, nothing)
(::Base.Colon)(start::T, ::Length{Nothing}, stop::T) where T<:Real = Base._range(start, nothing, stop, nothing)
(::Base.Colon)(::Length{Nothing}, start, stop) = Base._range(start, nothing, stop, nothing)

(::Base.Colon)(start, stop, step::Step) = start:step.step:stop
(::Base.Colon)(start::T, step::Step, stop::T) where T<:Real = start:step.step:stop
(::Base.Colon)(step::Step, start, stop) = start:step.step:stop

(::Base.Colon)(start, stop, ::Step{Nothing}) = start:stop
(::Base.Colon)(start::T, ::Step{Nothing}, stop::T) where T<:Real = start:stop
(::Base.Colon)(::Step{Nothing}, start, stop) = start:stop

Base.range(start, stop, length::Length) = Base._range(start, nothing, stop, length)
Base.range(start, length::Length, stop) = Base._range(start, nothing, stop, length)
Base.range(length::Length, start, stop) = Base._range(start, nothing, stop, length)

Base.range(start, stop, ::Length{Nothing}) = Base._range(start, nothing, stop, nothing)
Base.range(start, ::Length{Nothing}, stop) = Base._range(start, nothing, stop, nothing)
Base.range(::Length{Nothing}, start, stop) = Base._range(start, nothing, stop, nothing)

Base.range(start, stop, step::Step) = Base._range(start, step, stop, nothing)
Base.range(start, step::Step, stop) = Base._range(start, step, stop, nothing)
Base.range(step::Step, start, stop) = Base._range(start, step, stop, nothing)
Base.range(start::Number, step::Step, stop::T) where {T <:Integer} = Base._range(start, step, stop, nothing)

Base.range(start, stop; length::Length) = Base._range(start, nothing, stop, length)

# Use step and length methods to produce a new range with modifications

"""
    length(r::AbstractRange, len)

Create a range with the same `start` and `stop` but with the specified length, `len`
"""
Base.length(r::AbstractRange, len) = Base.range(first(r), last(r), length=len)
"""
    length(start => stop, len)

Create a range from start to stop and the specified length
"""
Base.length(r::Pair, len) = Base.range(first(r), last(r), length=len)

"""
    length(start, stop, length)

Create a range from start to stop (inclusive) with the specified length
"""
Base.length(start, stop, length) = Base.range(start, stop; length = length)

"""
    step(r::AbstractRange, step)

Return a range with the same `start` and `stop` but with the specified step
"""
Base.step(r::AbstractRange, s) = Base.range(first(r), last(r), step=s)

"""
    step(r::Pair, s)

Return a range with the `start` and `stop` but with the specified step
"""
Base.step(r::Pair, s) = Base.range(first(r), last(r), step=s)
Base.step(start, stop, step) = Base.range(start, stop; step=step)

"""
    first(r::AbstractRange, f)

Modify range with a new `start` element. Does not modify `end` or `step`.
"""
Base.first(r::AbstractRange, f) = Base.range(f, last(r), step=step(r))

"""
    last(r::AbstractRange, l)

Modify range with a new last element. Does not change `start` or `step`.
"""
Base.last(r::AbstractRange, l) = Base.range(first(r), l, step=step(r))


text:
#This script executes the reinforce algorithm from Williams

include("src/init.jl")

θa = init_θ_single_norm(params.structure_net_pl)
θv = init_θ_single_norm(params.structure_net_vf)
best_net = init_θ_single_norm(params.structure_net_pl)
res = zeros(params.n_iter)
prec = zeros(params.n_iter)
all_traj  = init_all_traj(params)

"""
REINFORCE
"""
function reinforce(params::Params, θa::Any, θv::Any, res::Array{Float64,1}, prec::Array{Float64,1}, all_traj::Traj_Data, best_net::Any)
    #Init Part:
    best_net = copy(θa)
    min_cost = Inf
    #------------------------
    #Initialize the Policy
    ĝ_a = copy(θa)
    opt_para_a = optimizers(θa,   Adam; lr = params.lr_pl)

    #Initialize the Value Function
    ĝ_v = copy(θv)
    opt_para_v = optimizers(θv,   Adam; lr = params.lr_vf)

    #Init Progress bar
    progress = Progress(params.n_iter,1)

    # Init rms
    rms_all = init_rms(params.s_init)
    rms_loc = init_rms(params.s_init)

    a_t_old = all_traj.traj[1].a_t
    weights_old = zeros(length(net_to_vec(θa)),2)

    for i_iter = 1:params.n_iter
        ĝ_a *= 0
        ĝ_v *= 0

        #Save data in intervalls
        if floor(i_iter/50) == ceil(i_iter/50)
            save_data(res, prec, θa, θv, all_traj, i_iter)
        end

        #-----------------------------------------------------------------------
        #Solve it in Parallel!
        res[i_iter], prec[i_iter] = parallel_episodes(params,θa,θv,all_traj,rms_all)

        #Apply simple baseline
        simple_baseline!(all_traj::Traj_Data, params::Params)

        #Reorganice Data in minibatches
        mini_batch = build_minibatches(all_traj,params);

        #Apply the gradient
        ĝ_a = parallel_gradient(params, θa, mini_batch)
        Knet.update!(θa, ĝ_a, opt_para_a)

        if min_cost > res[i_iter]
            min_cost = res[i_iter]
            best_net = copy(θa)
            print(" new_best: ", round(min_cost, digits = 3))
        end

        #Visualization_stuff
        #-----------------------------------------------------------------------
        ProgressMeter.update!(progress,i_iter) #update progress bar
        if floor(i_iter/50) == ceil(i_iter/50)
            a_t_old ,weights_old = plot_stuff(params, all_traj, a_t_old,
            weights_old, i_iter,θa, θv, res, prec/params.n_traj, ĝ_a, ĝ_v)
        end
        #-----------------------------------------------------------------------

        #clean the workspace from junk
        @everywhere GC.gc(true)
    end
    return res
end

reinforce(params, θa, θv, res, prec, all_traj, best_net)

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
function baseline_correct!(s::Spectrum, ppm_range::Tuple{Float64,Float64})
    s[:] .-= mean(s[ppm_range])
    s
end


text:
# I implement several simple but super loose bounds
# These only apply to closed tours

# the cost of a tour must be >= the sum over all vertices of
# the cost of the cheapest edge leaving that vertex
# likewise for the cheapest edge entering that vertex
# since we must go to and leave each vertex
function vertwise_bound(distmat::AbstractMatrix{T}) where {T<:Real}
	# the simple code below would tend to pick out the 0 costs on the diagonal
	# so make a doctored copy of the distance matrix with high costs on the diagonal
	m = maximum(distmat)
	distmat_nodiag = distmat + m * I
	leaving = sum(minimum(distmat_nodiag, dims = 2))
	entering = sum(minimum(distmat_nodiag, dims = 1))
	return maximum([leaving, entering])
end

# helper to get min spanning trees
# returns a (n-1) long Vector of Tuple{Int, Int} where each tuple is an edge in the MST
# and the total weight of the tree
# the matrix passed in must be symmetric or you won't get out the minimum spanning tree
function minspantree(dm::AbstractMatrix{T}) where {T<:Real} # accepts views
	mst_edges = Vector{Tuple{Int, Int}}()
	mst_cost = zero(T)
	n = size(dm, 1)

	# we keep a running list of the distance from each vertex to the partly formed tree
	# rather than using 0 for vertices already in the tree, we use a large value so that we
	# can find the closest non-tree vertex via call to Julia's `findmin`.
	bigval = maximum(dm) + one(T)
	tree_dists = dm[1,:] # distance to tree
	closest_tree_verts = ones(Int, n)
	tree_dists[1] = bigval # vert 1 is in tree now

	for _ in 1:(n-1) # need to add n - 1 other verts to tree
		cost, newvert = findmin(tree_dists)
		treevert = closest_tree_verts[newvert]
		# add costs and edges
		mst_cost += cost
		if treevert < newvert
			push!(mst_edges, (treevert, newvert))
		else
			push!(mst_edges, (newvert, treevert))
		end
		# update distances to tree
		tree_dists[newvert] = bigval
		for i in 1:n
			c = tree_dists[i]
			if c >= bigval # already in tree
				continue
			end
			# maybe this vertex is closer to the new vertex than the prior iteration's tree
			if c > dm[i, newvert]
				tree_dists[i] = dm[i, newvert]
				closest_tree_verts[i] = newvert
			end
		end
	end

	return mst_edges, mst_cost
end

# a simplified/looser version of Held-Karp bounds
# any tour is a spanning tree on (n-1) verts plus two edges
# from the left out vert
# so the maximum over all verts (as the left out vert) of
# MST cost on remaining vertices plus 2 cheapest edges from
# the left out vert is a lower bound
# for extra simplicity, the distance matrix is modified to be symmetric so we can treat
# the underlying graph as undirected. This also doesn't help the bound!
function hkinspired_bound(distmat::AbstractMatrix{T}) where {T<:Real}
	n = size(distmat, 1)
	if size(distmat, 2) != n
		error("Must pass square distance matrix to hkinspired_bound")
	end

	# make a symmetric copy of the distance matrix
	distmat = copy(distmat)
	for i in 1:n, j in 1:n
		if i == j
			continue
		end
		d = min(distmat[i,j], distmat[j,i])
		distmat[i,j] = d
		distmat[j,i] = d
	end

	# get a view of the distmat with one vertex deleted
	function del_vert(v)
		keep = [1:(v-1) ; (v+1):n]
		return view(distmat, keep, keep)
	end

	# make sure min(distmat[v,:]) doesn't pick diagonal elements
	m = maximum(distmat)
	distmat_nodiag = distmat + m * I

	# lower bound the optimal cost by leaving a single vertex out
	# forming spanning tree on the rest
	# connecting the left-out vertex
	function cost_leave_out(v)
		dmprime = del_vert(v)
		_, c = minspantree(dmprime)
		c += minimum(distmat_nodiag[v,:])
		c += minimum(distmat_nodiag[:,v])
		return c
	end

	return maximum(map(cost_leave_out, 1:n))
end

# best lower bound we have
"""
    lowerbound(distmat)

Lower bound the cost of the optimal TSP tour. At present, the bounds considered
are a simple bound based on the minimum cost of entering and exiting each city and
a slightly better bound inspired by the Held-Karp bounds; note that the implementation
here is simpler and less tight than proper HK bounds.

The Held-Karp-inspired bound requires computing many spanning trees. For a faster but
typically looser bound, use `TravelingSalesmanHeuristics.vertwise_bound(distmat)`.

!!! note
    The spanning tree bounds are only correct on symmetric problem instances and will not be
    used if the passed `distmat` is not symmetric. If you wish to use these bounds anyway,
    (e.g. for near-symmetric instances), use `TravelingSalesmanHeuristics.hkinspired_bound`.
"""
function lowerbound(distmat::AbstractMatrix{T} where {T<:Real})
	vb = vertwise_bound(distmat)
	if !issymmetric(distmat)
		return vb
	end
	return max(vb, hkinspired_bound(distmat))
end


text:
include("../utils.jl")

@testset "plot_data" begin
    # data / hasmissing / nobj
    x = [5.0, missing, missing, 3.0, 2.0]
    d, h, n = G.plot_data(x)
    @test h == true
    @test n == 1
    @test checkzip(d, hcat(1:length(x), x))

    y       = hcat([2.0, missing, 3.0, Inf, 2.0], [1.0, 2.0, NaN, 0, 1.0])
    d, h, n = G.plot_data(x, y)
    @test h == true
    @test n == 2
    @test checkzip(d, hcat(x, y))

    # -- errors --
    @test_throws ArgumentError G.plot_data([])
    @test_throws ArgumentError G.plot_data([], [1,2])
    @test_throws ArgumentError G.plot_data("hello")
    @test_throws ArgumentError G.plot_data([1,2], "hello")
    @test_throws DimensionMismatch G.plot_data([1,2],[3,4,5])
end

@testset "fill_data" begin
    x = 1:1:10
    y1 = @. exp(x)
    y2 = @. sin(x)
    d, = G.fill_data(x,y1,y2)
    @test checkzip(d, hcat(x,y1,y2))

    # -- errors --
    @test_throws ArgumentError G.fill_data(Real[], Real[], Real[])
    @test_throws DimensionMismatch G.fill_data([0.1,0.2],[0.1,0.2],[0.2])
end

@testset "hist_data" begin
    x = [1, 2, missing, 5]
    d, h, n, r = G.hist_data(x)
    @test checkzip(d, x)
    @test h == true
    @test n == 3
    @test r == (1.0, 5.0)

    @test_throws ArgumentError G.hist_data(Real[])
    @test_throws ArgumentError G.hist_data([missing, missing])
end

@testset "plot" begin
    p = plot([1,2],[3,4])
    @test p isa G.DrawingHandle{<:G.Scatter2D}
    s = p.drawing
    @test checkzip(s.data, [1 3; 2 4])
    @test s.nobj == 1
    @test s.hasmissing == false
end

@testset "fill_between" begin
end

@testset "hist" begin
end

@testset "bar" begin
end

@testset "boxplot" begin
end

@testset "heatmap_ticks" begin
end

@testset "heatmap" begin
end


text:
module TestPyVirtualenv

using PyVirtualenv: _leak, pycall_deps_jl, Py_IsInitialized, activate

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

N = 10000

@testset "_leak(Cstring, ...)" begin
    for i in 1:N
        x = String(rand('A':'z', rand(1:1000)))
        y = Base.unsafe_string(_leak(Cstring, x))
        @test x == y
    end
end

@testset "_leak(Cwstring, ...)" begin
    for i in 1:N
        x = String(rand('A':'z', rand(1:1000)))
        a = Base.cconvert(Cwstring, x)
        ptr = _leak(a)
        z = unsafe_wrap(Array, ptr, size(a))
        @test z[end] == 0
        y = transcode(String, z)[1:end-1]
        @test x == y
    end
end

@testset "activate" begin
    was_inited = Py_IsInitialized()
    pyprogramname = pycall_deps_jl().pyprogramname
    if was_inited
        @warn """
        Python interpreter is already initialized (PyCall is already imported?).
        not testing `activate`.  Use `Pkg.test("PyVirtualenv")` to run all the
        tests.
        """
    else
        @test activate(pyprogramname) == nothing
        @test activate(pyprogramname) == nothing  # it's idempotent
    end
    @test_throws ErrorException activate("non-existing-python",
                                         "non-existing-home")

    local sys_executable
    @eval using PyCall
    sys_executable = pyimport("sys")[:executable]
    if !was_inited
        if Sys.iswindows()
            @info "sys.executable = $sys_executable"
            @test_broken sys_executable == pyprogramname
        else
            @test sys_executable == pyprogramname
        end
    end
end

execdir = joinpath(@__DIR__, "exec")

if Sys.which("pipenv") == nothing
    @warn "`pipenv` command not found. Skipping `test/exec/*.jl`."
else
    @testset for file in sort([file for file in readdir(execdir)
                               if endswith(file, ".jl")])
        script = """
        $(Base.load_path_setup_code(false))
        include($(repr(joinpath(execdir, file))))
        """
        cmd = Base.julia_cmd()
        cmd = `$cmd --startup-file=no --project=$execdir -e $script`
        @debug "Running" cmd
        env = copy(ENV)
        env["PYVIRTUALENV_JL_TEST_TARGET"] = dirname(@__DIR__)
        mktempdir() do path
            cmd = setenv(cmd, env; dir=path)
            @test success(pipeline(cmd; stdout=stdout, stderr=stderr))
        end
    end
end

end  # module


text:
using SymPy

x, y, t = @vars x y t
T = x * (1-x) * y * (1-y)
# K = 1 + x^2 + x / (1+y^2)
K = 1.0
u = x * (1-x)* y * (1-y)
v = x * (1-x) * y * (1-y)
p = x * (1-x) * y * (1-y)
κ1 = 1.0
κ2 = 1.0
w1 = exp(-t) * x * (1-x) * y * (1-y)
w2 = exp(-t) * x * (1-x) * y * (1-y)

Tx = diff(T, x)
Ty = diff(T, y)
k_laplace_T = diff(K * Tx, x) + diff(K * Ty, y)
adv = u * Tx + v * Ty
Q = adv - k_laplace_T
s = sympy.julia_code(simplify(Q))
s = replace(replace(replace(s, ".*"=>"*"), ".^"=>"^"), "./"=>"/")

ffunc = u * diff(u, x) + v * diff(u, y) + diff(p, x) - 0.01 * (diff(diff(u, x), x) + diff(diff(u, y), y))
gfunc = u * diff(v, x) + v * diff(v, y) + diff(p, y) - 0.01 * (diff(diff(v, x), x) + diff(diff(v, y), y))
f = replace(replace(replace(sympy.julia_code(ffunc), ".*"=>"*"), ".^"=>"^"), "./"=>"/")
g = replace(replace(replace(sympy.julia_code(gfunc), ".*"=>"*"), ".^"=>"^"), "./"=>"/")

div_term = diff(u, x) + diff(v, y)
div_term = replace(replace(replace(sympy.julia_code(div_term), ".*"=>"*"), ".^"=>"^"), "./"=>"/")

w1t = diff(w1, t)
w2t = diff(w2, t)
q1_func = w1t  - κ1 * (u-w1)
q1 = replace(replace(replace(sympy.julia_code(q1_func), ".*"=>"*"), ".^"=>"^"), "./"=>"/")
q2_func = w2t  - κ2 * (v-w2)
q2 = replace(replace(replace(sympy.julia_code(q2_func), ".*"=>"*"), ".^"=>"^"), "./"=>"/")


text:
using LatBo.geometry: is_in_pipe, is_in_half_space, is_in_sphere

facts("Check pipe geometries") do
    N = 4  # Measure of the number of points to test

    context("2d") do
        # Direction
        d₀ = Float64[1, 1]
        d⟂ = Float64[1, -1]/√2
        center = Float64[5, 0]
        radius = 3.

        for x in rand(N) * 4., y in radius * [rand(N)..., (1. + rand(N))...]
            r = center + x * d₀ + y * d⟂
            @fact is_in_pipe(r, d₀, center, radius) --> (y <= radius)
        end
    end

    context("3d") do
        # Direction
        d₀ = Float64[1, 1, 1]
        d⟂⁰ = Float64[0, -1, 1]/√2
        d⟂¹ = Float64[-2, 1, 1]/√(4+1+1)
        center = Float64[5, 0, 0]
        radius = 3.
        in_pipe(r) = is_in_pipe(r, d₀, center, radius)

        ys = radius * [rand(N)...]
        zs = radius * [rand(N)..., (1. + rand(N))...]
        for x in rand(N) * 4., y in ys, z in zs
            if norm([z, y]) <= radius
                @fact center + x * d₀ + y * d⟂⁰ + z * d⟂¹ --> in_pipe
            else
                @fact center + x * d₀ + y * d⟂⁰ + z * d⟂¹ --> not(in_pipe)
            end
        end
    end
end

facts("Check half-plane geometries") do
    N = 4  # Measure of the number of points to test
    context("3d") do
        # Direction
        d₀ = Float64[-1, -1, -1]
        d⟂⁰ = Float64[0, -1, 1]/√2
        d⟂¹ = Float64[-2, 1, 1]/√(4+1+1)
        Ω₀ = Float64[4, 3, 1]

        xs = 4.0*[rand(N)..., -rand(N)...]
        ys = 4.0(rand(N) - 0.5)
        zs = 4.0(rand(N) - 0.5)
        for x in xs, y in ys, z in zs
            r = Ω₀ + x * d₀ + y * d⟂⁰ + z * d⟂¹
            inside = is_in_half_space
            @fact inside(r, d₀, Ω₀) --> (x >= 0)
        end
    end
end

facts("Check sphere geometries") do
    N = 4  # Measure of the number of points to test
    context("3d") do
        # Direction
        Ω₀ = Float64[4, 3, 1]
        radius = 3.

        xs = 0.5radius*[rand(N)..., 1. + rand(N)...]
        ys = 0.5radius*[rand(N)..., 1. + rand(N)...]
        zs = 0.5radius*[rand(N)..., 1. + rand(N)...]
        for x in xs, y in ys, z in zs
            r = [x, y, z] + Ω₀
            inside = is_in_sphere
            @fact inside(r, Ω₀, radius) --> (norm([x, y, z]) <= radius)
        end
    end
end


text:
# can define alternate method to change which types are blocked for setproperty! etc
_restricted_fields(::AbstractSample) = (:name, :metadata)

# Samples should at minimum have name and metadata

"""
    name(t::Union{AbstractSample, AbstractFeature})

Get the `name` field from an `AbstractSample` or `AbstractFeature`.
"""
name(as::AbstractSample) = as.name

"""
    metadata(t::AbstractSample)

Get the `metadata` field from an `AbstractSample`.
Note that this is not a copy, so modifications to the returned value
will update the parent `AbstractSample` as well.
"""
metadata(as::AbstractSample) = as.metadata

Base.:(==)(as1::AbstractSample, as2::AbstractSample) = name(as1) == name(as2)

@testset "Sample equality" begin
    as1 = MicrobiomeSample("test")
    as2 = deepcopy(as1)
    as3 = MicrobiomeSample("test2")
    set!(as1, :testvar, 1)
    @test as1 == as2
    @test as1 != as3
    @test as2 != as3
end

Base.String(as::AbstractSample) = name(as)
Base.String(af::AbstractFeature) = name(af)

"""
    getindex(as::AbstractSample, prop::Symbol)

Return the `prop` value in the metadata dictionary of `as`.
This enables using bracket syntax for access, eg `as[prop]`.
"""
function Base.getindex(as::AbstractSample, prop::Symbol)
    prop in _restricted_fields(as) && error("Do not use getindex to access $prop of $(typeof(as)). Use accessor function or getfield instead.")
    getindex(as.metadata, prop)
end

function Base.setindex!(as::AbstractSample, val, prop::Symbol)
    prop in _restricted_fields(as) && error("Do not use setindex! to change $prop of $(typeof(as)).")
    setindex!(as.metadata, val, prop)
    return as
end

"""
    getproperty(as::AbstractSample, prop::Symbol)

Return the `prop` value in the metadata dictionary of `as`.
This enables using dot syntax for access, eg `as.prop`.
"""
function Base.getproperty(as::AbstractSample, prop::Symbol)
    prop in _restricted_fields(as) ? getfield(as, prop) : as.metadata[prop]
end

function Base.setproperty!(as::AbstractSample, prop::Symbol, val)
    prop in _restricted_fields(as) && error("Cannot change $prop of $(typeof(as)) using setproperty!")
    setindex!(as.metadata, val, prop)
    return as
end

"""
    set!(as::AbstractSample, prop::Symbol, val)

Update or insert a value `val` to the metadata of sample `as` using a Symbol `prop`.
If you want an error to be thrown if the value already exists, use [`insert!`](@ref).
"""
function set!(as::AbstractSample, prop::Symbol, val)
    prop in _restricted_fields(as) && error("Cannot set! $prop for $(typeof(as)).")
    set!(as.metadata, prop, val)
    return as
end

"""
    unset!(as::AbstractSample, prop::Symbol)

Delete a metadata entry of sample `as` using the Symbol `prop`.
If you want an error to be thrown if the value does not exist, use [`delete!`](@ref).
"""
function unset!(as::AbstractSample, prop::Symbol)
    prop in _restricted_fields(as) && error("Cannot unset! $prop for $(typeof(as)).")
    unset!(as.metadata, prop)
    return as
end

"""
    insert!(as::AbstractSample, prop::Symbol, val)

Insert a value `val` to the metadata of sample `as` using a Symbol `prop`,
and it will throw an error if `prop` exists.
If you don't want an error to be thrown if the value exists, use [`set!`](@ref).
"""
function insert!(as::AbstractSample, prop::Symbol, val)
    prop in _restricted_fields(as) && error("Cannot insert! $prop for $(typeof(as)).")
    insert!(as.metadata, prop, val)
    return as
end

"""
    delete!(as::AbstractSample, prop::Symbol)

Delete a metadata entry of sample `as` using the Symbol `prop` if it exists, or throw an error otherwise.
If you don't want an error to be thrown if the value does not exist, use [`unset!`](@ref).
"""
function delete!(as::AbstractSample, prop::Symbol)
    prop in _restricted_fields(as) && error("Cannot delete! $prop for $(typeof(as)).")
    delete!(as.metadata, prop)
    return as
end

"""
    keys(as::AbstractSample)

Return an iterator over all keys of the metadata attached to sample `as`.
`collect(keys(as))` returns an array of keys.
"""
Base.keys(as::AbstractSample) = keys(metadata(as))

"""
    haskey(as::AbstractSample, key::Symbol)

Determine whether the metadata of sample `as` has a mapping for a given `key`.
Use `!haskey` to determine whether a sample `as` in a CommunityProfile doesn't have a mapping for a given `key`
"""
Base.haskey(as::AbstractSample, key::Symbol) = in(key, keys(as))

"""
    get(as::AbstractSample, key::Symbol, default)

Return the value of the metadata in the sample `as` stored for the given `key`, or the given `default` value if no mapping for the key is present.
"""
Base.get(as::AbstractSample, key::Symbol, default) = get(metadata(as), key, default)

function set!(as::AbstractSample, d::Union{NamedTuple, Dictionary{Symbol, <:Any}})
    for (key, value) in pairs(d)
        set!(as, key, value)
    end
    return as
end

function insert!(as::AbstractSample, d::Union{NamedTuple, Dictionary{Symbol, <:Any}})
    isempty(Set(keys(as)) ∩ Set(keys(d))) || throw(ArgumentError("Duplicate keys found. Use `set!` to overwrite"))
    for (key, value) in pairs(d)
        insert!(as, key, value)
    end
    return as
end

"""
    MicrobiomeSample(name::String, metadata::Dictionary{Symbol, T}) <: AbstractSample
    MicrobiomeSample(name::String; kwargs...)
    MicrobiomeSample(name::String)

Microbiome sample type that includes a name and a [`Dictionary`](https://github.com/andyferris/Dictionaries.jl)
of arbitrary metadata using `Symbol`s (other than `:name` or `:metadata`) as keys.

Metadata can be accessed using `getproperty` or `getindex` on the sample itself.

Samples can be instantiated with only a name, leaving the `metadata` `Dictionary` blank

Adding or changing metadata follows [the same rules](https://github.com/andyferris/Dictionaries.jl#accessing-dictionaries) as for the normal `Dictionary`.

"""
struct MicrobiomeSample <: AbstractSample
    name::String
    metadata::Dictionary{Symbol, T} where {T <: Any} # currently non-functional
end

MicrobiomeSample(n::AbstractString; kwargs...) = isempty(kwargs) ? MicrobiomeSample(n, Dictionary{Symbol, Any}()) : MicrobiomeSample(n, dictionary(kwargs))
MicrobiomeSample(n::AbstractString, d::Union{AbstractDict,NamedTuple}) = MicrobiomeSample(n; pairs(d)...)


text:
@testset "MatrixFactorisation" begin
    users = repmat(1:30,3)
    items = repmat(1:45,2)
    ratings = StatsBase.sample(1:5,90)
    block = Float64[users items ratings]
    mfmodel = MatrixFactorisationModel(block,20,90);

    gl = DataModel.getgrad(mfmodel)
    ll = DataModel.getllik(mfmodel)

    for i in 1:10
        x = randn(MatrixFactorisation.fetchnparams(mfmodel))
        @test MLUtilities.checkgrad_highdim(x,ll,gl) < 1e-2
    end
end


text:
module OptimalTransport

# MONGE PROBLEM
# -------------

using Combinatorics: permutations

X1 = [-3.334062779243765 -0.4501346178106358;
-17.089568765347632 13.122166467239909;
-14.347484681558905 -6.733978925249632;
4.023051337264342 -9.962304692857852;
-8.463617750898711 -7.272060919414795;
4.411632356371395 18.33766410251966;
0.33231514398722767 -8.943403560908134;
4.312466564646311 -4.769277598451736;
-0.1530439721692118 9.400009756002328;
-5.9829675943508995 0.07220200778660099];

X2 = [-15.776290315139038 -5.26704453200824;
5.035133784625716 6.388525577583841;
-1.302878737641465 -14.087978411169559;
-17.05014605333986 10.189492174622178;
-10.102710138685817 1.6112603598816224;
-17.880416477734084 -5.15108311036524;
-18.49092111472347 -5.193032093476489;
1.2972905950753137 34.72586483583929;
1.3580992685186684 7.334727533023157;
-7.482909025521323 6.326317440680775];

function monge_brute_force(C)
  n, m = size(C)
  @assert n == m "C should be square"
  # loop over all permutations and to find the
  # matching with the lowest cost

  best_cost = typemax(eltype(C))
  for p in permutations(1:n)
      cost = zero(eltype(C))
      for (i, j) in enumerate(p)
        cost += C[i,j]
      end
      if cost < best_cost
        best_cost = cost
        global best_perm = p
      end
    end
  return best_perm, best_cost
end

export X1, X2, monge_brute_force

# SINKHORN
# -------

using LinearAlgebra

function sinkhorn(C::Matrix, a::Vector, b::Vector; λ=1.0, ϵ=1e-8)
    n, m = size(C)
    @assert n == length(a) && m == length(b) throw(DimensionMismatch("a and b do not match"))
    @assert sum(a) ≈ sum(b) "a and b don't have equal sums"
    u, v = copy(a), copy(b)
    M = exp.(-λ * C)
    # normalize this matrix
    while maximum(abs.(a .- Diagonal(u) * (M * v))) > ϵ
        u .= a ./ (M * v)
        v .= b ./ (M' * u)
      end
    return Diagonal(u) * M * Diagonal(v)
  end

export sinkhorn

end


text:
using Cbc

include("MPB_wrapper.jl")
include("MOI_wrapper.jl")


text:
######################################################################
# projection.jl
#
# Handles projections and descents
######################################################################

# returns xdn and udn, the feasible projected trajectory
function project(em::ErgodicManager, tm::TrajectoryManager, K::VMF, xd::VVF, ud::VVF, zd::VVF, vd::VVF, step_size::Float64)
	xdn = [xd[1]]
	udn = Array{VF}(0)

	# perform descent
	alpha = VVF(tm.N + 1)
	for n = 0:tm.N
		alpha[n+1] = xd[n+1] + step_size * zd[n+1]
	end

	# perform the projection
	for n = 1:tm.N
		push!(udn, ud[n] + step_size*vd[n] + K[n]*(alpha[n] - xdn[n]))
		push!(xdn, integrate(tm, xdn[n], udn[n]) )
	end
	return xdn, udn
end

# A projection for LTI systems
function project2(em::ErgodicManager, tm::TrajectoryManager, K::VMF, xd::VVF, ud::VVF, zd::VVF, vd::VVF, step_size::Float64)
	xdn = [xd[1]]
    udn = Array{VF}(0)

	xdn = VVF(0)
	udn = VVF(0)

	# perform the projection
	# Shouldn't need to even integrate...
	for n = 1:tm.N
		push!(udn, ud[n] + step_size*vd[n])
		push!(xdn, xd[n] + step_size*zd[n])
		#push!(xdn, integrate(tm, xdn[n], udn[n]) )
		#push!(xdn, symplectic_euler(tm, xdn[n], udn[n]) )
	end
	push!(xdn, xd[tm.N+1] + step_size*zd[tm.N+1])
	return xdn, udn
end

# No projection and performed in place
function descend!(xd::VVF, ud::VVF, zd::VVF, vd::VVF, step_size::Float64)
	N = length(ud)
	for ni = 1:N
		ud[ni] += step_size * vd[ni]
		xd[ni] += step_size * zd[ni]
	end
	xd[N+1] += step_size * zd[N+1]
end


text:
using Flux, MLDatasets, Statistics, Random, BSON
using Flux.Optimise: update!
using Flux: logitbinarycrossentropy

batchSize, latentDim = 500, 100
epochs = 40
etaD, etaG = 0.0002, 0.0002

images, _ = MLDatasets.MNIST.traindata(Float32)
imageTensor = reshape(@.(2f0 * images - 1f0), 28, 28, 1, :)
data = [imageTensor[:, :, :, r] for r in Iterators.partition(1:60000, batchSize)]

dscr = Chain(Conv((4,4),1=>64;stride=2,pad=1),x->leakyrelu.(x,0.2f0),
        Dropout(0.25),Conv((4,4),64=>128;stride=2,pad=1),x->leakyrelu.(x,0.2f0),
        Dropout(0.25), x->reshape(x, 7 * 7 * 128, :), Dense(7 * 7 * 128, 1))
gen =  Chain(Dense(latentDim,7*7*256),BatchNorm(7*7*256,relu),
        x->reshape(x,7,7,256,:),ConvTranspose((5,5),256=>128;stride=1,pad=2),
        BatchNorm(128,relu),ConvTranspose((4,4),128=>64;stride=2,pad=1),
        BatchNorm(64,relu),ConvTranspose((4,4),64=>1,tanh;stride=2,pad=1))

dLoss(realOut,fakeOut) =    mean(logitbinarycrossentropy.(realOut,1f0)) +
                            mean(logitbinarycrossentropy.(fakeOut,0f0))
gLoss(u) = mean(logitbinarycrossentropy.(u, 1f0))

function updateD!(gen, dscr, x, opt_dscr)
    noise = randn!(similar(x, (latentDim, batchSize)))
    fakeInput = gen(noise)
    ps = Flux.params(dscr)
    loss, back = Flux.pullback(()->dLoss(dscr(x), dscr(fakeInput)), ps)
    grad = back(1f0)
    update!(opt_dscr, ps, grad)
    return loss
end

function updateG!(gen, dscr, x, optGen)
    noise = randn!(similar(x, (latentDim, batchSize)))
    ps = Flux.params(gen)
    loss, back = Flux.pullback(()->gLoss(dscr(gen(noise))),ps)
    grad = back(1f0)
    update!(optGen, ps, grad)
    return loss
end

optDscr, optGen = ADAM(etaD), ADAM(etaG)
cd(@__DIR__)
@time begin
    for ep in 1:epochs
        for (bi,x) in enumerate(data)
            lossD = updateD!(gen, dscr, x, optDscr)
            lossG = updateG!(gen, dscr, x, optGen)
            @info "Epoch $ep, batch $bi, D loss = $(lossD), G loss = $(lossG)"
        end
        @info "Saving generator for epcoh $ep"
        BSON.@save "../data/mnistGAN$(ep).bson" genParams=cpu.(params(gen))
    end
end

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
## simple workspace browser for julia, not the prettiest of things.
using PySide
using DataFrames

## Some functions to work with a module
function get_names(m::Module)
    sort!(map(string, names(m)))
end

unique_id(v::Symbol, m::Module) = isdefined(m,v) ? unique_id(eval(m,v)) : ""
unique_id(x) = string(object_id(x))

## short_summary
## can customize description here
short_summary(x) = summary(x)
short_summary(x::String) = "A string"

## update ids, returning false if the same, true if not
__ids__ = Array(String, 0)
function update_ids(m::Module)
    global __ids__
    nms = get_names(m)
    nms = filter(u -> u != "__ids__", nms)
    a_ids = map(u -> unique_id(symbol(u), m), nms)
    if __ids__ == a_ids
        false
    else
        __ids__ = a_ids
        true
    end
end

MaybeFunction = Union(Nothing, Function)

## get array of names and summaries
## m module
## name_filter Function to filter names by  u -> ismatch(r"^A", u)
## obj_filter Function to filter objects by u -> isa(u, Function)
function get_names_summaries(m::Module, name_filter::MaybeFunction, obj_filter::MaybeFunction)
    nms = get_names(m)

    ## filter by name
    if !isa(name_filter, Nothing)
        nms = filter(name_filter, nms)
    end

    if !isa(obj_filter, Nothing)
        f(nm) = begin
            if !isdefined(m, symbol(nm)) return(false) end
            obj = eval(m, symbol(nm))
            obj_filter(obj)
        end
        nms = filter(f, nms)
    end

    summaries = map(u -> isdefined(m, symbol(u)) ? short_summary(eval(m,symbol(u))) : "undefined", nms)

    if length(nms) == length(summaries)
        return [nms summaries]
    else
        return nothing                  #  didn't match
    end
end

## we use globals to record the filters
global nm_filter = nothing
global obj_filter = nothing

function make_dataframe()
    d = get_names_summaries(Main, nm_filter, obj_filter)
    d = DataFrames.DataFrame(Variable=d[:,1],
                  Variable__BackgroundRole= DataFrames.rep("goldenrod", size(d)[1]),
                  Description=d[:,2])
    d
end

## Simple interface
w = Widget()
w[:resize](500, 400)
setWindowTitle(w, "Workspace browser")
lyt = VBoxLayout(w)
setLayout(w, lyt)

## View of objects
view = TableView(w)
addWidget(lyt, view)

m = DataFrameModel(make_dataframe(), view)
setModel(view, m)

header = view[:horizontalHeader]()
header[:setStretchLastSection](true)

header = view[:verticalHeader]()
header[:hide]()

view[:resizeColumnToContents](0)

function update(view)
    m = DataFrameModel(make_dataframe(), view)
    setModel(view, m)
end

## Name Filter combobox
ed = LineEdit(w)
qconnect(ed, :editingFinished) do
  val = text(ed)
  global nm_filter
  nm_filter = length(val) > 0 ? u -> ismatch(Regex(val), u) : nothing
  update(view)
end
## Object Filter combobox
cb = ComboBox(w)
set_items(cb, ["none", "Function", "PyObject", "Module"])
qconnect(cb, :activated) do idx
   di = Dict([("none", nothing),
              ("Function", u -> isa(u, Function)),
              ("PyObject", u -> isa(u, PyCall.PyObject)),
              ("Module", u -> isa(u, Module))]
         )

    val = get_value(cb)

    global obj_filter
    obj_filter = di[val]
    update(view)
end

flyt = FormLayout(nothing)
flyt[:setAlignment](qt_enum(["AlignLeft", "AlignTop"], how="|"))
addLayout(lyt, flyt)

addRow(flyt, "Name filter:", ed)
addRow(flyt, "Type filter:", cb)

raise(w)

## Or run in timer
## run_timer(::Int32) = run_timer()
## run_timer() = update(view)

## global timer
## timer = Base.TimeoutAsyncWork(run_timer)
## Base.start_timer(timer,int64(200),int64(200))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "TracyProfiler"
version = v"0.7.8"
sources = [
    ArchiveSource("https://github.com/wolfpld/tracy/archive/refs/tags/v$(version).tar.gz",
		  "4021940a2620570ac767eee84e58d572a3faf1570edfaf5309c609752146e950"),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/10.15/MacOSX10.15.sdk.tar.xz",
                  "2408d07df7f324d3beea818585a6d990ba99587c218a3969f924dfcc4de93b62"),
    DirectorySource("./bundled")
]

script = raw"""
mkdir -vp $bindir
mkdir -vp $libdir

cd ${WORKSPACE}/srcdir/tracy*

# Apply patches to disable forcing -march
atomic_patch -p1 ../patches/unix_library_release.patch
atomic_patch -p1 ../patches/unix_common_make.patch
atomic_patch -p1 ../patches/library_extension.patch

# Need full c++17 support so upgrade min osx version and install newer SDK
if [[ "${target}" == x86_64-apple-darwin* ]]; then
    echo "Installing newer MacOS 10.15 SDK"

    pushd $WORKSPACE/srcdir/MacOSX10.*.sdk
    rm -rf /opt/${target}/${target}/sys-root/System
    cp -ra usr/* "/opt/${target}/${target}/sys-root/usr/."
    cp -ra System "/opt/${target}/${target}/sys-root/."
    popd

    # append flag to mkfile for macos target min version
    echo "CFLAGS += -mmacosx-version-min=10.15" >> common/unix-release.mk
    # Disable link-time optimization, we have some mess in our macOS toolchain.
    sed -i 's/ -flto//' {profiler,update,capture,csvexport,import-chrome}/build/unix/release.mk
fi

# Build / install the library
make -j${nproc} -C library/unix release
cp -v "./library/unix/libtracy-release.${dlext}" ${libdir}

# Build / install the profiler GUI
make -j${nproc} -C profiler/build/unix release
cp -v ./profiler/build/unix/Tracy-release* $bindir

# Build / install update utility
make -j${nproc} -C update/build/unix release
cp -v ./update/build/unix/update-release* $bindir

# Build / install capture utility
make -j${nproc} -C capture/build/unix release
cp -v ./capture/build/unix/capture-release* $bindir

# Build / install csvexport utility
make -j${nproc} -C csvexport/build/unix release
cp -v ./csvexport/build/unix/csvexport-release* $bindir

# Build / install import-chrome utility
make -j${nproc} -C import-chrome/build/unix release
cp -v ./import-chrome/build/unix/import-chrome-release* $bindir

install_license LICENSE
"""

# Only supports x86_64 builds for Linux / Unix
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
]
platforms = expand_cxxstring_abis(platforms)

products = Product[
    LibraryProduct("libtracy-release", :libtracy),
    ExecutableProduct("Tracy-release", :tracy_profiler_exe),
    ExecutableProduct("update-release", :tracy_update_exe),
    ExecutableProduct("capture-release", :tracy_capture_exe),
    ExecutableProduct("csvexport-release", :tracy_csvexport_exe),
    ExecutableProduct("import-chrome-release", :tracy_import_chrome_exe),
]

dependencies = [
    Dependency("GLFW_jll"),
    Dependency("FreeType2_jll"),
    Dependency("Capstone_jll"),
    # GTK3 is only needed for the system file dialog on Linux
    Dependency("GTK3_jll", platforms=filter(Sys.islinux, platforms)),
    BuildDependency("Xorg_xorgproto_jll"),
]

# requires std-c++17, full support in gcc 7+, clang 8+
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    julia_compat="1.6",
    preferred_gcc_version=v"8",
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function Coherence(B::Array{Float64,2})
    u=Inf
    for k=1:size(B,2)
        for j=1:(k-1)
            len=(sum(B[:,k].^2).^.5)*(sum(B[:,j].^2).^.5)
            innerProd=abs(sum(B[:,k].*B[:,j]))
            u= minimum([u,innerProd/len])
        end
    end
    return u
end

Coherence(B) = Coherence(convert(Array{Float64,2},B))

max_stars_repo_name: dahlend/CompressedSensing


text:
using Documenter, PMIx

makedocs(sitename="PMIx.jl Documentation")

max_issues_repo_name: JuliaParallel/PMIx

max_forks_repo_name: JuliaParallel/PMIx


text:
#File: efficfron.jl
#Author:
#Date: 29-June-2020

#import XL data
[names,~,returns] = getXLData('hfma_matlab_data','10 Hedge Funds');

N = size(returns,2); ## data points
am1 = fMean(returns,1); #monthly mean
covar = cov(returns); #covariance matrix

#plot mean-variance efficient frontier
portopt(am1,covar,N)


text:
using ChaosTools
using Test, StaticArrays

test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax
let
# Known fixed point locations from publication:
o2x = [0.0
0.0
3.1416
3.1416
1.3346
4.9486]
o2y = [0.0
3.1416
3.1416
0.0
2.6693
3.6139]
o3x = [0.0
0.8121
5.471
0.8121
3.1416
5.471
3.1416
1.7774
4.5058
0.0
0.0
1.7774
6.2832
4.5058
6.2832]
o3y = [0.0
1.6243
2.3295
3.9537
2.3295
4.6589
3.9537
1.7774
4.5058
1.7774
4.5058
3.5548
1.7774
2.7284
4.5058]
ox = [[1.0], o2x, o3x]
oy = [[1.0], o2y, o3y]
println("\nTesting finding stable & unstable fixed points...")
@testset "Standard Map" begin

    ds = Systems.standardmap()
    xs = range(0, stop = 2π, length = 11); ys = copy(xs)
    ics = [SVector{2}(x,y) for x in xs for y in ys]

    # All permutations of [±1, ±1]:
    singss = [[+1, +1], [-1, -1], [+1, -1], [-1, +1]]
    # I know from personal research I only need this `inds`:
    indss = [[1,2]] # <- must be container of vectors!!!
    λs = 0.005 # <- only this allowed to not be vector (could also be vector)

    @testset "order = $o" for o in [2, 3]
        FP = periodicorbits(ds, o, ics, λs, indss, singss; roundtol = 4)
        for fp in FP
            @test round(fp[1], digits = 4) ∈ ox[o]
            @test round(fp[2], digits = 4) ∈ oy[o]
        end
    end
end
end


text:
using MeshCatMechanisms
include(joinpath(pwd(), "src/models/biped_alt.jl"))

# Visualize
include(joinpath(pwd(), "src/models/visualize.jl"))
vis = Visualizer()
open(vis)

urdf = joinpath(pwd(), "src/models/biped/urdf/biped_left_pinned.urdf")
mechanism = parse_urdf(urdf, floating=false)
mvis = MechanismVisualizer(mechanism,
    URDFVisuals(urdf, package_path=[dirname(dirname(urdf))]), vis)

ϵ = 1.0e-8
θ = 10 * pi / 180
h = model.l2 + model.l1 * cos(θ)
ψ = acos(h / (model.l1 + model.l2))
stride = sin(θ) * model.l1 + sin(ψ) * (model.l1 + model.l2)
x1 = [π - θ, π + ψ, θ, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
xT = [π + ψ, π - θ - ϵ, 0.0, θ, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
kinematics(model, x1)[1]
kinematics(model, xT)[1] * 2.0

kinematics(model, x1)[2]
kinematics(model, xT)[2]

q1 = transformation_to_urdf_left_pinned(model, x1[1:5])

set_configuration!(mvis, q1)

qT = transformation_to_urdf_left_pinned(model, xT[1:5])
set_configuration!(mvis, qT)

ζ = 11
xM = [π, π - ζ * pi / 180.0, 0.0, 2.0 * ζ * pi / 180.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
qM = transformation_to_urdf_left_pinned(model, xM[1:5])
set_configuration!(mvis, qM)
kinematics(model, xM)[1]
kinematics(model, xM)[2]

x1_foot_des = kinematics(model, x1)[1]
xT_foot_des = kinematics(model, xT)[1]
xc = 0.5 * (x1_foot_des + xT_foot_des)

x1_foot_des * -1.0 + xT_foot_des
# r1 = x1_foot_des - xc
r1 = xT_foot_des - xc
r2 = 0.1

zM_foot_des = r2

function z_foot_traj(x)
    sqrt((1.0 - ((x - xc)^2.0) / (r1^2.0)) * (r2^2.0))
end

foot_x_ref = range(x1_foot_des, stop = xT_foot_des, length = T)
foot_z_ref = z_foot_traj.(foot_x_ref)

# plot(foot_x_ref,foot_z_ref)

@assert norm(Δ(xT)[1:5] - x1[1:5]) < 1.0e-5

# Horizon
T = 21

tf0 = 0.5
h0 = tf0 / (T-1)

# Bounds
ul, uu = control_bounds(model, T, [-20.0; 0.0], [20.0; 5.0 * h0])

xl, xu = state_bounds(model, T,
    x1 = [Inf * ones(5); x1[1:5]],
    xT = [Inf * ones(5); xT[1:5]])

function c_stage!(c,x,u,t,model)
    c[1] = kinematics(model,x[1:5])[2]
    nothing
end
#
# function c_stage!(c,x,t,model)
#     c[1] = kinematics(model,x[1:5])[2]
#     nothing
# end

m_stage = 1
stage_ineq = (1:m_stage)

include("../src/loop_delta.jl")

# Objective
qq = 0.1*[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]
Q = [t < T ? Diagonal(qq) : Diagonal(qq) for t = 1:T]
R = [Diagonal(1.0e-1*ones(model.nu)) for t = 1:T-1]
c = 1.0
obj = QuadraticTrackingObjective(Q,R,c,
    [xT for t=1:T],[zeros(model.nu) for t=1:T-1])
penalty_obj = PenaltyObjective(1000.0,foot_z_ref,[t for t = 1:T])
multi_obj = MultiObjective([obj,penalty_obj])

# Problem
prob = init_problem(model.nx,model.nu,T,model,multi_obj,
                    xl=xl_traj,
                    xu=xu_traj,
                    ul=[ul*ones(model.nu) for t=1:T-1],
                    uu=[uu*ones(model.nu) for t=1:T-1],
                    hl=[hl for t=1:T-1],
                    hu=[hu for t=1:T-1],
                    general_constraints=true,
                    m_general=model.nx,
                    general_ineq=(1:0),
                    stage_constraints=true,
                    m_stage=[t==1 ? 0 : m_stage for t = 1:T-1],
                    stage_ineq=[t==1 ? (1:0) : stage_ineq for t = 1:T-1])

# MathOptInterface problem
prob_moi = init_MOI_Problem(prob)

# Trajectory initialization
X0 = linear_interp(x1,xT,T) # linear interpolation on state
U0 = [0.001*rand(model.nu) for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0,U0,h0,prob)

# Solve nominal problem
@time Z_nominal_step = solve(prob_moi,copy(Z0),nlp=:SNOPT,max_iter=100,time_limit=120)

# Unpack solutions
X_nominal_step, U_nominal_step, H_nominal_step = unpack(Z_nominal_step,prob)

Q_nominal_step = [X_nominal_step[t][1:5] for t = 1:T]
plot(hcat(Q_nominal_step...)')
foot_traj_nom = [kinematics(model,Q_nominal_step[t]) for t = 1:T]

foot_x_nom = [foot_traj_nom[t][1] for t=1:T]
foot_z_nom = [foot_traj_nom[t][2] for t=1:T]

plt_ft_nom = plot(foot_x_nom,foot_z_nom,aspect_ratio=:equal,xlabel="x",ylabel="z",width=2.0,
    title="Foot 1 trajectory",label="",color=:red)

plot(hcat(U_nominal_step...)',linetype=:steppost)

plot(foot_x_nom)
plot(foot_z_nom)

sum(H_nominal_step)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
module FastIOBuffers

export
    FastWriteBuffer,
    FastReadBuffer

export
    setdata!

include("fastwritebuffer.jl")
include("fastreadbuffer.jl")

end # module


text:
export BrickletRemoteSwitchIdentity
struct BrickletRemoteSwitchIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletRemoteSwitch
"""
Controls remote mains switches
"""
mutable struct BrickletRemoteSwitch <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletRemoteSwitch(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_remote_switch")
        deviceInternal = package.BrickletRemoteSwitch(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export switch_socket
"""
    $(SIGNATURES)

This function is deprecated, use :func:`Switch Socket A` instead.
"""
function switch_socket(device::BrickletRemoteSwitch, house_code, receiver_code, switch_to)
    device.deviceInternal.switch_socket(house_code, receiver_code, switch_to)
end

export get_switching_state
"""
    $(SIGNATURES)

Returns the current switching state. If the current state is busy, the
Bricklet is currently sending a code to switch a socket. It will not
accept any requests to switch sockets until the state changes to ready.

How long the switching takes is dependent on the number of repeats, see
:func:`Set Repeats`.
"""
function get_switching_state(device::BrickletRemoteSwitch)
    return device.deviceInternal.get_switching_state()
end

export set_repeats
"""
    $(SIGNATURES)

Sets the number of times the code is sent when one of the switch socket
functions is called. The repeats basically correspond to the amount of time
that a button of the remote is pressed.

Some dimmers are controlled by the length of a button pressed,
this can be simulated by increasing the repeats.
"""
function set_repeats(device::BrickletRemoteSwitch, repeats)
    device.deviceInternal.set_repeats(repeats)
end

export get_repeats
"""
    $(SIGNATURES)

Returns the number of repeats as set by :func:`Set Repeats`.
"""
function get_repeats(device::BrickletRemoteSwitch)
    return device.deviceInternal.get_repeats()
end

export switch_socket_a
"""
    $(SIGNATURES)

To switch a type A socket you have to give the house code, receiver code and the
state (on or off) you want to switch to.

A detailed description on how you can figure out the house and receiver code
can be found :ref:`here <remote_switch_bricklet_type_a_house_and_receiver_code>`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function switch_socket_a(device::BrickletRemoteSwitch, house_code, receiver_code, switch_to)
    device.deviceInternal.switch_socket_a(house_code, receiver_code, switch_to)
end

export switch_socket_b
"""
    $(SIGNATURES)

To switch a type B socket you have to give the address, unit and the state
(on or off) you want to switch to.

To switch all devices with the same address use 255 for the unit.

A detailed description on how you can teach a socket the address and unit can
be found :ref:`here <remote_switch_bricklet_type_b_address_and_unit>`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function switch_socket_b(device::BrickletRemoteSwitch, address, unit, switch_to)
    device.deviceInternal.switch_socket_b(address, unit, switch_to)
end

export dim_socket_b
"""
    $(SIGNATURES)

To control a type B dimmer you have to give the address, unit and the
dim value you want to set the dimmer to.

A detailed description on how you can teach a dimmer the address and unit can
be found :ref:`here <remote_switch_bricklet_type_b_address_and_unit>`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function dim_socket_b(device::BrickletRemoteSwitch, address, unit, dim_value)
    device.deviceInternal.dim_socket_b(address, unit, dim_value)
end

export switch_socket_c
"""
    $(SIGNATURES)

To switch a type C socket you have to give the system code, device code and the
state (on or off) you want to switch to.

A detailed description on how you can figure out the system and device code
can be found :ref:`here <remote_switch_bricklet_type_c_system_and_device_code>`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function switch_socket_c(device::BrickletRemoteSwitch, system_code, device_code, switch_to)
    device.deviceInternal.switch_socket_c(system_code, device_code, switch_to)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletRemoteSwitch)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletRemoteSwitch, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using DSGE
using Test, BenchmarkTools
using JLD2

import DSGE: klein_transition_matrices, n_model_states, n_backward_looking_states

# What do you want to do?
check_steady_state = true
check_jacobian = true
check_solution = true

path = dirname(@__FILE__)

m = RealBondMkup()

# Steady-state computation
if check_steady_state
    steadystate!(m)
    # @btime steadystate!(m)

    file = jldopen("$path/reference/steady_state.jld2", "r")
    saved_ell  = read(file, "ell")
    saved_c    = read(file, "c")
    saved_η    = read(file, "eta")
    saved_μ    = read(file, "mu")
    saved_β    = read(file, "beta")
    saved_χss  = read(file, "chi_ss")
    close(file)

    @testset "Check steady state outputs" begin
        @test saved_ell ≈ m[:lstar].value
        @test saved_c   ≈ m[:cstar].value
        @test saved_η   ≈ m[:ηstar].value
        @test saved_μ   ≈ m[:μstar].value
        @test saved_β   ≈ m[:βstar].value
        @test saved_χss ≈ m[:χstar].value
    end
end

# Jacobian computation
if check_jacobian
    m.testing = true        # So that it will test against the unnormalized Jacobian
    JJ = DSGE.jacobian(m)
    # @btime JJ = DSGE.jacobian(m)

    file = jldopen("$path/reference/jacobian.jld2", "r")
    saved_JJ  = read(file, "JJ")

    MUP = read(file, "MUP")
    LIIP = read(file, "LIIP")
    ZP  = read(file, "ZP")
    MONP = read(file, "MONP")
    MKPP = read(file, "MKPP")
    ELLP = read(file, "ELLP")
    RRP = read(file, "RRP")
    IIP = read(file, "IIP")
    WWP = read(file, "WWP")
    PIP = read(file, "PIP")
    TTP = read(file, "TTP")
    MU  = read(file, "MU")
    LII = read(file, "LII")
    Z   = read(file, "Z")
    MON = read(file, "MON")
    MKP = read(file, "MKP")
    ELL = read(file, "ELL")
    RR  = read(file, "RR")
    II  = read(file, "II")
    WW  = read(file, "WW")
    PI  = read(file, "PI")
    TT  = read(file, "TT")

    F1 = read(file, "F1")
    F2 = read(file, "F2")
    F3 = read(file, "F3")
    F4 = read(file, "F4")
    F5 = read(file, "F5")
    F6 = read(file, "F6")
    F7 = read(file, "F7")
    F8 = read(file, "F8")
    F9 = read(file, "F9")
    F10 = read(file, "F10")
    F11 = read(file, "F11")
    close(file)

   #= jldopen("$path/reference/jacobian.jld2", "w") do file
        file["JJ"]   = JJ
        file["MUP"]  = MUP
        file["LIIP"] = LIIP
        file["ZP"]   = ZP
        file["MONP"] = MONP
        file["MKPP"] = MKPP
        file["ELLP"] = ELLP
        file["RRP"]  = RRP
        file["IIP"]  = IIP
        file["WWP"]  = WWP
        file["PIP"]  = PIP
        file["TTP"]  = TTP
        file["MU"]   = MU
        file["LII"]  = LII
        file["Z"]    = Z
        file["MON"]  = MON
        file["MKP"]  = MKP
        file["ELL"]  = ELL
        file["RR"]   = RR
        file["II"]   = II
        file["WW"]   = WW
        file["PI"]   = PI
        file["TT"]   = TT

        file["F1"]   = F1
        file["F2"]   = F2
        file["F3"]   = F3
        file["F4"]   = F4
        file["F5"]   = F5
        file["F6"]   = F6
        file["F7"]   = F7
        file["F8"]   = F8
        file["F9"]   = F9
        file["F10"]  = F10
        file["F11"]  = F11
    end =#

    @testset "Check jacobian outputs" begin
        @testset "Euler Equation" begin
            @test saved_JJ[F1, ELLP] ≈ JJ[F1, ELLP]
            @test saved_JJ[F1, RRP]  ≈ JJ[F1, RRP]
            @test saved_JJ[F1, WWP]  ≈ JJ[F1, WWP]
            @test saved_JJ[F1, TTP]  ≈ JJ[F1, TTP]
            @test saved_JJ[F1, ELL]  ≈ JJ[F1, ELL]
            @test saved_JJ[F1, RR]   ≈ JJ[F1, RR]
            @test saved_JJ[F1, WW]   ≈ JJ[F1, WW]
            @test saved_JJ[F1, TT]   ≈ JJ[F1, TT]
        end

        @testset "Kolmogorov Forward Equation" begin
            @test saved_JJ[F2, MUP] ≈ JJ[F2, MUP]
            @test saved_JJ[F2, MU]  ≈ JJ[F2, MU]
            @test saved_JJ[F2, ELL] ≈ JJ[F2, ELL]
            @test saved_JJ[F2, WW]  ≈ JJ[F2, WW]
            @test saved_JJ[F2, RR]  ≈ JJ[F2, RR]
            @test saved_JJ[F2, TT]  ≈ JJ[F2, TT]
        end

        @testset "Market clearing condition" begin
            @test saved_JJ[F3, MU]  ≈ JJ[F3, MU]
            @test saved_JJ[F3, Z]   ≈ JJ[F3, Z]
            @test saved_JJ[F3, ELL] ≈ JJ[F3, ELL]
            @test saved_JJ[F3, WW]  ≈ JJ[F3, WW]
            @test saved_JJ[F3, RR]  ≈ JJ[F3, RR]
            @test saved_JJ[F3, TT]  ≈ JJ[F3, TT]
        end

        @testset "Technology process" begin
            @test saved_JJ[F4, ZP] ≈ JJ[F4, ZP]
            @test saved_JJ[F4, Z] ≈ JJ[F4, Z]
        end

        @testset "Phillips equation" begin
            @test saved_JJ[F5, PIP] ≈ JJ[F5, PIP]
            @test saved_JJ[F5, Z]   ≈ JJ[F5, Z]
            @test saved_JJ[F5, WW]  ≈ JJ[F5, WW]
            @test saved_JJ[F5, PI]  ≈ JJ[F5, PI]
        end

        @testset "Taylor equation" begin
            @test saved_JJ[F6, II]  ≈ JJ[F6, II]
            @test saved_JJ[F6, PI]  ≈ JJ[F6, PI]
            @test saved_JJ[F6, MON] ≈ JJ[F6, MON]
        end

        @testset "Fisher equation" begin
            @test saved_JJ[F7, II]  ≈ JJ[F7, II]
            @test saved_JJ[F7, PIP] ≈ JJ[F7, PIP]
            @test saved_JJ[F7, RR]  ≈ JJ[F7, RR]
        end

        @testset "Transfers equation" begin
            @test saved_JJ[F8, TT] ≈ JJ[F8, TT]
            @test saved_JJ[F8, Z]  ≈ JJ[F8, Z]
            @test saved_JJ[F8, WW] ≈ JJ[F8, WW]
        end

        @testset "Monetary policy shock equation" begin
            @test saved_JJ[F9, MONP] ≈ JJ[F9, MONP]
            @test saved_JJ[F9, MON]  ≈ JJ[F9, MON]
        end

        @testset "Markup shock equation" begin
            @test saved_JJ[F10, MKPP] ≈ JJ[F10, MKPP]
            @test saved_JJ[F10, MKP]  ≈ JJ[F10, MKP]
        end

        @testset "Lagged monetary policy equation" begin
            @test saved_JJ[F11, LIIP] ≈ JJ[F11, LIIP]
            @test saved_JJ[F11, II]   ≈ JJ[F11, II]
        end
    end
end

# Solve
if check_solution
    m.testing = false      # So the Jacobian will be normalized within the klein solution
    gx, hx = klein(m)
    # @btime klein(m)

    #=JLD2.jldopen("$path/reference/solve.jld2", "w") do file
        file["gx"] = gx
        file["hx"] = hx
    end=#

    file = JLD2.jldopen("$path/reference/solve.jld2", "r")
    saved_gx   = read(file, "gx")
    saved_hx   = read(file, "hx")
    close(file)

    @testset "Check solve outputs" begin
        @test saved_gx  ≈ gx
        @test saved_hx  ≈ hx
    end
end


text:
# push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))
using JLD2
filepath = pwd()
# jlfile = jldopen("../data/nearly_local.jld2", "a+")
jlfile = jldopen("data/nonlocal_less_diffusivity.jld2", "a+")
# jlfile = jldopen("data/nonlocal.jld2", "a+")
# filepath = pwd()
# jlfile = jldopen("data/nearly_local.jld2", "a+")
ψ¹ = jlfile["streamfunction"]["ψ¹"]
ψ² = jlfile["streamfunction"]["ψ²"]
ψ³ = jlfile["streamfunction"]["ψ³"]
ψ⁴ = jlfile["streamfunction"]["ψ⁴"]
z = jlfile["grid"]["z"][:]
x = jlfile["grid"]["x"][:]

## Plot it
using GLMakie
fig = Figure(resolution = (1800, 1300), title = "Local Operators")
titlestring = "ψ¹"
ax1 = Axis(fig[1, 1], title = titlestring, titlesize = 30)
titlestring = "ψ²"
ax2 = Axis(fig[1, 3], title = titlestring, titlesize = 30)
titlestring = "ψ³"
ax3 = Axis(fig[2, 1], title = titlestring, titlesize = 30)
titlestring = "ψ⁴"
ax4 = Axis(fig[2, 3], title = titlestring, titlesize = 30)

colormap = :balance
hm1 = heatmap!(ax1, x, z, ψ¹, colormap = colormap, interpolate = true)

hm2 = heatmap!(ax2, x, z, ψ², colormap = colormap, interpolate = true)

hm3 = heatmap!(ax3, x, z, ψ³, colormap = colormap, interpolate = true)

hm4 = heatmap!(ax4, x, z, ψ⁴, colormap = colormap, interpolate = true)

Colorbar(fig[1, 2], hm1, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize = 25, tickalign = 1,)
Colorbar(fig[2, 2], hm2, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize = 25, tickalign = 1,)
Colorbar(fig[1, 4], hm3, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize = 25, tickalign = 1,)
Colorbar(fig[2, 4], hm4, height = Relative(3 / 4), width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize = 25, tickalign = 1,)
display(fig)

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
function printsubtypes(T, indent=0)
    sT = subtypes(T)
    println(" "^indent, T, isempty(sT) ? "" : ":")
    for S in sT
        printsubtypes(S, indent + 1)
    end
end

function supertypes(T)
    print(T)
    if T != Any
        print(" <: ")
        S = supertype(T)
        supertypes(S)
    end
end

max_stars_repo_path: Chapter06/05. Understanding the structure of Julia numeric types with introspection/types.jl


text:
using Distributions
import Distributions.logpdf
using Base.Cartesian
using Base.Threads
using FillArrays

export logpdf
export rand

export For
struct For{F,T,D,X}
    f :: F
    θ :: T
end

#########################################################
# T <: NTuple{N,J} where {J <: Integer}
#########################################################

For(f, θ::J...) where {J <: Integer} = For(f,θ)

function For(f::F, θ::T) where {F, N, J <: Integer, T <: NTuple{N,J}}
    d = f.(Ones{Int}(N)...)
    D = typeof(d)
    X = eltype(d)
    For{F, NTuple{N,J}, D, X}(f,θ)
end

@inline function logpdf(d::For{F,T,D,X1},xs::AbstractArray{X2,N}) where {F, N, J <: Integer, T <: NTuple{N,J}, D,  X1,  X2}
    s = 0.0
    @inbounds @simd for θ in CartesianIndices(d.θ)
        s += logpdf(d.f(Tuple(θ)...), xs[θ])
    end
    s
end

function Base.rand(dist::For)
    map(CartesianIndices(dist.θ)) do I
        (rand ∘ dist.f)(Tuple(I)...)
    end
end

#########################################################
# T <: NTuple{N,J} where {J <: AbstractUnitRange}
#########################################################

For(f, θ::J...) where {J <: AbstractUnitRange} = For(f,θ)

function For(f::F, θ::T) where {F, N, J <: AbstractRange, T <: NTuple{N,J}}
    d = f.(ones(Int, N)...)
    D = typeof(d)
    X = eltype(d)
    For{F, NTuple{N,J}, D, X}(f,θ)
end

@inline function logpdf(d::For{F,T,D,X1},xs::AbstractArray{X2,N}) where {F, N, J <: AbstractRange,  T <: NTuple{N,J}, D, X1, X2}
    s = 0.0
    @inbounds @simd for θ in CartesianIndices(d.θ)
        s += logpdf(d.f(Tuple(θ)...), xs[θ])
    end
    s
end

function Base.rand(dist::For{F,T}) where {F,  N, J <: AbstractRange, T <: NTuple{N,J}}
    map(CartesianIndices(dist.θ)) do I
        (rand ∘ dist.f)(Tuple(I)...)
    end
end

#########################################################
# T <: Base.Generator
#########################################################

function For(f::F, θ::T) where {F, T <: Base.Generator}
    d = f(θ.f(θ.iter[1]))
    D = typeof(d)
    X = eltype(d)
    For{F, T, D, X}(f,θ)
end

@inline function logpdf(d :: For{F,T}, x) where {F,T <: Base.Generator}
    s = 0.0
    for (θj, xj) in zip(d.θ, x)
        s += logpdf(d.f(θj), xj)
    end
    s
end

@inline function rand(d :: For{F,T,D,X}) where {F,T <: Base.Generator, D, X}
    rand.(Base.Generator(d.θ.f, d.θ.iter))
end

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

export logpdf2
@inline function logpdf2(d::For{F,N,X1},xs) where {F,N, X1, X2}
    results = zeros(eltype(xs), nthreads())

    θ = CartesianIndices(d.θ)

    total = Threads.Atomic{Float64}(0.0)
    @threads for tid in 1:nthreads()
        # split work
        start = 1 + ((tid - 1) * length(xs)) ÷ nthreads()
        stop = (tid * length(xs)) ÷ nthreads()
        domain = start:stop

        s = 0.0
        for j in domain
            @inbounds θj = θ[j]
            @inbounds s += logpdf(d.f(Tuple(θj)...), xs[θj])
        end

        Threads.atomic_add!(total, s)
    end

    total.value
end

@inline function importanceSample(p::For{F1,N,X}, q::For{F2,N,X}) where {F1,F2,N,X}
    _ℓ = 0.0
    x = Array{X, N}(undef, length.(q.θ))
    @inbounds @simd for θ in CartesianIndices(q.θ)
        I = Tuple(θ)
        ℓx = importanceSample(p.f(I...), q.f(I...))
        _ℓ += ℓx.ℓ
        x[θ] = ℓx.val
    end
    Weighted(_ℓ,x)
end

# using Transducers
# using Transducers: @next, complete
# function Transducers.__foldl__(rf, val, d::For)
#     for θ in d.θ
#         val = @next(rf, val, f(θ))
#     end
#     return complete(rf, val)
# end


text:
function row2data(row)
    i, a = match(r"(.+) \(contains (.+)\)", row).captures
    return split(i, " "), split(a, ", ")
end

data = open("allergen.txt") do f
    map(r -> row2data(r), readlines(f))
end

allergens = reduce(vcat, map(r -> r[2], data))
ingredients = reduce(vcat, map(r -> r[1], data))

ingredient_candidates = []
for a in unique(allergens)
    i = intersect(map(r -> r[1], filter(r -> a in r[2], data))...)
    push!(ingredient_candidates, a => i)
end

ingredient_counts = Dict([(i, count(x -> x == i, ingredients)) for i in unique(ingredients)])

for i in unique(Iterators.flatten(map(x -> x[2], ingredient_candidates)))
    delete!(ingredient_counts, i)
end
println(sum(values(ingredient_counts)))

while !all(x -> length(x[2]) == 1, ingredient_candidates)
    solved = collect(Iterators.flatten(map(x -> x[2], filter(x -> length(x[2]) == 1, ingredient_candidates))))
    for (idx, i) in enumerate(ingredient_candidates)
        if length(i[2]) == 1
            continue
        end
        ingredient_candidates[idx] = i[1] => filter(x -> x ∉ solved, i[2])
    end
end

println(join(map(x -> x[2][1], sort(ingredient_candidates, by=x -> x[1])),","))


text:
using Documenter, WrapperSimplification

makedocs(;
    modules=[WrapperSimplification],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/oxinabox/WrapperSimplification.jl/blob/{commit}{path}#L{line}",
    sitename="WrapperSimplification.jl",
    authors="Lyndon White",
    assets=String[],
)

deploydocs(;
    repo="github.com/oxinabox/WrapperSimplification.jl",
)


text:
using ExcelFiles
using Mimi
using Random
using Statistics
using Test

include("../../wip/load_empirical_dist.jl")

filename = joinpath(@__DIR__, "RB-ECS-distribution.xls")
d = load_empirical_dist(filename, "Sheet1!A2:A1001", "Sheet1!B2:B1001")

# Set the seed to get repeatable results (with some caveats...)
Random.seed!(1234567)

# Get the statistical outline of the distribution
q = quantile.(Ref(d), [0.01, 0.10, 0.25, 0.50, 0.75, 0.90, 0.99])

expected = [1.47, 1.92, 2.33, 3.01, 4.16, 5.85, 9.16]
# @info "quantiles: $q\n expected: $expected"

@test isapprox(q, expected, atol=0.01)

# Test that the EmpiricalDistribution gets saved as a SampleStore and values
# get re-used across multiple scenarios

num_scenarios = 4
num_trials = 5
output_dir = "./out/"

results = zeros(num_scenarios, num_trials)

_values = collect(1:10)

# try providing probabilities and also default uniform distribution of probabilities
for _probs in [0.1 * ones(10), nothing]

    _probs = 0.1 * ones(10)

    sd = @defsim begin
        sampling(LHSData)
        test.p = EmpiricalDistribution(_values, _probs)
    end

    @defcomp test begin
        p = Parameter(default = 5)
        function run_timestep(p, v, d, t) end
    end

    scenario_args = [
        :num => collect(1:num_scenarios)
    ]

    function scenario_func(sim_inst::SimulationInstance{T}, tup::Tuple) where T <: AbstractSimulationData
        nothing
    end

    function post_trial_func(sim_inst::SimulationInstance{T}, trialnum::Int, ntimesteps::Int,
                            tup::Tuple)  where T <: AbstractSimulationData
        m, = sim_inst.models
        scenario_num, = tup
        results[scenario_num, trialnum] = m[:test, :p]
    end

    m = Model()
    set_dimension!(m, :time, 2000:2001)
    add_comp!(m, test)

    si = run(sd, m, num_trials;
        results_output_dir = output_dir,
        scenario_args = scenario_args,
        scenario_func = scenario_func,
        post_trial_func = post_trial_func
    )

    for rv in values(si.sim_def.rvdict)
        @test rv.dist isa Mimi.SampleStore
    end

    for i = 1:num_scenarios
        @test results[i, :] == results[1, :]
    end

    rm(output_dir, recursive = true)

    si1 = run(sd, m, num_trials)
    trial1 = copy(collect(values(si1.sim_def.rvdict))[1].dist.values)

    for rv in values(si1.sim_def.rvdict)
        @test rv.dist isa Mimi.SampleStore
    end

    si2 = run(sd, m, num_trials)
    trial2 = copy(collect(values(si2.sim_def.rvdict))[1].dist.values)

    @test trial1!=trial2
end


text:
# This file is a part of ForwardDiffPullbacks.jl, licensed under the MIT License (MIT).

using ForwardDiffPullbacks
using Test

using ChainRulesCore

include("testfuncs.jl")

@testset "ChainRulesCore" begin
    @test @tinferred(ChainRulesCore.rrule(fwddiff(f), xs...)) isa Tuple{Tuple, Function}
    @test @tinferred((ChainRulesCore.rrule(fwddiff(f), xs...)[2])(ΔΩ)) isa Tuple{NoTangent, ForwardDiffPullbacks.FwdDiffPullbackThunk, ForwardDiffPullbacks.FwdDiffPullbackThunk,ForwardDiffPullbacks.FwdDiffPullbackThunk}
    @test @tinferred(map(unthunk, (ChainRulesCore.rrule(fwddiff(f), xs...)[2])(ΔΩ))) == (NoTangent(), 280, (600, 1040), SVector(1600, 2280, 3080))
end


text:
using Test, BalancedIncompleteBlockDesigns

n = 3
b = n*n+n+1
v = b
k = n+1
r = n+1
l = 1

A = prime_plane(3)
@test BIBD_check(A) == (b,v,r,k,l)