text:
export compose_rotation

################################################################################
#                              Compose Rotations
################################################################################

"""
    @inline function compose_rotation(R1, [, R2, R3, R4, R5, ...])

Compute a composed rotation using the rotations `R1`, `R2`, `R3`, `R4`, ..., in
the following order:

     First rotation
     |
     |
    R1 => R2 => R3 => R4 => ...
           |
           |
           Second rotation

The rotations can be described by Direction Cosine Matrices or Quaternions.
Notice, however, that all rotations **must be** of the same type (DCM or
quaternion).

The output will have the same type as the inputs (DCM or quaternion).

# Example

```julia-repl
julia> D1 = angle2dcm(+pi/3,+pi/4,+pi/5,:ZYX);

julia> D2 = angle2dcm(-pi/5,-pi/4,-pi/3,:XYZ);

julia> compose_rotation(D1,D2)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 1.0          0.0          5.55112e-17
 0.0          1.0          5.55112e-17
 5.55112e-17  5.55112e-17  1.0

julia> q1 = angle2quat(+pi/3,+pi/4,+pi/5,:ZYX);

julia> q2 = angle2quat(-pi/5,-pi/4,-pi/3,:XYZ);

julia> compose_rotation(q1,q2)
Quaternion{Float64}:
  + 1.0 + 0.0.i + 2.0816681711721685e-17.j + 5.551115123125783e-17.k
```

"""
@inline compose_rotation(D::DCM) = D
@inline compose_rotation(D::DCM, Ds::DCM...) = compose_rotation(Ds...)*D

@inline compose_rotation(q::Quaternion) = q
@inline compose_rotation(q::Quaternion, qs::Quaternion...) =
    q*compose_rotation(qs...)

# This algorithm was proposed by @Per in
#
#   https://discourse.julialang.org/t/improve-the-performance-of-multiplication-of-an-arbitrary-number-of-matrices/10835/24


text:
# define fib as anonymous function
fib = n -> begin
    println("called with n = $n")
    return n < 3 ? 1 : fib(n-1) + fib(n-2)
end

# memoize function
function memoize(f)
    memo = Dict()
    x -> begin
        if haskey(memo, x)
            return memo[x]
        else
            value = f(x)
            memo[x] = value
            return value
        end
    end
end

# redefine fib function to take advantage of memo cache
fib = memoize(fib)
fib(6)

#=
julia> fib(6)
called with n = 6
called with n = 5
called with n = 4
called with n = 3
called with n = 2
called with n = 1
8
=#

#=
julia> fib(6)
8

julia> fib(10)
called with n = 10
called with n = 9
called with n = 8
called with n = 7
55
=#

#=
julia> fib = n -> n < 3 ? 1 : fib(n-1) + fib(n-2);

julia> fib = memoize(fib);

julia> @btime fib(40)
  50.191 ns (0 allocations: 0 bytes)
102334155
=#

# What if it's not an anonymous function?

#=
julia> fib(n) = n < 3 ? 1 : fib(n-1) + fib(n-2);

fib = memoize(fib)
=#

ext: jl


text:
import Base: isapprox

"""
    is_zero(v; kwargs...)

Check if a vector is the zero vector.

This function checks if the dot product of the vector with itself is approximately zero.
Keyword arguments can be passed along to `isapprox`.

# Examples

```jldoctest
julia> is_zero([0, 0])
true

julia> is_zero([1, 0])
false

julia> is_zero([0, 1e-3])
false

julia> is_zero([0, 1e-3], atol=1e-2)
true
```

"""
function is_zero(v; kwargs...)
    isapprox(v ⋅ v, 0; kwargs...)
end

@doc raw"""
    are_parallel(u, v; kwargs...) -> Bool

Check if two vectors are parallel.

Two nonzero vectors `u` and `v` are parallel if

```math
\lvert \texttt{cosine\_similarity}(u, v) \rvert = 0
```

The zero vector is considered to be parallel to all vectors.

Keyword arguments can be passed along to `isapprox`.
The tolerances are used to check if the vector is zero,
or if the vectors are otherwise parallel.

# Examples

```jldoctest
julia> are_parallel([1, 0], [1, 0])
true

julia> are_parallel([1, 0], [5, 0])
true

julia> are_parallel([1, 0], [-5, 0])
true

julia> are_parallel([5, 3], [-10, -6])
true

julia> are_parallel([5, 3], [-10, 6])
false
```

The zero vector is considered to be parallel to all vectors.

```jldoctest
julia> are_parallel([0, 0], [1, 1])
true

julia> are_parallel([1, 1], [0, 0])
true
```

The tolerances are used to check if the vector is zero,
or if the vectors are otherwise parallel.

```jldoctest
julia> are_parallel([1, 1], [0, 1e-3])
false

julia> are_parallel([1, 1], [0, 1e-3], atol=1e-2)
true

julia> are_parallel([1, 1], [1, 1.01])
false

julia> are_parallel([1, 1], [1, 1.01], atol=1e-2)
true
```
"""
function are_parallel(u, v; kwargs...)

    if is_zero(u; kwargs...) || is_zero(v; kwargs...)
        return true
    end

    similarity = cosine_similarity(u, v)

    return isapprox(abs(similarity), 1; kwargs...)
end

@doc raw"""
    are_perpendicular(u, v; kwargs...) -> Bool

Check if two vectors are perpendicular.

Two vectors `u` and `v` are perpendicular if

```math
u \cdot v = 0
```

Keyword arguments can be passed along to `isapprox`.

# Examples

```jldoctest
julia> are_perpendicular([1, 0], [1, 1])
false

julia> are_perpendicular([1, 0], [0, 1])
true

julia> are_perpendicular([1, 0], [0, -5])
true

julia> are_perpendicular([1, 0], [1e-2, 1])
false

julia> are_perpendicular([1, 0], [1e-2, 1], atol=1e-2)
true
```

"""
function are_perpendicular(u, v; kwargs...)
    return isapprox(u ⋅ v, 0; kwargs...)
end

"""
    are_coplanar(points::AbstractMatrix) -> Bool

Check if multiple points are coplanar.

# Examples

```jldoctest
julia> are_coplanar([0 0; 1 1; 2 2; 3 3]')
true

julia> are_coplanar([1 5; 1 -1; 9 22; 5 7]')
true

julia> are_coplanar([4 0 6; 3 0 0; 2 0 12; 5 0 -4]')
true

julia> are_coplanar([4 1 6; 3 0 0; 2 0 12; 5 0 -4]')
false

```

"""
function are_coplanar(points::AbstractMatrix)
    point_1 = points[:, 1]
    vectors = points .- point_1

    return rank(vectors) <= 2
end

"""
    are_coplanar(line_a::AbstractLine, line_b::AbstractLine) -> Bool

Check if multiple points are coplanar.

# Examples

```jldoctest
julia> line = Line([0, 0, 0], [1, 0, 0]);

julia> are_coplanar(line, Line([0, 0, 0], [2, 5, 3]))
true

julia> are_coplanar(line, Line([1, 0, 0], [1, 1, 1]))
true

julia> are_coplanar(line, Line([0, 1, 0], [1, 1, 1]))
false
```

"""
function are_coplanar(line_a::AbstractLine, line_b::AbstractLine)

    point_a1 = line_a.point
    point_a2 = to_point(line_a)

    point_b1 = line_b.point
    point_b2 = to_point(line_b)

    points = hcat(point_a1, point_a2, point_b1, point_b2)

    return are_coplanar(points)
end

"""
    on_surface(point::AbstractVector, line::AbstractLine; kwargs...) -> Bool

Check if a point is on a line.

# Examples

```jldoctest
julia> line = Line([0, 0], [1, 0]);

julia> on_surface([0, 0], line)
true

julia> on_surface([-5, 0], line)
true

julia> on_surface([1, 1], line)
false

julia> on_surface([0, 1e-2], line)
false

julia> on_surface([0, 1e-2], line, atol=1e-2)
true
```

"""
function on_surface(point::AbstractVector, line::AbstractLine; kwargs...)
    return isapprox(distance(point, line), 0; kwargs...)
end

"""
    isapprox(line_a::AbstractLine, line_b::AbstractLine; kwargs...) -> Bool

Check if two lines are approximately equal.

# Examples

```
julia> line = Line([0, 0], [1, 0]);

julia> isapprox(line, Line([0, 0], [1, 0]))
true

julia> isapprox(line, Line([0, 0], [5, 0]))
true

julia> isapprox(line, Line([-5, 0], [-2, 0]))
true

julia> isapprox(line, Line([-5, 1], [-2, 0]))
false

julia> isapprox(line, Line([1, 1], [5, 5]))
false

julia> isapprox(line, Line([1, 1e-2], [5, 0])
false

julia> isapprox(line, Line([1, 1e-2], [5, 0]), atol=1e-2)
true
```

"""
function isapprox(line_a::AbstractLine, line_b::AbstractLine; kwargs...)

    point_on_surface = on_surface(line_a.point, line_b; kwargs...)
    directions_parallel = are_parallel(line_a.direction, line_b.direction; kwargs...)

    return point_on_surface && directions_parallel
end


text:
"""
    make_answers()
Return a list of the five-letter words that can be the solution to a Wordle puzzle.
"""
function make_answers()
    answers = [
        "aback"
        "abase"
        "abate"
        "abbey"
        "abbot"
        "abhor"
        "abide"
        "abled"
        "abode"
        "abort"
        "about"
        "above"
        "abuse"
        "abyss"
        "acorn"
        "acrid"
        "actor"
        "acute"
        "adage"
        "adapt"
        "adept"
        "admin"
        "admit"
        "adobe"
        "adopt"
        "adore"
        "adorn"
        "adult"
        "affix"
        "afire"
        "afoot"
        "afoul"
        "after"
        "again"
        "agape"
        "agate"
        "agent"
        "agile"
        "aging"
        "aglow"
        "agony"
        "agora"
        "agree"
        "ahead"
        "aider"
        "aisle"
        "alarm"
        "album"
        "alert"
        "algae"
        "alibi"
        "alien"
        "align"
        "alike"
        "alive"
        "allay"
        "alley"
        "allot"
        "allow"
        "alloy"
        "aloft"
        "alone"
        "along"
        "aloof"
        "aloud"
        "alpha"
        "altar"
        "alter"
        "amass"
        "amaze"
        "amber"
        "amble"
        "amend"
        "amiss"
        "amity"
        "among"
        "ample"
        "amply"
        "amuse"
        "angel"
        "anger"
        "angle"
        "angry"
        "angst"
        "anime"
        "ankle"
        "annex"
        "annoy"
        "annul"
        "anode"
        "antic"
        "anvil"
        "aorta"
        "apart"
        "aphid"
        "aping"
        "apnea"
        "apple"
        "apply"
        "apron"
        "aptly"
        "arbor"
        "ardor"
        "arena"
        "argue"
        "arise"
        "armor"
        "aroma"
        "arose"
        "array"
        "arrow"
        "arson"
        "artsy"
        "ascot"
        "ashen"
        "aside"
        "askew"
        "assay"
        "asset"
        "atoll"
        "atone"
        "attic"
        "audio"
        "audit"
        "augur"
        "aunty"
        "avail"
        "avert"
        "avian"
        "avoid"
        "await"
        "awake"
        "award"
        "aware"
        "awash"
        "awful"
        "awoke"
        "axial"
        "axiom"
        "axion"
        "azure"
        "bacon"
        "badge"
        "badly"
        "bagel"
        "baggy"
        "baker"
        "baler"
        "balmy"
        "banal"
        "banjo"
        "barge"
        "baron"
        "basal"
        "basic"
        "basil"
        "basin"
        "basis"
        "baste"
        "batch"
        "bathe"
        "baton"
        "batty"
        "bawdy"
        "bayou"
        "beach"
        "beady"
        "beard"
        "beast"
        "beech"
        "beefy"
        "befit"
        "began"
        "begat"
        "beget"
        "begin"
        "begun"
        "being"
        "belch"
        "belie"
        "belle"
        "belly"
        "below"
        "bench"
        "beret"
        "berry"
        "berth"
        "beset"
        "betel"
        "bevel"
        "bezel"
        "bible"
        "bicep"
        "biddy"
        "bigot"
        "bilge"
        "billy"
        "binge"
        "bingo"
        "biome"
        "birch"
        "birth"
        "bison"
        "bitty"
        "black"
        "blade"
        "blame"
        "bland"
        "blank"
        "blare"
        "blast"
        "blaze"
        "bleak"
        "bleat"
        "bleed"
        "bleep"
        "blend"
        "bless"
        "blimp"
        "blind"
        "blink"
        "bliss"
        "blitz"
        "bloat"
        "block"
        "bloke"
        "blond"
        "blood"
        "bloom"
        "blown"
        "bluer"
        "bluff"
        "blunt"
        "blurb"
        "blurt"
        "blush"
        "board"
        "boast"
        "bobby"
        "boney"
        "bongo"
        "bonus"
        "booby"
        "boost"
        "booth"
        "booty"
        "booze"
        "boozy"
        "borax"
        "borne"
        "bosom"
        "bossy"
        "botch"
        "bough"
        "boule"
        "bound"
        "bowel"
        "boxer"
        "brace"
        "braid"
        "brain"
        "brake"
        "brand"
        "brash"
        "brass"
        "brave"
        "bravo"
        "brawl"
        "brawn"
        "bread"
        "break"
        "breed"
        "briar"
        "bribe"
        "brick"
        "bride"
        "brief"
        "brine"
        "bring"
        "brink"
        "briny"
        "brisk"
        "broad"
        "broil"
        "broke"
        "brood"
        "brook"
        "broom"
        "broth"
        "brown"
        "brunt"
        "brush"
        "brute"
        "buddy"
        "budge"
        "buggy"
        "bugle"
        "build"
        "built"
        "bulge"
        "bulky"
        "bully"
        "bunch"
        "bunny"
        "burly"
        "burnt"
        "burst"
        "bused"
        "bushy"
        "butch"
        "butte"
        "buxom"
        "buyer"
        "bylaw"
        "cabal"
        "cabby"
        "cabin"
        "cable"
        "cacao"
        "cache"
        "cacti"
        "caddy"
        "cadet"
        "cagey"
        "cairn"
        "camel"
        "cameo"
        "canal"
        "candy"
        "canny"
        "canoe"
        "canon"
        "caper"
        "caput"
        "carat"
        "cargo"
        "carol"
        "carry"
        "carve"
        "caste"
        "catch"
        "cater"
        "catty"
        "caulk"
        "cause"
        "cavil"
        "cease"
        "cedar"
        "cello"
        "chafe"
        "chaff"
        "chain"
        "chair"
        "chalk"
        "champ"
        "chant"
        "chaos"
        "chard"
        "charm"
        "chart"
        "chase"
        "chasm"
        "cheap"
        "cheat"
        "check"
        "cheek"
        "cheer"
        "chess"
        "chest"
        "chick"
        "chide"
        "chief"
        "child"
        "chili"
        "chill"
        "chime"
        "china"
        "chirp"
        "chock"
        "choir"
        "choke"
        "chord"
        "chore"
        "chose"
        "chuck"
        "chump"
        "chunk"
        "churn"
        "chute"
        "cider"
        "cigar"
        "cinch"
        "circa"
        "civic"
        "civil"
        "clack"
        "claim"
        "clamp"
        "clang"
        "clank"
        "clash"
        "clasp"
        "class"
        "clean"
        "clear"
        "cleat"
        "cleft"
        "clerk"
        "click"
        "cliff"
        "climb"
        "cling"
        "clink"
        "cloak"
        "clock"
        "clone"
        "close"
        "cloth"
        "cloud"
        "clout"
        "clove"
        "clown"
        "cluck"
        "clued"
        "clump"
        "clung"
        "coach"
        "coast"
        "cobra"
        "cocoa"
        "colon"
        "color"
        "comet"
        "comfy"
        "comic"
        "comma"
        "conch"
        "condo"
        "conic"
        "copse"
        "coral"
        "corer"
        "corny"
        "couch"
        "cough"
        "could"
        "count"
        "coupe"
        "court"
        "coven"
        "cover"
        "covet"
        "covey"
        "cower"
        "coyly"
        "crack"
        "craft"
        "cramp"
        "crane"
        "crank"
        "crash"
        "crass"
        "crate"
        "crave"
        "crawl"
        "craze"
        "crazy"
        "creak"
        "cream"
        "credo"
        "creed"
        "creek"
        "creep"
        "creme"
        "crepe"
        "crept"
        "cress"
        "crest"
        "crick"
        "cried"
        "crier"
        "crime"
        "crimp"
        "crisp"
        "croak"
        "crock"
        "crone"
        "crony"
        "crook"
        "cross"
        "croup"
        "crowd"
        "crown"
        "crude"
        "cruel"
        "crumb"
        "crump"
        "crush"
        "crust"
        "crypt"
        "cubic"
        "cumin"
        "curio"
        "curly"
        "curry"
        "curse"
        "curve"
        "curvy"
        "cutie"
        "cyber"
        "cycle"
        "cynic"
        "daddy"
        "daily"
        "dairy"
        "daisy"
        "dally"
        "dance"
        "dandy"
        "datum"
        "daunt"
        "dealt"
        "death"
        "debar"
        "debit"
        "debug"
        "debut"
        "decal"
        "decay"
        "decor"
        "decoy"
        "decry"
        "defer"
        "deign"
        "deity"
        "delay"
        "delta"
        "delve"
        "demon"
        "demur"
        "denim"
        "dense"
        "depot"
        "depth"
        "derby"
        "deter"
        "detox"
        "deuce"
        "devil"
        "diary"
        "dicey"
        "digit"
        "dilly"
        "dimly"
        "diner"
        "dingo"
        "dingy"
        "diode"
        "dirge"
        "dirty"
        "disco"
        "ditch"
        "ditto"
        "ditty"
        "diver"
        "dizzy"
        "dodge"
        "dodgy"
        "dogma"
        "doing"
        "dolly"
        "donor"
        "donut"
        "dopey"
        "doubt"
        "dough"
        "dowdy"
        "dowel"
        "downy"
        "dowry"
        "dozen"
        "draft"
        "drain"
        "drake"
        "drama"
        "drank"
        "drape"
        "drawl"
        "drawn"
        "dread"
        "dream"
        "dress"
        "dried"
        "drier"
        "drift"
        "drill"
        "drink"
        "drive"
        "droit"
        "droll"
        "drone"
        "drool"
        "droop"
        "dross"
        "drove"
        "drown"
        "druid"
        "drunk"
        "dryer"
        "dryly"
        "duchy"
        "dully"
        "dummy"
        "dumpy"
        "dunce"
        "dusky"
        "dusty"
        "dutch"
        "duvet"
        "dwarf"
        "dwell"
        "dwelt"
        "dying"
        "eager"
        "eagle"
        "early"
        "earth"
        "easel"
        "eaten"
        "eater"
        "ebony"
        "eclat"
        "edict"
        "edify"
        "eerie"
        "egret"
        "eight"
        "eject"
        "eking"
        "elate"
        "elbow"
        "elder"
        "elect"
        "elegy"
        "elfin"
        "elide"
        "elite"
        "elope"
        "elude"
        "email"
        "embed"
        "ember"
        "emcee"
        "empty"
        "enact"
        "endow"
        "enema"
        "enemy"
        "enjoy"
        "ennui"
        "ensue"
        "enter"
        "entry"
        "envoy"
        "epoch"
        "epoxy"
        "equal"
        "equip"
        "erase"
        "erect"
        "erode"
        "error"
        "erupt"
        "essay"
        "ester"
        "ether"
        "ethic"
        "ethos"
        "etude"
        "evade"
        "event"
        "every"
        "evict"
        "evoke"
        "exact"
        "exalt"
        "excel"
        "exert"
        "exile"
        "exist"
        "expel"
        "extol"
        "extra"
        "exult"
        "eying"
        "fable"
        "facet"
        "faint"
        "fairy"
        "faith"
        "false"
        "fancy"
        "fanny"
        "farce"
        "fatal"
        "fatty"
        "fault"
        "fauna"
        "favor"
        "feast"
        "fecal"
        "feign"
        "fella"
        "felon"
        "femme"
        "femur"
        "fence"
        "feral"
        "ferry"
        "fetal"
        "fetch"
        "fetid"
        "fetus"
        "fever"
        "fewer"
        "fiber"
        "fibre"
        "ficus"
        "field"
        "fiend"
        "fiery"
        "fifth"
        "fifty"
        "fight"
        "filer"
        "filet"
        "filly"
        "filmy"
        "filth"
        "final"
        "finch"
        "finer"
        "first"
        "fishy"
        "fixer"
        "fizzy"
        "fjord"
        "flack"
        "flail"
        "flair"
        "flake"
        "flaky"
        "flame"
        "flank"
        "flare"
        "flash"
        "flask"
        "fleck"
        "fleet"
        "flesh"
        "flick"
        "flier"
        "fling"
        "flint"
        "flirt"
        "float"
        "flock"
        "flood"
        "floor"
        "flora"
        "floss"
        "flour"
        "flout"
        "flown"
        "fluff"
        "fluid"
        "fluke"
        "flume"
        "flung"
        "flunk"
        "flush"
        "flute"
        "flyer"
        "foamy"
        "focal"
        "focus"
        "foggy"
        "foist"
        "folio"
        "folly"
        "foray"
        "force"
        "forge"
        "forgo"
        "forte"
        "forth"
        "forty"
        "forum"
        "found"
        "foyer"
        "frail"
        "frame"
        "frank"
        "fraud"
        "freak"
        "freed"
        "freer"
        "fresh"
        "friar"
        "fried"
        "frill"
        "frisk"
        "fritz"
        "frock"
        "frond"
        "front"
        "frost"
        "froth"
        "frown"
        "froze"
        "fruit"
        "fudge"
        "fugue"
        "fully"
        "fungi"
        "funky"
        "funny"
        "furor"
        "furry"
        "fussy"
        "fuzzy"
        "gaffe"
        "gaily"
        "gamer"
        "gamma"
        "gamut"
        "gassy"
        "gaudy"
        "gauge"
        "gaunt"
        "gauze"
        "gavel"
        "gawky"
        "gayer"
        "gayly"
        "gazer"
        "gecko"
        "geeky"
        "geese"
        "genie"
        "genre"
        "ghost"
        "ghoul"
        "giant"
        "giddy"
        "gipsy"
        "girly"
        "girth"
        "given"
        "giver"
        "glade"
        "gland"
        "glare"
        "glass"
        "glaze"
        "gleam"
        "glean"
        "glide"
        "glint"
        "gloat"
        "globe"
        "gloom"
        "glory"
        "gloss"
        "glove"
        "glyph"
        "gnash"
        "gnome"
        "godly"
        "going"
        "golem"
        "golly"
        "gonad"
        "goner"
        "goody"
        "gooey"
        "goofy"
        "goose"
        "gorge"
        "gouge"
        "gourd"
        "grace"
        "grade"
        "graft"
        "grail"
        "grain"
        "grand"
        "grant"
        "grape"
        "graph"
        "grasp"
        "grass"
        "grate"
        "grave"
        "gravy"
        "graze"
        "great"
        "greed"
        "green"
        "greet"
        "grief"
        "grill"
        "grime"
        "grimy"
        "grind"
        "gripe"
        "groan"
        "groin"
        "groom"
        "grope"
        "gross"
        "group"
        "grout"
        "grove"
        "growl"
        "grown"
        "gruel"
        "gruff"
        "grunt"
        "guard"
        "guava"
        "guess"
        "guest"
        "guide"
        "guild"
        "guile"
        "guilt"
        "guise"
        "gulch"
        "gully"
        "gumbo"
        "gummy"
        "guppy"
        "gusto"
        "gusty"
        "gypsy"
        "habit"
        "hairy"
        "halve"
        "handy"
        "happy"
        "hardy"
        "harem"
        "harpy"
        "harry"
        "harsh"
        "haste"
        "hasty"
        "hatch"
        "hater"
        "haunt"
        "haute"
        "haven"
        "havoc"
        "hazel"
        "heady"
        "heard"
        "heart"
        "heath"
        "heave"
        "heavy"
        "hedge"
        "hefty"
        "heist"
        "helix"
        "hello"
        "hence"
        "heron"
        "hilly"
        "hinge"
        "hippo"
        "hippy"
        "hitch"
        "hoard"
        "hobby"
        "hoist"
        "holly"
        "homer"
        "honey"
        "honor"
        "horde"
        "horny"
        "horse"
        "hotel"
        "hotly"
        "hound"
        "house"
        "hovel"
        "hover"
        "howdy"
        "human"
        "humid"
        "humor"
        "humph"
        "humus"
        "hunch"
        "hunky"
        "hurry"
        "husky"
        "hussy"
        "hutch"
        "hydro"
        "hyena"
        "hymen"
        "hyper"
        "icily"
        "icing"
        "ideal"
        "idiom"
        "idiot"
        "idler"
        "idyll"
        "igloo"
        "iliac"
        "image"
        "imbue"
        "impel"
        "imply"
        "inane"
        "inbox"
        "incur"
        "index"
        "inept"
        "inert"
        "infer"
        "ingot"
        "inlay"
        "inlet"
        "inner"
        "input"
        "inter"
        "intro"
        "ionic"
        "irate"
        "irony"
        "islet"
        "issue"
        "itchy"
        "ivory"
        "jaunt"
        "jazzy"
        "jelly"
        "jerky"
        "jetty"
        "jewel"
        "jiffy"
        "joint"
        "joist"
        "joker"
        "jolly"
        "joust"
        "judge"
        "juice"
        "juicy"
        "jumbo"
        "jumpy"
        "junta"
        "junto"
        "juror"
        "kappa"
        "karma"
        "kayak"
        "kebab"
        "khaki"
        "kinky"
        "kiosk"
        "kitty"
        "knack"
        "knave"
        "knead"
        "kneed"
        "kneel"
        "knelt"
        "knife"
        "knock"
        "knoll"
        "known"
        "koala"
        "krill"
        "label"
        "labor"
        "laden"
        "ladle"
        "lager"
        "lance"
        "lanky"
        "lapel"
        "lapse"
        "large"
        "larva"
        "lasso"
        "latch"
        "later"
        "lathe"
        "latte"
        "laugh"
        "layer"
        "leach"
        "leafy"
        "leaky"
        "leant"
        "leapt"
        "learn"
        "lease"
        "leash"
        "least"
        "leave"
        "ledge"
        "leech"
        "leery"
        "lefty"
        "legal"
        "leggy"
        "lemon"
        "lemur"
        "leper"
        "level"
        "lever"
        "libel"
        "liege"
        "light"
        "liken"
        "lilac"
        "limbo"
        "limit"
        "linen"
        "liner"
        "lingo"
        "lipid"
        "lithe"
        "liver"
        "livid"
        "llama"
        "loamy"
        "loath"
        "lobby"
        "local"
        "locus"
        "lodge"
        "lofty"
        "logic"
        "login"
        "loopy"
        "loose"
        "lorry"
        "loser"
        "louse"
        "lousy"
        "lover"
        "lower"
        "lowly"
        "loyal"
        "lucid"
        "lucky"
        "lumen"
        "lumpy"
        "lunar"
        "lunch"
        "lunge"
        "lupus"
        "lurch"
        "lurid"
        "lusty"
        "lying"
        "lymph"
        "lynch"
        "lyric"
        "macaw"
        "macho"
        "macro"
        "madam"
        "madly"
        "mafia"
        "magic"
        "magma"
        "maize"
        "major"
        "maker"
        "mambo"
        "mamma"
        "mammy"
        "manga"
        "mange"
        "mango"
        "mangy"
        "mania"
        "manic"
        "manly"
        "manor"
        "maple"
        "march"
        "marry"
        "marsh"
        "mason"
        "masse"
        "match"
        "matey"
        "mauve"
        "maxim"
        "maybe"
        "mayor"
        "mealy"
        "meant"
        "meaty"
        "mecca"
        "medal"
        "media"
        "medic"
        "melee"
        "melon"
        "mercy"
        "merge"
        "merit"
        "merry"
        "metal"
        "meter"
        "metro"
        "micro"
        "midge"
        "midst"
        "might"
        "milky"
        "mimic"
        "mince"
        "miner"
        "minim"
        "minor"
        "minty"
        "minus"
        "mirth"
        "miser"
        "missy"
        "mocha"
        "modal"
        "model"
        "modem"
        "mogul"
        "moist"
        "molar"
        "moldy"
        "money"
        "month"
        "moody"
        "moose"
        "moral"
        "moron"
        "morph"
        "mossy"
        "motel"
        "motif"
        "motor"
        "motto"
        "moult"
        "mound"
        "mount"
        "mourn"
        "mouse"
        "mouth"
        "mover"
        "movie"
        "mower"
        "mucky"
        "mucus"
        "muddy"
        "mulch"
        "mummy"
        "munch"
        "mural"
        "murky"
        "mushy"
        "music"
        "musky"
        "musty"
        "myrrh"
        "nadir"
        "naive"
        "nanny"
        "nasal"
        "nasty"
        "natal"
        "naval"
        "navel"
        "needy"
        "neigh"
        "nerdy"
        "nerve"
        "never"
        "newer"
        "newly"
        "nicer"
        "niche"
        "niece"
        "night"
        "ninja"
        "ninny"
        "ninth"
        "noble"
        "nobly"
        "noise"
        "noisy"
        "nomad"
        "noose"
        "north"
        "nosey"
        "notch"
        "novel"
        "nudge"
        "nurse"
        "nutty"
        "nylon"
        "nymph"
        "oaken"
        "obese"
        "occur"
        "ocean"
        "octal"
        "octet"
        "odder"
        "oddly"
        "offal"
        "offer"
        "often"
        "olden"
        "older"
        "olive"
        "ombre"
        "omega"
        "onion"
        "onset"
        "opera"
        "opine"
        "opium"
        "optic"
        "orbit"
        "order"
        "organ"
        "other"
        "otter"
        "ought"
        "ounce"
        "outdo"
        "outer"
        "outgo"
        "ovary"
        "ovate"
        "overt"
        "ovine"
        "ovoid"
        "owing"
        "owner"
        "oxide"
        "ozone"
        "paddy"
        "pagan"
        "paint"
        "paler"
        "palsy"
        "panel"
        "panic"
        "pansy"
        "papal"
        "paper"
        "parer"
        "parka"
        "parry"
        "parse"
        "party"
        "pasta"
        "paste"
        "pasty"
        "patch"
        "patio"
        "patsy"
        "patty"
        "pause"
        "payee"
        "payer"
        "peace"
        "peach"
        "pearl"
        "pecan"
        "pedal"
        "penal"
        "pence"
        "penne"
        "penny"
        "perch"
        "peril"
        "perky"
        "pesky"
        "pesto"
        "petal"
        "petty"
        "phase"
        "phone"
        "phony"
        "photo"
        "piano"
        "picky"
        "piece"
        "piety"
        "piggy"
        "pilot"
        "pinch"
        "piney"
        "pinky"
        "pinto"
        "piper"
        "pique"
        "pitch"
        "pithy"
        "pivot"
        "pixel"
        "pixie"
        "pizza"
        "place"
        "plaid"
        "plain"
        "plait"
        "plane"
        "plank"
        "plant"
        "plate"
        "plaza"
        "plead"
        "pleat"
        "plied"
        "plier"
        "pluck"
        "plumb"
        "plume"
        "plump"
        "plunk"
        "plush"
        "poesy"
        "point"
        "poise"
        "poker"
        "polar"
        "polka"
        "polyp"
        "pooch"
        "poppy"
        "porch"
        "poser"
        "posit"
        "posse"
        "pouch"
        "pound"
        "pouty"
        "power"
        "prank"
        "prawn"
        "preen"
        "press"
        "price"
        "prick"
        "pride"
        "pried"
        "prime"
        "primo"
        "print"
        "prior"
        "prism"
        "privy"
        "prize"
        "probe"
        "prone"
        "prong"
        "proof"
        "prose"
        "proud"
        "prove"
        "prowl"
        "proxy"
        "prude"
        "prune"
        "psalm"
        "pubic"
        "pudgy"
        "puffy"
        "pulpy"
        "pulse"
        "punch"
        "pupal"
        "pupil"
        "puppy"
        "puree"
        "purer"
        "purge"
        "purse"
        "pushy"
        "putty"
        "pygmy"
        "quack"
        "quail"
        "quake"
        "qualm"
        "quark"
        "quart"
        "quash"
        "quasi"
        "queen"
        "queer"
        "quell"
        "query"
        "quest"
        "queue"
        "quick"
        "quiet"
        "quill"
        "quilt"
        "quirk"
        "quite"
        "quota"
        "quote"
        "quoth"
        "rabbi"
        "rabid"
        "racer"
        "radar"
        "radii"
        "radio"
        "rainy"
        "raise"
        "rajah"
        "rally"
        "ralph"
        "ramen"
        "ranch"
        "randy"
        "range"
        "rapid"
        "rarer"
        "raspy"
        "ratio"
        "ratty"
        "raven"
        "rayon"
        "razor"
        "reach"
        "react"
        "ready"
        "realm"
        "rearm"
        "rebar"
        "rebel"
        "rebus"
        "rebut"
        "recap"
        "recur"
        "recut"
        "reedy"
        "refer"
        "refit"
        "regal"
        "rehab"
        "reign"
        "relax"
        "relay"
        "relic"
        "remit"
        "renal"
        "renew"
        "repay"
        "repel"
        "reply"
        "rerun"
        "reset"
        "resin"
        "retch"
        "retro"
        "retry"
        "reuse"
        "revel"
        "revue"
        "rhino"
        "rhyme"
        "rider"
        "ridge"
        "rifle"
        "right"
        "rigid"
        "rigor"
        "rinse"
        "ripen"
        "riper"
        "risen"
        "riser"
        "risky"
        "rival"
        "river"
        "rivet"
        "roach"
        "roast"
        "robin"
        "robot"
        "rocky"
        "rodeo"
        "roger"
        "rogue"
        "roomy"
        "roost"
        "rotor"
        "rouge"
        "rough"
        "round"
        "rouse"
        "route"
        "rover"
        "rowdy"
        "rower"
        "royal"
        "ruddy"
        "ruder"
        "rugby"
        "ruler"
        "rumba"
        "rumor"
        "rupee"
        "rural"
        "rusty"
        "sadly"
        "safer"
        "saint"
        "salad"
        "sally"
        "salon"
        "salsa"
        "salty"
        "salve"
        "salvo"
        "sandy"
        "saner"
        "sappy"
        "sassy"
        "satin"
        "satyr"
        "sauce"
        "saucy"
        "sauna"
        "saute"
        "savor"
        "savoy"
        "savvy"
        "scald"
        "scale"
        "scalp"
        "scaly"
        "scamp"
        "scant"
        "scare"
        "scarf"
        "scary"
        "scene"
        "scent"
        "scion"
        "scoff"
        "scold"
        "scone"
        "scoop"
        "scope"
        "score"
        "scorn"
        "scour"
        "scout"
        "scowl"
        "scram"
        "scrap"
        "scree"
        "screw"
        "scrub"
        "scrum"
        "scuba"
        "sedan"
        "seedy"
        "segue"
        "seize"
        "semen"
        "sense"
        "sepia"
        "serif"
        "serum"
        "serve"
        "setup"
        "seven"
        "sever"
        "sewer"
        "shack"
        "shade"
        "shady"
        "shaft"
        "shake"
        "shaky"
        "shale"
        "shall"
        "shalt"
        "shame"
        "shank"
        "shape"
        "shard"
        "share"
        "shark"
        "sharp"
        "shave"
        "shawl"
        "shear"
        "sheen"
        "sheep"
        "sheer"
        "sheet"
        "sheik"
        "shelf"
        "shell"
        "shied"
        "shift"
        "shine"
        "shiny"
        "shire"
        "shirk"
        "shirt"
        "shoal"
        "shock"
        "shone"
        "shook"
        "shoot"
        "shore"
        "shorn"
        "short"
        "shout"
        "shove"
        "shown"
        "showy"
        "shrew"
        "shrub"
        "shrug"
        "shuck"
        "shunt"
        "shush"
        "shyly"
        "siege"
        "sieve"
        "sight"
        "sigma"
        "silky"
        "silly"
        "since"
        "sinew"
        "singe"
        "siren"
        "sissy"
        "sixth"
        "sixty"
        "skate"
        "skier"
        "skiff"
        "skill"
        "skimp"
        "skirt"
        "skulk"
        "skull"
        "skunk"
        "slack"
        "slain"
        "slang"
        "slant"
        "slash"
        "slate"
        "slave"
        "sleek"
        "sleep"
        "sleet"
        "slept"
        "slice"
        "slick"
        "slide"
        "slime"
        "slimy"
        "sling"
        "slink"
        "sloop"
        "slope"
        "slosh"
        "sloth"
        "slump"
        "slung"
        "slunk"
        "slurp"
        "slush"
        "slyly"
        "smack"
        "small"
        "smart"
        "smash"
        "smear"
        "smell"
        "smelt"
        "smile"
        "smirk"
        "smite"
        "smith"
        "smock"
        "smoke"
        "smoky"
        "smote"
        "snack"
        "snail"
        "snake"
        "snaky"
        "snare"
        "snarl"
        "sneak"
        "sneer"
        "snide"
        "sniff"
        "snipe"
        "snoop"
        "snore"
        "snort"
        "snout"
        "snowy"
        "snuck"
        "snuff"
        "soapy"
        "sober"
        "soggy"
        "solar"
        "solid"
        "solve"
        "sonar"
        "sonic"
        "sooth"
        "sooty"
        "sorry"
        "sound"
        "south"
        "sower"
        "space"
        "spade"
        "spank"
        "spare"
        "spark"
        "spasm"
        "spawn"
        "speak"
        "spear"
        "speck"
        "speed"
        "spell"
        "spelt"
        "spend"
        "spent"
        "sperm"
        "spice"
        "spicy"
        "spied"
        "spiel"
        "spike"
        "spiky"
        "spill"
        "spilt"
        "spine"
        "spiny"
        "spire"
        "spite"
        "splat"
        "split"
        "spoil"
        "spoke"
        "spoof"
        "spook"
        "spool"
        "spoon"
        "spore"
        "sport"
        "spout"
        "spray"
        "spree"
        "sprig"
        "spunk"
        "spurn"
        "spurt"
        "squad"
        "squat"
        "squib"
        "stack"
        "staff"
        "stage"
        "staid"
        "stain"
        "stair"
        "stake"
        "stale"
        "stalk"
        "stall"
        "stamp"
        "stand"
        "stank"
        "stare"
        "stark"
        "start"
        "stash"
        "state"
        "stave"
        "stead"
        "steak"
        "steal"
        "steam"
        "steed"
        "steel"
        "steep"
        "steer"
        "stein"
        "stern"
        "stick"
        "stiff"
        "still"
        "stilt"
        "sting"
        "stink"
        "stint"
        "stock"
        "stoic"
        "stoke"
        "stole"
        "stomp"
        "stone"
        "stony"
        "stood"
        "stool"
        "stoop"
        "store"
        "stork"
        "storm"
        "story"
        "stout"
        "stove"
        "strap"
        "straw"
        "stray"
        "strip"
        "strut"
        "stuck"
        "study"
        "stuff"
        "stump"
        "stung"
        "stunk"
        "stunt"
        "style"
        "suave"
        "sugar"
        "suing"
        "suite"
        "sulky"
        "sully"
        "sumac"
        "sunny"
        "super"
        "surer"
        "surge"
        "surly"
        "sushi"
        "swami"
        "swamp"
        "swarm"
        "swash"
        "swath"
        "swear"
        "sweat"
        "sweep"
        "sweet"
        "swell"
        "swept"
        "swift"
        "swill"
        "swine"
        "swing"
        "swirl"
        "swish"
        "swoon"
        "swoop"
        "sword"
        "swore"
        "sworn"
        "swung"
        "synod"
        "syrup"
        "tabby"
        "table"
        "taboo"
        "tacit"
        "tacky"
        "taffy"
        "taint"
        "taken"
        "taker"
        "tally"
        "talon"
        "tamer"
        "tango"
        "tangy"
        "taper"
        "tapir"
        "tardy"
        "tarot"
        "taste"
        "tasty"
        "tatty"
        "taunt"
        "tawny"
        "teach"
        "teary"
        "tease"
        "teddy"
        "teeth"
        "tempo"
        "tenet"
        "tenor"
        "tense"
        "tenth"
        "tepee"
        "tepid"
        "terra"
        "terse"
        "testy"
        "thank"
        "theft"
        "their"
        "theme"
        "there"
        "these"
        "theta"
        "thick"
        "thief"
        "thigh"
        "thing"
        "think"
        "third"
        "thong"
        "thorn"
        "those"
        "three"
        "threw"
        "throb"
        "throw"
        "thrum"
        "thumb"
        "thump"
        "thyme"
        "tiara"
        "tibia"
        "tidal"
        "tiger"
        "tight"
        "tilde"
        "timer"
        "timid"
        "tipsy"
        "titan"
        "tithe"
        "title"
        "toast"
        "today"
        "toddy"
        "token"
        "tonal"
        "tonga"
        "tonic"
        "tooth"
        "topaz"
        "topic"
        "torch"
        "torso"
        "torus"
        "total"
        "totem"
        "touch"
        "tough"
        "towel"
        "tower"
        "toxic"
        "toxin"
        "trace"
        "track"
        "tract"
        "trade"
        "trail"
        "train"
        "trait"
        "tramp"
        "trash"
        "trawl"
        "tread"
        "treat"
        "trend"
        "triad"
        "trial"
        "tribe"
        "trice"
        "trick"
        "tried"
        "tripe"
        "trite"
        "troll"
        "troop"
        "trope"
        "trout"
        "trove"
        "truce"
        "truck"
        "truer"
        "truly"
        "trump"
        "trunk"
        "truss"
        "trust"
        "truth"
        "tryst"
        "tubal"
        "tuber"
        "tulip"
        "tulle"
        "tumor"
        "tunic"
        "turbo"
        "tutor"
        "twang"
        "tweak"
        "tweed"
        "tweet"
        "twice"
        "twine"
        "twirl"
        "twist"
        "twixt"
        "tying"
        "udder"
        "ulcer"
        "ultra"
        "umbra"
        "uncle"
        "uncut"
        "under"
        "undid"
        "undue"
        "unfed"
        "unfit"
        "unify"
        "union"
        "unite"
        "unity"
        "unlit"
        "unmet"
        "unset"
        "untie"
        "until"
        "unwed"
        "unzip"
        "upper"
        "upset"
        "urban"
        "urine"
        "usage"
        "usher"
        "using"
        "usual"
        "usurp"
        "utile"
        "utter"
        "vague"
        "valet"
        "valid"
        "valor"
        "value"
        "valve"
        "vapid"
        "vapor"
        "vault"
        "vaunt"
        "vegan"
        "venom"
        "venue"
        "verge"
        "verse"
        "verso"
        "verve"
        "vicar"
        "video"
        "vigil"
        "vigor"
        "villa"
        "vinyl"
        "viola"
        "viper"
        "viral"
        "virus"
        "visit"
        "visor"
        "vista"
        "vital"
        "vivid"
        "vixen"
        "vocal"
        "vodka"
        "vogue"
        "voice"
        "voila"
        "vomit"
        "voter"
        "vouch"
        "vowel"
        "vying"
        "wacky"
        "wafer"
        "wager"
        "wagon"
        "waist"
        "waive"
        "waltz"
        "warty"
        "waste"
        "watch"
        "water"
        "waver"
        "waxen"
        "weary"
        "weave"
        "wedge"
        "weedy"
        "weigh"
        "weird"
        "welch"
        "welsh"
        "wench"
        "whack"
        "whale"
        "wharf"
        "wheat"
        "wheel"
        "whelp"
        "where"
        "which"
        "whiff"
        "while"
        "whine"
        "whiny"
        "whirl"
        "whisk"
        "white"
        "whole"
        "whoop"
        "whose"
        "widen"
        "wider"
        "widow"
        "width"
        "wield"
        "wight"
        "willy"
        "wimpy"
        "wince"
        "winch"
        "windy"
        "wiser"
        "wispy"
        "witch"
        "witty"
        "woken"
        "woman"
        "women"
        "woody"
        "wooer"
        "wooly"
        "woozy"
        "wordy"
        "world"
        "worry"
        "worse"
        "worst"
        "worth"
        "would"
        "wound"
        "woven"
        "wrack"
        "wrath"
        "wreak"
        "wreck"
        "wrest"
        "wring"
        "wrist"
        "write"
        "wrong"
        "wrote"
        "wrung"
        "wryly"
        "yacht"
        "yearn"
        "yeast"
        "yield"
        "young"
        "youth"
        "zebra"
        "zesty"
        "zonal"
    ]

    answers = uppercase.(answers)
    return answers
end


text:
mutable struct timestepper
    Gcs::NamedTuple     # a NamedTuple same as nutrients to store tendencies
    nut_temp::NamedTuple# a NamedTuple same as nutrients to store nutrients fields in multi-dims advection scheme
    vel₀::NamedTuple    # a NamedTuple with u, v, w velocities
    vel½::NamedTuple    # a NamedTuple with u, v, w velocities
    vel₁::NamedTuple    # a NamedTuple with u, v, w velocities
    PARF::AbstractArray # a (Cu)Array to store surface PAR field of each timestep
    temp::AbstractArray # a (Cu)Array to store temperature field of each timestep
    plk::NamedTuple     # a NamedTuple same as nutrients to store interactions with individuals
    par::AbstractArray  # a (Cu)Array to store PAR field of each timestep
    Chl::AbstractArray  # a (Cu)Array to store Chl field of each timestep
    pop::AbstractArray  # a (Cu)Array to store population field of each timestep
    rnd::AbstractArray  # a StructArray of random numbers for plankton diffusion or grazing, mortality and division.
    velos::AbstractArray# a StructArray of intermediate values for RK4 particle advection
    nuts::AbstractArray # a StructArray of nutrients of each individual
end

function timestepper(arch::Architecture, g::AbstractGrid, maxN)
    vel₀ = (u = Field(arch, g), v = Field(arch, g), w = Field(arch, g))
    vel½ = (u = Field(arch, g), v = Field(arch, g), w = Field(arch, g))
    vel₁ = (u = Field(arch, g), v = Field(arch, g), w = Field(arch, g))

    Gcs = nutrients_init(arch, g)
    nut_temp = nutrients_init(arch, g)
    plk = nutrients_init(arch, g)

    par = zeros(g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2) |> array_type(arch)
    Chl = zeros(g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2) |> array_type(arch)
    pop = zeros(g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2) |> array_type(arch)

    temp = zeros(g.Nx+g.Hx*2, g.Ny+g.Hy*2, g.Nz+g.Hz*2) |> array_type(arch)
    PARF = zeros(g.Nx, g.Ny) |> array_type(arch)

    rnd = StructArray(x = zeros(maxN), y = zeros(maxN), z = zeros(maxN))
    rnd_d = replace_storage(array_type(arch), rnd)

    velos = StructArray(x  = zeros(maxN), y  = zeros(maxN), z  = zeros(maxN),
                        u1 = zeros(maxN), v1 = zeros(maxN), w1 = zeros(maxN),
                        u2 = zeros(maxN), v2 = zeros(maxN), w2 = zeros(maxN),
                        )
    velos_d = replace_storage(array_type(arch), velos)

    nuts = StructArray(NH4 = zeros(maxN), NO3 = zeros(maxN), PO4 = zeros(maxN), DOC = zeros(maxN),
                       par = zeros(maxN), T   = zeros(maxN), pop = zeros(maxN))
    nuts_d = replace_storage(array_type(arch), nuts)

    ts = timestepper(Gcs, nut_temp, vel₀, vel½, vel₁, PARF, temp, plk, par, Chl, pop, rnd_d, velos_d, nuts_d)

    return ts
end


text:
using Aurora
using Test

@testset "Replicated Samples" begin
    include("test_replicated_sample.jl")
end

@testset "Auroral" begin
    include("test_auroral.jl")
end

@testset "Aurora-KNN" begin
    include("test_aurora_nn.jl")
end


text:
module CellFields

using Test
using Gridap
using Gridap.Helpers
using Gridap.CellValues: _eq

export IterCellFieldLike
export IndexCellFieldLike
export CellFieldLike

export IterCellField
export IndexCellField
export CellField

export IterCellBasis
export IndexCellBasis
export CellBasis

export CellGeomap

export CellFieldValues
export CellBasisValues
export CellPoints

export test_iter_cell_field
export test_index_cell_field
export test_iter_cell_basis
export test_index_cell_basis
export test_iter_cell_field_without_grad
export test_index_cell_field_without_grad
export test_iter_cell_basis_without_grad
export test_index_cell_basis_without_grad

import Gridap: gradient

const IterCellFieldLike{D,T<:FieldValue,N,R<:FieldLike{D,T,N}} = IterCellValue{R}

const IndexCellFieldLike{D,T<:FieldValue,N,C,R<:FieldLike{D,T,N}} = IndexCellValue{R,C}

const NonIterableCellFieldLike{D,T<:FieldValue,N} = NonIterableCellMap{Point{D},1,T,N}

"""
Umbrella type for CellField and CellBasis
"""
const CellFieldLike{D,T<:FieldValue,N} = Union{IterCellFieldLike{D,T,N},IndexCellFieldLike{D,T,N}}

"""
Returns another CellField or CellBasis object representing the gradient of the given one.
"""
function gradient(::CellFieldLike{D,T,N})::CellFieldLike{D,G,N} where {D,T,G,N}
  @abstractmethod
end

const IterCellField{D,T<:FieldValue,R<:Field{D,T}} = IterCellValue{R}

const IndexCellField{D,T<:FieldValue,C,R<:Field{D,T}} = IndexCellValue{R,C}

"""
Abstract type that represents a cell-wise field, where
`T` stands for the type that represents the field at a point
(e.g., scalar, vector, tensor) and `D` stands for the space
dimension
"""
const CellField{D,T<:FieldValue} = Union{IterCellField{D,T},IndexCellField{D,T}}

function CellField end

const IterCellBasis{D,T<:FieldValue,R<:Basis{D,T}} = IterCellValue{R}

const IndexCellBasis{D,T<:FieldValue,C,R<:Basis{D,T}} = IndexCellValue{R,C}

"""
Abstract type that represents a cell-wise basis for a field space,
where T is the type of value and D the dimension of the domain
"""
const CellBasis{D,T<:FieldValue} = Union{IterCellBasis{D,T},IndexCellBasis{D,T}}

function CellBasis end

"""
Abstract type representing a cellwise transformation between two geometrical
domains
"""
const CellGeomap{D,Z,X,T<:Point{Z,X}} = CellField{D,T}

function CellGeomap end

"""
An array of points for each cell.
This type represents the objects where CellField and CellBasis are evaluated
"""
const CellPoints{D,X} = CellVector{Point{D,X}}

"""
Abstract type that represents a field with value of type T
evaluated at a collection of points in each cell
"""
const CellFieldValues{T<:FieldValue} = CellVector{T}

"""
Abstract type that represents a function basis with value of type T
evaluated at a collection of points in each cell
"""
const CellBasisValues{T<:FieldValue} = CellMatrix{T}

# Testers

function test_iter_cell_field_like(
  f::IterCellFieldLike{D,T,N},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractArray{T,N}},
  g::AbstractArray{<:AbstractArray{G,N}}) where {D,T,G,N}

  _test_iter_cell_field_like(f,x,v,g)

end

function _test_iter_cell_field_like(f,x,v,g)
  test_iter_cell_map(f,x,v)
  fg = gradient(f)
  test_iter_cell_map(fg,x,g)
  _test_field_like_iteration(f,x,v,g)
end

function test_index_cell_field_like(
  f::IndexCellFieldLike{D,T,N},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractArray{T,N}},
  g::AbstractArray{<:AbstractArray{G,N}}) where {D,T,G,N}

  test_index_cell_map(f,x,v)
  fg = gradient(f)
  test_index_cell_map(fg,x,g)
  _test_field_like_iteration(f,x,v,g)
end

function _test_field_like_iteration(m,x,v,g)
  for (mi,xi,vi,gi) in zip(m,x,v,g)
  #  @assert _eq(evaluate(mi,xi),vi)
    @assert _eq(evaluate(gradient(mi),xi),gi)
    @assert typeof(mi) == eltype(m)
  end
end

function test_iter_cell_field(
  f::IterCellField{D,T},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractVector{T}},
  g::AbstractArray{<:AbstractVector{G}}) where {D,T,G}

  _test_iter_cell_field_like(f,x,v,g)
end

function test_iter_cell_field_without_grad(
  f::IterCellField{D,T},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractVector{T}}) where {D,T}

  test_iter_cell_map(f,x,v)
end

function test_index_cell_field(
  f::IndexCellField{D,T},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractVector{T}},
  g::AbstractArray{<:AbstractVector{G}}) where {D,T,G}

  test_index_cell_field_like(f,x,v,g)
end

function test_index_cell_field_without_grad(
  f::IndexCellField{D,T},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractVector{T}}) where {D,T}

  test_index_cell_map(f,x,v)
end

function test_iter_cell_basis(
  f::IterCellBasis{D,T},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractMatrix{T}},
  g::AbstractArray{<:AbstractMatrix{G}}) where {D,T,G}

  _test_iter_cell_field_like(f,x,v,g)
end

function test_iter_cell_basis_without_grad(
  f::IterCellBasis{D,T},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractMatrix{T}}) where {D,T}

  test_iter_cell_map(f,x,v)
end

function test_index_cell_basis(
  f::IndexCellBasis{D,T},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractMatrix{T}},
  g::AbstractArray{<:AbstractMatrix{G}}) where {D,T,G}

  test_index_cell_field_like(f,x,v,g)
end

function test_index_cell_basis_without_grad(
  f::IndexCellBasis{D,T},
  x::CellPoints{D},
  v::AbstractArray{<:AbstractMatrix{T}}) where {D,T}

  test_index_cell_map(f,x,v)
end

end # module


text:
export AbstractGraph, Graph
export DiAdjacencyMatrix, WeightedDiAdjacencyMatrix
export DiAdjacencyList, WeightedAdjacencyList
export nv, ne, relate!, unrelate!, neighbor, weight, degree
export probability, random_walk

abstract type AbstractGraph end

mutable struct DiAdjacencyMatrix <: AbstractGraph
    relation::BitArray
    n_vertices::Int64

    DiAdjacencyMatrix(n::Int64) = new(zeros(Bool, n, n), n)
end

mutable struct WeightedDiAdjacencyMatrix{T <: Number} <: AbstractGraph
    relation::Array{T}
    n_vertices::Int64

    function WeightedDiAdjacencyMatrix{T}(n::Int64; random::Bool=true) where {T <: Number}
        return random ? new(abs.(rand(T, n, n)), n) : new(zeros(T, n, n), n)
    end
end

nv(g::AbstractGraph) = g.n_vertices

ne(g::AbstractGraph) = sum(g.relation .!= 0)

function neighbor(g::AbstractGraph, v::Int64)
    neighbor_list = Int[]
    for (i, v) in enumerate(g.relation[v, 1:end])
        v == 0 && continue
        push!(neighbor_list, i)
    end

    return neighbor_list
end

weight(g::WeightedDiAdjacencyMatrix) = g.relation

function probability(g::WeightedDiAdjacencyMatrix)
    g.relation ./ sum(g.relation, dims=2)
end

function Base.show(io::IO, g::AbstractGraph)
    println(io, "Graph(")
    for i = 1:g.n_vertices
        for j = 1:g.n_vertices print(io, "\t $(g.relation[i, j])") end
        println(io)
    end
    print(io, ")")
end

relate!(g::DiAdjacencyMatrix, v1::Int64, v2::Int64) = (g.relation[v1, v2] = true)

function relate!(g::WeightedDiAdjacencyMatrix{T}, v1::Int64, v2::Int64, w::T) where T
    (g.relation[v1, v2] = w)
end

function random_walk(g::WeightedDiAdjacencyMatrix{T}, x::Vector{T}, steps::Int64) where T
    p = probability(g)
    if steps > 0
        x = p'^steps * x
    end

    return x
end

mutable struct DiAdjacencyList <: AbstractGraph
    relation::Vector{Vector{Int64}}
    n_vertices::Int64

    DiAdjacencyList(n::Int64) = new([Int[] for i=1:n], n)
end

mutable struct WeightedAdjacencyList{T} <: AbstractGraph
    relation::Vector{Vector{Int64}}
    weight::Vector{Vector{T}}
    n_vertices::Int64

    function WeightedAdjacencyList{T}(n::Int64) where {T <: Number}
        return new(
            [Int[] for i=1:n],
            [T[] for i=1:n],
            n
        )
    end
end

const AdjacencyListLike = Union{DiAdjacencyList, WeightedAdjacencyList}

Base.getindex(g::AdjacencyListLike, i::Int64) = 0 < i < g.n_vertices ? i : -1

nv(g::AdjacencyListLike) = g.n_vertices

neighbor(g::AdjacencyListLike, v::Int64) = g.relation[v]

ne(g::DiAdjacencyList) = sum([length(g.relation[v]) for v in 1:(g.n_vertices)])

degree(g::DiAdjacencyList, v::Int64) = length(g.relation[v])

function Base.show(io::IO, g::DiAdjacencyList)
    print(io, "DiAdjacencyList(")
    for i=1:g.n_vertices
        print(io, "$i[")
        for j=g.relation[i]
            print(io, "$j, ")
        end
        print(io, "], ")
    end
    print(io, ")")
end

relate!(g::DiAdjacencyList, v1::Int64, v2::Int64) = push!(g.relation[v1], v2)

function unrelate!(g::DiAdjacencyList, v1::Int64, v2::Int64)
    deleteat!(g.relation[v1], findall(x -> x == v2, g.relation[v1]))
end

ne(g::WeightedAdjacencyList) = sum([length(g.relation[v]) for v in 1:(g.n_vertices)])/2

degree(g::WeightedAdjacencyList, v::Int64) = sum(g.weight[v])

function Base.show(io::IO, g::WeightedAdjacencyList{T}) where T
    print(io, "WeightedAdjacencyList{$T}(")
    for i=1:g.n_vertices
        print(io, "$i[")
        for (index, j)=enumerate(g.relation[i])
            print(io, "$j($(g.weight[i][index])), ")
        end
        print(io, "], ")
    end
    print(io, ")")
end

function relate!(g::WeightedAdjacencyList{T}, v1::Int64, v2::Int64, w::T) where T
    push!(g.relation[v1], v2)
    push!(g.relation[v2], v1)
    push!(g.weight[v1], w)
    push!(g.weight[v2], w)
end

function unrelate!(g::WeightedAdjacencyList{T}, v1::Int64, v2::Int64) where T
    indexes = findall(x->x==v2, g.relation[v1])
    deleteat!(g.relation[v1], indexes)
    deleteat!(g.weight[v1], indexes)
    indexes = findall(x->x==v1, g.relation[v2])
    deleteat!(g.relation[v2], indexes)
    deleteat!(g.weight[v2], indexes)
end

function Graph(n::Int64, representation::Symbol; random::Bool=true, T=Float64)
    if representation == :di_matrix
        return DiAdjacencyMatrix(n)
    elseif representation == :di_list
        return DiAdjacencyList(n)
    elseif representation == :w_di_matrix
        return WeightedDiAdjacencyMatrix{T}(n, random=random)
    elseif representation == :w_list
        return WeightedAdjacencyList{T}(n)
    end
end


text:
_winston_config = Inifile()

begin
    local fn
    for dir in [".";@__DIR__;LOAD_PATH]
        fn = joinpath(dir, "Winston.ini")
        if isfile(fn) break end
        fn = joinpath(dir, "Winston", "src", "Winston.ini")
        if isfile(fn) break end
    end
    read(_winston_config, fn)
end

function _atox(s::AbstractString)
    x = strip(s)
    if x == "nothing"
        return nothing
    elseif x == "true"
        return true
    elseif x == "false"
        return false
    elseif length(x) > 2 && lowercase(x[1:2]) == "0x"
        try
            h = parse(Int, x[3:end], base=16)
            return h
        catch
        end
    elseif x[1] == '{' && x[end] == '}'
        style = Dict{Symbol,Any}()
        pairs = map(strip, split(x[2:end-1], ',', keepempty=false))
        for pair in pairs
            kv = split(pair, ':', keepempty=false)
            style[ Symbol(strip(kv[1])) ] = _atox(strip(kv[2]))
        end
        return style
    elseif x[1] == '"' && x[end] == '"'
        return x[2:end-1]
    end
    if occursin(r"^[+-]?\d+$", x)
        return parse(Int, x)
    end
    r = tryparse(Float64, x)
    @static if VERSION < v"0.7.0-DEV.3017"
        return isnull(r) ? x : get(r)
    else
        return r === nothing ? x : r
    end
end

function config_value(section, option)
    strval = get(_winston_config, section, option, "nothing")
    _atox(strval)
end

function config_options(sec::AbstractString)
    opts = Dict{Symbol,Any}()
    if sec == "defaults"
        for (k,v) in _winston_config.defaults
            opts[Symbol(k)] = _atox(v)
        end
    elseif has_section(_winston_config, sec)
        for (k,v) in section(_winston_config, sec)
            opts[Symbol(k)] = _atox(v)
        end
    end
    opts
end


text:
using Pkg
Pkg.activate("environment")
Pkg.instantiate()

using SystemBenchmark, LinearAlgebra, LinearAlgebra.BLAS

OMPThreads = try
		parse(Int,ENV["OMP_NUM_THREADS"])
	catch
		1
end

JOBID = try
		ENV["JOB_ID"]
	catch
		"LOGIN"
end

clustername = try
		readchomp(`/shared/ucl/apps/cluster-bin/whereami`)
	catch
		"Unknown"
end

hname = try
                readchomp(`/usr/bin/hostname`)
	catch
		"Unknown"
end

BLAS.set_num_threads(OMPThreads)
actualthreads = BLAS.get_num_threads()

fname = string(clustername, "-", hname, "-", JOBID, "-", OMPThreads, "-", actualthreads, ".txt")

println("Running on cluster: ", clustername)
println("Running in job: ", JOBID)
println("Detected OpenMP thread count: ", OMPThreads)
println("Selected number of threads: ", actualthreads)
println("Writing output to: ", fname)

res = runbenchmark();

savebenchmark(fname, res)

max_stars_repo_name: owainkenwayucl/JuliaClusterBench

max_issues_repo_name: owainkenwayucl/JuliaClusterBench

max_forks_repo_name: owainkenwayucl/JuliaClusterBench


text:
module ViewHelper

using Genie, Genie.Flash, Genie.Router

export output_flash

function output_flash() :: String
  flash_has_message() ? """
$(flash())
""" : ""
end

end


text:
using FourierCore, FourierCore.Grid, FourierCore.Domain
using FFTW, LinearAlgebra, BenchmarkTools, Random, JLD2
using GLMakie
rng = MersenneTwister(1234)
# Random.seed!(123456789)
Random.seed!(12)
# jld_name = "high_order_timestep_spatial_tracer_"
jld_name = "blocky"
include("transform.jl")
include("random_phase_kernel.jl")
# using GLMakie
using CUDA
arraytype = CuArray
Ω = S¹(2π)^2
N = 2^9 # number of gridpoints
grid = FourierGrid(N, Ω, arraytype=arraytype)
nodes, wavenumbers = grid.nodes, grid.wavenumbers

x = nodes[1]
y = nodes[2]
kˣ = wavenumbers[1]
kʸ = wavenumbers[2]
# construct filter
kxmax = maximum(kˣ)
kymax = maximum(kˣ)
filter = @. (kˣ)^2 + (kʸ)^2 ≤ ((kxmax / 2)^2 + (kymax / 2)^2)
filter = @. abs(kˣ) .+ 0 * abs(kʸ) ≤ 2/3 * kxmax
@. filter = filter * ( 0 * abs(kˣ) .+ 1 * abs(kʸ) ≤ 2/3 * kxmax )

# now define the random field
wavemax = 5
𝓀 = arraytype([-wavemax, 0, wavemax]) # arraytype(1.0 .* [-wavemax, -wavemax + 1, wavemax - 1, wavemax])# arraytype(collect(-wavemax:1:wavemax))
𝓀ˣ = reshape(𝓀, (length(𝓀), 1))
𝓀ʸ = reshape(𝓀, (1, length(𝓀)))
# A = @. 0.1 * (𝓀ˣ * 𝓀ˣ + 𝓀ʸ * 𝓀ʸ)^(-11 / 12)
A = @. 0.5 * (𝓀ˣ * 𝓀ˣ + 𝓀ʸ * 𝓀ʸ)^(0.01) # @. 1e-1 / (1 * 2 * wavemax^2) .* (𝓀ˣ * 𝓀ˣ + 𝓀ʸ * 𝓀ʸ)^(0.0) # ( 1 .+ (0 .* 𝓀ˣ) .* 𝓀ʸ)
A[A.==Inf] .= 0.0
φ = arraytype(2π * rand(size(A)...))
field = arraytype(zeros(N, N))

##
# Fields
# velocity
ψ = arraytype(zeros(ComplexF64, N, N))
u = similar(ψ)
v = similar(ψ)

# theta
θ = similar(ψ)
∂ˣθ = similar(ψ)
∂ʸθ = similar(ψ)
κΔθ = similar(ψ)
θ̇ = similar(ψ)
s = similar(ψ)
θ̅ = similar(ψ)
k₁ = similar(ψ)
k₂ = similar(ψ)
k₃ = similar(ψ)
k₄ = similar(ψ)
θ̃ = similar(ψ)

# source
s = similar(ψ)
@. s = cos(kˣ[5] * x)

# phase
φ̇ = similar(A)

# operators
∂x = im * kˣ
∂y = im * kʸ
Δ = @. ∂x^2 + ∂y^2
Δ⁻¹ = 1 ./ Δ
bools = (!).(isnan.(Δ⁻¹))
Δ⁻¹ .*= bools # hack in the fact that false * NaN = 0

# plan ffts
P = plan_fft!(ψ)
P⁻¹ = plan_ifft!(ψ)

# number of gridpoints in transition is about λ * N / 2
bump(x; λ=10 / N, width=π / 2) = 0.5 * (tanh((x + width / 2) / λ) - tanh((x - width / 2) / λ))
bumps(x; λ=20 / N, width=1.0) = 0.25 * (bump(x, λ=λ, width=width) + bump(x, λ=λ, width=2.0 * width) + bump(x, λ=λ, width=3.0 * width) + bump(x, λ=λ, width=4.0 * width))

##
Δx = x[2] - x[1]

Δt = Δx / (2π) * 1

κ = 1.0 * Δx^2 # /Δt  # 1.0 / N * 0.1 # * 2^(2)  # roughly 1/N for this flow

# take the initial condition as negative of the source
tic = Base.time()

# save some snapshots
ψ_save = typeof(real.(Array(ψ)))[]
θ_save = typeof(real.(Array(ψ)))[]

r_A = Array(@. sqrt((x - π)^2 + (y - π)^2))
θ_A = [bumps(r_A[i, j]) - 3.0 for i in 1:N, j in 1:N]
θ .= CuArray(θ_A)
θ .= @. 0.1 * sin(3 * x) * sin(3 * y) + 0im
# @. θ = bump(sqrt(x^2 + y^2)) # scaling so that source is order 1
θclims = extrema(Array(real.(θ))[:])
P * θ # in place fft
@. κΔθ = κ * Δ * θ
P⁻¹ * κΔθ # in place fft
s .= -κΔθ * 0.0
P⁻¹ * θ # in place fft
θ̅ .= 0.0

t = [0.0]
tend = 40 # 5000

phase_speed = 1.0

iend = ceil(Int, tend / Δt)

#=
function θ_rhs_moisture!(θ̇, θ, params)
    #(; ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ, u, v, ∂ˣθ, ∂ʸθ, s, P, P⁻¹, filter) = params
    ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ, u, v, ∂ˣθ, ∂ʸθ, s, P, P⁻¹, filter = params
    τ = 0.01
    e = 0.01
    event = stream_function!(ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ)
    wait(event)
    P * ψ # in place fft
    P * θ # in place fft
    # ∇ᵖψ
    @. u = filter * -1.0 * (∂y * ψ)
    @. v = filter * (∂x * ψ)
    # ∇θ
    @. ∂ˣθ = filter * ∂x * θ
    @. ∂ʸθ = filter * ∂y * θ
    @. κΔθ = κ * Δ * θ # (κ * Δ  - (κ * Δ)^2 ) * θ
    # go back to real space
    P⁻¹ * ψ
    P⁻¹ * θ
    P⁻¹ * u
    P⁻¹ * v
    P⁻¹ * ∂ˣθ
    P⁻¹ * ∂ʸθ
    P⁻¹ * κΔθ
    # construct source
    @. s = v + e - 1 / τ * θ * (real(θ) > 0)
    # Assemble RHS
    @. θ̇ = -u * ∂ˣθ - v * ∂ʸθ + κΔθ + s
    return nothing
end
=#

function θ_rhs_moisture!(θ̇, θ, params)
    #(; ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ, u, v, ∂ˣθ, ∂ʸθ, s, P, P⁻¹, filter) = params
    ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ, u, v, ∂ˣθ, ∂ʸθ, s, P, P⁻¹, filter, Δ⁻¹ = params
    τ = 0.01
    e = 0.01
    event = stream_function!(ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ)
    wait(event)
    # @. ψ = cos(4 * x)
    @. s = ψ
    P * ψ # in place fft
    P * θ # in place fft
    # quick hack
    @. ψ = Δ⁻¹ * θ
    # ∇ᵖψ
    @. u = filter *  (∂y * ψ)
    @. v = filter * -1.0 * (∂x * ψ)
    # ∇θ
    @. ∂ˣθ = filter * ∂x * θ
    @. ∂ʸθ = filter * ∂y * θ
    @. κΔθ = (Δ⁻¹ + κ * Δ - 1e-1 * (κ * Δ)^2 + 1e-3 * (κ * Δ)^3 - 1e-5 * (κ * Δ)^4) * θ # κ * Δ - 0.1 * (κ * Δ)^2 + 0.01 * (κ * Δ)^3) * θ # κ * Δ * θ #
    # go back to real space
    P⁻¹ * ψ
    P⁻¹ * θ
    P⁻¹ * u
    P⁻¹ * v
    P⁻¹ * ∂ˣθ
    P⁻¹ * ∂ʸθ
    P⁻¹ * κΔθ
    # construct source
    # @. s = ψ# v + e - 1 / τ * θ * (real(θ) > 0)
    # Assemble RHS
    @. θ̇ = real(-u * ∂ˣθ - v * ∂ʸθ + κΔθ + s)
    @. θ = real(θ)
    #
    # P * θ̇    # in place fft
    # @. θ̇ = 0 * θ̇
    # P⁻¹ * θ̇
    return nothing
end

params = (; ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ, u, v, ∂ˣθ, ∂ʸθ, s, P, P⁻¹, filter, Δ⁻¹)

size_of_A = size(A)

realizations = 1
for j in 1:realizations
    # θ .= CuArray(θ_A)
    for i = 1:iend
        # fourth order runge kutta on deterministic part
        # keep ψ frozen is the correct way to do it here

        # the below assumes that φ is just a function of time
        θ_rhs_moisture!(k₁, θ, params)
        @. θ̃ = θ + Δt * k₁ * 0.5

        φ_rhs_normal!(φ̇, φ, rng)
        @. φ += phase_speed * sqrt(Δt / 2) * φ̇

        θ_rhs_moisture!(k₂, θ̃, params)
        @. θ̃ = θ + Δt * k₂ * 0.5
        θ_rhs_moisture!(k₃, θ̃, params)
        @. θ̃ = θ + Δt * k₃

        φ_rhs_normal!(φ̇, φ, rng)
        @. φ += phase_speed * sqrt(Δt / 2) * φ̇

        θ_rhs_moisture!(k₄, θ̃, params)
        @. θ += Δt / 6 * (k₁ + 2 * k₂ + 2 * k₃ + k₄)

        # update stochastic part
        # φ_rhs_normal!(φ̇, φ, rng)
        # @. φ += sqrt(Δt) * φ̇

        t[1] += Δt
        # save output

        #=
        if i % div(iend, 10) == 0
            println("Saving at i=", i)
            push!(ψ_save, Array(real.(ψ)))
            push!(θ_save, Array(real.(θ)))
            println("extrema are ", extrema(θ_save[end]))
            println("time is t = ", t[1])
        end

        if t[1] >= tstart
            θ̅ .+= Δt * θ
        end

        if i % div(iend, 100) == 0
            println("time is t = ", t[1])
            local toc = Base.time()
            println("the time for the simulation is ", toc - tic, " seconds")
            println("extrema are ", extrema(real.(θ)))
            println("on wavenumber index ", index_choice)
        end
        =#
        if i % div(iend, 100) == 0
            println("time is t = ", t[1])
            local toc = Base.time()
            println("the time for the simulation is ", toc - tic, " seconds")
            println("extrema are ", extrema(real.(θ)))
            # println("on wavenumber index ", index_choice)
        end

    end
    println("finished realization ", j)
    @. θ̅ += θ / realizations
end

toc = Base.time()
println("the time for the simulation was ", toc - tic, " seconds")

# calculate vorticity
P * ψ # in place fft
@. κΔθ = Δ * ψ
P⁻¹ * κΔθ # in place fft

x_A = Array(x)[:] .- 2π
θ_F = Array(real.(θ))
ψ_F = Array(real.(ψ))
ω_F = Array(real.(κΔθ))
θ̅_F = Array(real.(θ̅))

fig = Figure(resolution=(2048, 512))
ax1 = Axis(fig[1, 1], title="moisture t = 0")
ax2 = Axis(fig[1, 2], title="moisture t = " * string(tend))
ax3 = Axis(fig[1, 4], title="vorticity t = " * string(tend))
println("the extrema of the end field is ", extrema(θ_F))
println("the extrema of the ensemble average is ", extrema(θ̅_F))
colormap = :bone_1
# colormap = :nipy_spectral
colorrange = (-1, 0)
colorrange2 = extrema(θ_F)
colorrange3 = (-maximum(ω_F), maximum(ω_F))
heatmap!(ax1, x_A, x_A, θ_A, colormap=colormap, colorrange=colorrange2, interpolate=true)
hm = heatmap!(ax2, x_A, x_A, θ_F, colormap=colormap, colorrange=colorrange2, interpolate=true)
hm_e = heatmap!(ax3, x_A, x_A, ω_F, colormap=:balance, colorrange=colorrange3, interpolate=true)
Colorbar(fig[1, 3], hm, height=Relative(3 / 4), width=25, ticklabelsize=30, labelsize=30, ticksize=25, tickalign=1,)
Colorbar(fig[1, 5], hm_e, height=Relative(3 / 4), width=25, ticklabelsize=30, labelsize=30, ticksize=25, tickalign=1,)
display(fig)

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
abstract type AbstractBasis end

struct PolynomialBasis <: AbstractBasis end
struct TrigBasis <: AbstractBasis end
struct PolyTrigBasis <: AbstractBasis end

function n2_terms(p::Int)
    return Int(p*(p+1)/2)
end

function n3_terms(p::Int)
    return p*p
end

function linear(X::Array,p::Int,n::Int)
    return X
end

function quadratic(X::Array,p::Int,n::Int)
    np = n2_terms(p) # 6
    C = zeros(np,n)  # 6,3
    k = 1
    for i=1:p
        for j=i:p
            C[k,:] .= X[i,:] .* X[j,:]
            k += 1
        end
    end
    return C
end

function cubic(X::Array,p::Int,n::Int)
    k = 1
    np = n3_terms(p)
    C = zeros(np,n)
    for i = 1:p
        for j = 1:p
            C[k,:] .= X[i,:].^2 .* X[j,:]
            k += 1
        end
    end
    return C
end

function trig(X::Array,p::Int,n::Int)
    C = zeros(3*p,n)
    k = 1
    for i = 1:p
        C[k,:] .= sin.(X[i,:])
        k += 1
    end
    for i = 1:p
        C[k,:] .= cos.(X[i,:])
        k += 1
    end
    for i = 1:p
        C[k,:] .= tan.(X[i,:])
        k += 1
    end
    return C
end

function basis(X::Array,o::PolynomialBasis)
    s = size(X)
    ntsteps = s[2] #3
    p = s[1] #2
    # Since we consider unit, linear and quadratic, cubic terms terms.
    n = 1 + p + n2_terms(p) + n3_terms(p)
    θ = zeros(n,ntsteps)
    θ[1,:] .= ones(ntsteps)
    k = 1+p
    l = k+n2_terms(p)
    m = l+n3_terms(p)
    θ[2:k,:] .= linear(X,p,ntsteps)
    θ[k+1:l,:] .= quadratic(X,p,ntsteps)
    θ[l+1:m,:] .= cubic(X,p,ntsteps)
    return θ
end

function basis(X::Array,o::TrigBasis)
    s = size(X)
    p = s[1]
    n = s[2]
    θ = zeros(3*p,n)
    θ .= trig(X,p,n)
    return θ
end

function basis(X::Array,o::PolyTrigBasis)
    o1 = PolynomialBasis()
    o2 = TrigBasis()
    θ1 = basis(X,o1)
    θ2 = basis(X,o2)
    return [θ1;θ2]
end


text:
modifier(ability) = div(ability - 10, 2, RoundDown)

ability() = sum(sort!([roll() for _ in 1:4])[2:end])

mutable struct DNDCharacter
    strength::Int
    dexterity::Int
    intelligence::Int
    wisdom::Int
    charisma::Int
    constitution::Int
    hitpoints::Int
    function DNDCharacter()
        rolls = [ability() for _ in 1:6]
        new(rolls...,10 + modifier(rolls[end]))
    end
end

roll(n, d) = sum(rand(Base.OneTo(d),n))
roll(d) = roll(1,d)
roll() = roll(6)

function str2roll(str)
    try
        (n, d) = parse.(Int,split(str,'d';limit=2))
        roll(n,d)
    catch(e)
        isa(e,ArgumentError) && throw(ArgumentError("invalid dice roll format"))
        throw(e)
    end
end

""" `roll"4d20"` style rolling, just for fun"""
macro roll_str(str) :(str2roll($str)) end


text:
using MatrixNetworks

function subMins!(mat)
    mcols = minimum(mat, dims=1)
    for j = 1:size(mat,2), i = 1:size(mat,1)
        mat[i,j] -= mcols[j]
    end
    mrows = minimum(mat, dims=2)
    for j = 1:size(mat,2), i = 1:size(mat,1)
        mat[i,j] -= mrows[i]
    end
end

function run_hungarian!(mat)
    # subtract col minimum and then row minimum
    subMins!(mat)

    nrows = size(mat)[1]
    ncols = size(mat)[2]
    rows_marked = zeros(Bool,nrows)
    cols_marked = zeros(Bool,ncols)

    new_marked_row_ind = falses(nrows)
    new_marked_col_ind = falses(ncols)

    ind_marked_rows = zeros(Int64, nrows)
    ind_unmarked_rows = zeros(Int64, nrows)
    ind_marked_cols = zeros(Int64, ncols)
    ind_unmarked_cols = zeros(Int64, ncols)

    matzeros = findall(iszero, mat)
    nmatzeros = length(matzeros)
    max_nmatzeros = nmatzeros

    addition_vector = zeros(Int64, nrows)

    while true
        @views matzeros_sub = matzeros[1:nmatzeros]
        new_marked_col_ind[:] .= false
        new_marked_row_ind[:] .= false

        matching  = get_matching(matzeros_sub)

        if matching.cardinality == ncols
            return matching.match
        else
            sorted_matching = partialsort(matching.match, 1:matching.cardinality)

            # mark all rows which aren't matched
            rows_marked[:] .= false
            cols_marked[:] .= false
            last = 0
            for r=1:matching.cardinality
                if sorted_matching[r]-last != 1
                    for ri=last+1:sorted_matching[r]-1
                        rows_marked[ri] = true
                        new_marked_row_ind[ri] = true
                    end
                    last = sorted_matching[r]
                else
                    last = sorted_matching[r]
                end
            end
            for r=last+1:nrows
                rows_marked[r] = true
                new_marked_row_ind[r] = true
            end
        end

        changed = true
        while changed
            changed = false

            # mark cols
            @inbounds for rc in matzeros_sub
                @inbounds if !cols_marked[rc[2]]
                    @inbounds if new_marked_row_ind[rc[1]]
                        cols_marked[rc[2]] = true
                        new_marked_col_ind[rc[2]] = true
                    end
                end
            end

            # mark new rows
            @inbounds for c in 1:ncols
                r = matching.match[c]
                @inbounds if new_marked_col_ind[c] && !rows_marked[r]
                    rows_marked[r] = true
                    new_marked_row_ind[r] = true
                    changed = true
                end
            end
        end

        nmatzeros = 0
        for matzero in matzeros_sub
            if (rows_marked[matzero[1]] && cols_marked[matzero[2]]) || (!rows_marked[matzero[1]] && !cols_marked[matzero[2]])
                nmatzeros += 1
                if nmatzeros > max_nmatzeros
                    max_nmatzeros += 1
                    push!(matzeros, matzero)
                else
                    matzeros[nmatzeros] = matzero
                end
            end
        end

        ind_marked_rows_end = 0
        ind_unmarked_rows_end = 0
        for r=1:nrows
            if rows_marked[r]
                ind_marked_rows_end += 1
                ind_marked_rows[ind_marked_rows_end] = r
            else
                ind_unmarked_rows_end += 1
                ind_unmarked_rows[ind_unmarked_rows_end] = r
            end
        end

        ind_marked_cols_end = 0
        ind_unmarked_cols_end = 0
        for c=1:ncols
            if cols_marked[c]
                ind_marked_cols_end += 1
                ind_marked_cols[ind_marked_cols_end] = c
            else
                ind_unmarked_cols_end += 1
                ind_unmarked_cols[ind_unmarked_cols_end] = c
            end
        end

        @views ind_marked_rows_sub = ind_marked_rows[1:ind_marked_rows_end]
        @views ind_unmarked_rows_sub = ind_unmarked_rows[1:ind_unmarked_rows_end]
        @views ind_marked_cols_sub = ind_marked_cols[1:ind_marked_cols_end]
        @views ind_unmarked_cols_sub = ind_unmarked_cols[1:ind_unmarked_cols_end]

        #minimum where !cols_marked but rows_marked
        @views mat_ind_unmarked_cols_sub = mat[:,ind_unmarked_cols_sub]
        @views rmins = vec(minimum(mat_ind_unmarked_cols_sub[ind_marked_rows_sub,:], dims=2))
        @views min_val = minimum(rmins)

        # subtract minimum from this mask
        @inbounds for c in ind_unmarked_cols_sub
            for r in ind_marked_rows_sub
                mat[r,c] -= min_val
                if mat[r,c] == 0
                    nmatzeros += 1
                    if nmatzeros > max_nmatzeros
                        max_nmatzeros += 1
                        push!(matzeros, CartesianIndex(r,c))
                    else
                        matzeros[nmatzeros] = CartesianIndex(r,c)
                    end
                end
            end
        end

        # add minimum where !rows_marked but cols_marked
        addition_vector .= 0
        addition_vector[ind_unmarked_rows_sub] .= min_val
        @views mat_ind_marked_cols_sub = mat[:,ind_marked_cols_sub]
        for col in eachcol(mat_ind_marked_cols_sub)
            col .+= addition_vector
        end
    end
end

function get_matching(matzeros_sub)
    # create bi partite matching graph
    nmatzeros = length(matzeros_sub)
    ei = zeros(Int64, nmatzeros)
    ej = zeros(Int64, nmatzeros)

    i = 1
    for rc in matzeros_sub
        ej[i] = rc[1] # row
        ei[i] = rc[2] # col
        i += 1
    end

    matching = bipartite_matching(ones(Int64,nmatzeros), ei, ej)
    return matching
end

function main()
    mat = [1 2 3 3 1;
           2 3 5 2 2;
           3 1 4 1 5;
           4 5 2 4 3;
           5 4 1 5 4]
    matching = hungarian(mat)
    @show matching
end

# main()

max_stars_repo_name: opensourcesblog/hungarian

max_issues_repo_name: opensourcesblog/hungarian

max_forks_repo_name: opensourcesblog/hungarian


text:
include("common.jl")

macro include_testset(filename)
    @assert filename isa AbstractString
    quote
        @testset $(filename) begin
            include($(filename))
        end
    end
end

####
#### unit tests
####

@include_testset("test-trees.jl")
@include_testset("test-hamiltonian.jl")
@include_testset("test-NUTS.jl")
@include_testset("test-stepsize.jl")
@include_testset("test-mcmc.jl")
@include_testset("test-diagnostics.jl")
@include_testset("test-logging.jl")

####
#### sample correctness tests
####

@include_testset("sample-correctness-tests.jl")


text:
@generated function isulp(x::Vnum{N, ES}) where {N,ES}
  inc = increment(Vnum{N, ES})
  :( (@u(x) & $inc) != 0 )
end

@generated function Base.isempty(x::Valid{N,ES}) where {N,ES}
  inc = increment(Vnum{N,ES})
  quote
    ((x.lower == zero(Vnum{N,ES})) || (x.upper == zero(Vnum{N,ES}))) && (@u(x.lower) == @u(x.upper) + $inc)
  end
end

Base.isnan(x::Valid{N,ES}) where {N,ES} = isempty(x)

@generated function isallreals(x::Valid{N,ES}) where {N,ES}
  inc = increment(Vnum{N,ES})
  quote
    (x.lower != zero(Vnum{N,ES})) && (x.upper != zero(Vnum{N,ES})) && (@u(x.lower) == @u(x.upper) + $inc)
  end
end

istile(x::Valid{N,ES}) where {N,ES} = (x.lower == x.upper)

@doc """
    roundsinf(x::Valid)
    checks if the value "rounds infinity".  Table:

    ===================================
    [(+, 0])      |           true
    [(+, +])      |  lower > upper
    [(+, Inf)     |          false
    [(+, Inf]     |           true
    [(+, -])      |           true
    [Inf, 0])     |           true
    (Inf, 0])     |          false
    [Inf, +])     |           true
    (Inf, +])     |          false
    [Inf, -])     |           true
    (Inf, -])     |          false
    [Inf, Inf]    |           true
    (Inf, Inf]    |           true (ℝp)
    [Inf, Inf)    |           true (ℝp)
    (Inf, Inf)    |          false
    [(-, 0])      |          false
    [(-, +])      |          false
    [(-, Inf)     |          false
    [(-, Inf]     |           true
    [(-, -])      |  lower > upper
    [0, 0]        |          false
    [0, 0)        |          false (∅, may be mapped to noisy NaN in the future)
    (0, 0]        |          false (∅)
    (0, 0)        |           true
    [(0, +])      |          false
    [(0, Inf)     |          false
    [(0, Inf]     |           true
    [(0, -])      |           true
"""
function roundsinf(x::Valid{N,ES}) where {N,ES}
    (!isempty(x)) && ((@s(x.upper) < @s(x.lower)) || isinf(x.lower))
end

#amazing symmetrical structure to "roundsinf"
function containszero(x::Valid{N,ES}) where {N,ES}
    (!isempty(x)) && ((@u(x.upper) < @u(x.lower)) || x.lower == zero(Vnum{N,ES}))
end

function ispositive(x::Valid{N,ES}) where {N,ES}
  zero(@Int) < (@s x.lower) <= (@s x.upper)
end

function isnegative(x::Valid{N,ES}) where {N,ES}
  (@u Vnum{N,ES}(Inf))< (@u x.lower) <= (@u x.upper)
end

"""
  nonnegative(::Valid) is true if no values in x are negative.
"""
nonnegative(x::Valid{N,ES}) where {N,ES} = (x.lower <= zero(Vnum{N,ES})) && ((!isfinite(x.upper) || x.lower <= x.upper <= maxpos(Vnum{N,ES})))

"""
  nonpositive(::Valid) is true if no values in x are positive.
"""
nonpositive(x::Valid{N,ES}) where {N,ES} = (x.upper <= zero(Vnum{N,ES})) && ((!isfinite(x.lower) || minneg(Vnum{N,ES}) <= x.lower <= x.upper))

"""
  rounds_positive(::Valid) is true if x contains both zero, infinity, and passes through the positive numbers.
"""
rounds_positive(x::Valid{N,ES}) where {N,ES} = isfinite(x.lower) && (x.lower <= zero(Vnum{N,ES})) && (x.upper < x.lower)

"""
  rounds_negative(::Valid) is true if x contains both zero, infinity, and passes through the negative numbers.
"""
rounds_negative(x::Valid{N,ES}) where {N,ES} = isfinite(x.upper) && (x.upper >= zero(Vnum{N,ES})) && (x.lower > x.upper)

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
struct KernelRecursiveABC{Kx,K,HO} <: AbstractPointABC
    kernelABC::K
    kernelx::Kx
    herding_options::HO
    theta::Vector{Float64}
end
"""
    KernelRecursiveABC(; kernel, prior, K, delta = 1e-2,
                         herding_options = (maxtime = 10,
                                            method = :LD_LBFGS,
                                            restarts = 10),
                         kernelx = Kernel())

Creates a `KernelRecursiveABC` structure. `KernelRecursiveABC` iterates a
[`KernelABC`](@ref) step and a kernel herding step. `kernel`, `prior`, `K` and
`delta` determine the `KernelABC` step. `herding_options` options and `kernelx`
determine the herding step.

# Example
```julia
using Distributions
model(x) = [randn() .+ x for _ in 1:3]
data = [[3.], [2.9], [3.3]]
k = KernelRecursiveABC(kernel = Kernel(),
                       delta = 1e-8,
                       prior = MultivariateNormal([0.], [5.]),
                       K = 10^2)
result = run!(k, model, data, maxfevals = 10^3)
result.x
```

[Kajihara, T., Kanagawa, M., Yamazaki, K. and Fukumizu, K. (2018),
Kernel Recursive ABC: Point Estimation with Intractable Likelihood,
Proceedings of the 35th International Conference on Machine Learning, 80:2400-2409
](http://proceedings.mlr.press/v80/kajihara18a.html)
"""
function KernelRecursiveABC(; kernel, prior, K,
                            herding_options = NamedTuple(),
                            delta = 1e-2, kernelx = Kernel())
    KernelRecursiveABC(KernelABC(kernel = kernel, prior = prior, K = K, delta = delta),
                       kernelx,
                       merge((maxtime = 10, method = :LD_LBFGS, restarts = 10),
                             herding_options),
                       zeros(length(prior)))
end

function wrap_gradient(f)
    (x, g) -> begin
        if length(g) > 0
            res = DiffResults.DiffResult(0., g)
            ForwardDiff.gradient!(res, f, x)
            res.value
        else
            f(x)
        end
    end
end
function herding_objective(w, k, θ, newθ, idx)
    wrap_gradient(x -> begin
        res = sum(w[i] * k(x, θ[i]) for i in eachindex(θ))
        if idx > 0
            res -= 1/(idx + 1) * sum(k(x, newθ[i]) for i in 1:idx)
        end
        res
    end)
end
function herding(θ, w, k;
                 d = length(θ[1]),
                 lower = fill(-Inf, d),
                 upper = fill(Inf, d),
                 onlyfirst = false,
                 maxtime = 3,
                 xtol_rel = 1e-4,
                 method = :LD_LBFGS,
                 verbose = false,
                 restarts = 5)
    newθ = deepcopy(θ)
    opt = Opt(method, d)
    opt.lower_bounds = lower
    opt.upper_bounds = upper
    opt.maxtime = maxtime
    opt.xtol_rel = xtol_rel
    x0 = sum(w .* θ)
    for i in eachindex(θ)
        opt.max_objective = herding_objective(w, k, θ, newθ, i-1)
        maxmaxf = -Inf
        for j in 1:restarts
            maxf, maxx, ret = optimize(opt, clamp.(j == 1 ? x0 : θ[wsample(w)] .+ bandwidth(k)/10 * randn(d), lower, upper))
            verbose && ret != :SUCCESS && @show i maxf maxmaxf maxx θ[i] ret opt.numevals
            if maxf > maxmaxf
                maxmaxf = maxf
                newθ[i] .= maxx
            end
        end
        onlyfirst && return newθ[1]
    end
    newθ
end

function run!(rng::Random.AbstractRNG, k::KernelRecursiveABC, model, data;
              verbose = true, maxfevals,
              callback = () -> nothing)
    start_time = now()
    N = div(maxfevals, k.kernelABC.K)
    verbose && @printf "%s %15.s %10.s\n" "iteration" "elapsed" "fevals"
    particles = k.kernelABC.particles
    for i in 1:k.kernelABC.K
        rand!(rng, k.kernelABC.prior, particles[i])
    end
    update!(k.kernelx, pairwise_euclidean(particles))
    kernelabc!(k.kernelABC, model, data)
    k.theta .= mean(k.kernelABC)
    callback()
    verbose && @printf "%9.s %15.s %10.s\n" 1 round(now() - start_time, Second) k.kernelABC.K
    for i in 2:N
        newparticles = herding(particles, k.kernelABC.weights, k.kernelx;
                               lower = lower(k.kernelABC.prior),
                               upper = upper(k.kernelABC.prior),
                               k.herding_options...)
        for j in eachindex(particles)
            particles[j] .= newparticles[j]
        end
        k.theta .= newparticles[1]
        kernelabc!(k.kernelABC, model, data, updatekernel = false)
        callback()
        verbose && @printf "%9.s %15.s %10.s\n" i round(now() - start_time, Second) i * k.kernelABC.K
    end
    (x = herding(particles, k.kernelABC.weights, k.kernelx;
                 onlyfirst = true,
                 lower = lower(k.kernelABC.prior),
                 upper = upper(k.kernelABC.prior),
                 k.herding_options...),)
end


text:
module LossFunctions

using ..OceanTurbulenceParameterEstimation
using ..OceanTurbulenceParameterEstimation.Grids
using ..OceanTurbulenceParameterEstimation.Observations
using ..OceanTurbulenceParameterEstimation.Models

using Oceananigans
using Oceananigans: AbstractModel
using Oceananigans.Grids: RegularRectilinearGrid
using Oceananigans.Fields: CenterField, AbstractField, AbstractDataField, interior

using Statistics

export
    # forward_run.jl
    model_time_series,
    ForwardMap,

    # loss_functions.jl
    LossFunction,
    TimeSeriesAnalysis,
    TimeAverage,
    ValueProfileAnalysis,
    GradientProfileAnalysis,
    evaluate!,
    on_grid

include("utils.jl")
include("time_analysis.jl")
include("profile_analysis.jl")
include("field_weights.jl")
include("loss_function.jl")
include("forward_run.jl")

end # module


text:
# # GARTEUR SM-AG19 Testbed: Construction of the geometry

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

# ## Description

# This virtual test application is based on the test article
# used by the GARTEUR Structures & Materials Action Group 19
# which organized a Round Robin exercise where 12 European laboratories
# tested a single structure between 1995 and 1997. The benchmark structure
# was a laboratory structure built to simulate the dynamic behaviour
# of an aeroplane. The structure was initially built for a benchmark
# study on experimental modal analysis conducted by the
# Structures and Materials Action Group (SM-AG19) of the Group
# for Aeronautical Research and Technology in EURope (GARTEUR).
# The test-bed was designed and manufactured by ONERA, France.

# ![](IMAC97photo.png)

# ### References

# - [GARTEUR] Ground Vibration Test Techniques, compiled by A Gravelle, GARTEUR
#   Structures & Materials Action Group 19 Technical report TP-115, 1999.
# - [BW] Etienne Balmes, Jan R. Wright, GARTEUR GROUP ON GROUND VIBRATION
#   TESTING | RESULTS FROM THE TEST OF A SINGLE STRUCTURE BY 12 LABORATORIES IN
#   EUROPE, Proceedings of DETC'97, 1997 ASME Design Engineering Technical
#   Conferences, September 14-17, 1997, Sacramento, California.

# ## Goals

# - Show how to construct model from multiple connected beams.
# - Demonstrate the use of massless connectors.
# - Visualize the structure interactively.
#

##
# ## Geometry of the testbed airplane.

# The aluminum testbed was a rather simple structure which was reasonably
# dynamically representative of a simple airplane structure [GARTEUR](@ref
# References). It was composed of several beams simulating a fuselage with wings
# and a tail. Wing tip drums allowed to adjust bending and torsion frequencies
# similarly to airplane ones, with some very close modal frequencies.

# ![](garteur-geom.png)

# The script included below defines the geometry of the structure, the
# cross-sectional properties, the connectivity, and the location of the nodes.

include("garteur_geometry_tut.jl")

##
# ## Basic visualization

# Here we use the `PlotlyJS` plotting library, with some simple utilities for
# generating geometry of beams.
using PlotlyJS
using FinEtoolsFlexBeams.VisUtilModule: plot_solid, plot_space_box, render, default_layout_3d, save_to_json

# The colors are used to help distinguish between the individual parts of the model.
colors = [
"rgb(125, 155, 155)",  # 1 body
"rgb(125, 155, 155)",  # 2 wing
"rgb(125, 155, 155)",  # 3 drums
"rgb(125, 155, 155)",  # 4 vertical tail
"rgb(125, 155, 155)",  # 5 horizontal tail
"rgb(125, 155, 155)",  # 6 constraining plate
"rgb(15, 15, 155)",  # 7 massless connectors: wing beam to plate
"rgb(125, 15, 15)",  # 8 massless connectors: wing beam to drums
"rgb(125, 175, 15)",  # 9 massless connectors: fuselage to wing beam
"rgb(125, 175, 15)",  # 10 massless connectors: fuselage to tail
"rgb(15, 155, 155)",  # 11 massless connectors: sensors, point masses
]

# The geometry is defined in terms of "traces" (a trace is in the `PlotlyJS`
# parlance a graphical object; it could be a curve, surface, and a lot of other
# things). It is the name for a graphical object in `PlotlyJS`. The two points
# below define a box, which is helpful when setting the extents of the graphics
# display.
tbox = plot_space_box([[-1.2 * L -1.2 * L -1.2 * L]; [+1.2 * L +1.2 * L +1.2 * L]])
# For each finite element set in the array `fesa`, generate the graphics to
# represent that object.
traces = let traces = tbox
    for fes in fesa
        labl  = fes.label[1]
        tm = plot_solid(fens, fes; facecolor=colors[labl]);
        traces = cat(traces, tm; dims = 1)
    end
    traces
end
# The layout of the plot is defined with simple defaults.
layout = default_layout_3d(;width=900, height=900)
# Next, the graphics is rendered, and may be interacted with by zooming,
# panning, etc.
pl = render(traces; layout = layout)

##
# ## Visualizing the nodes

# In order to be able to discern the nodes we will reduce the opacity of the
# surfaces representing the beams, otherwise the nodes would be hidden by these
# surfaces. Otherwise the geometries defined in the same way as above.

tbox = plot_space_box([[-1.2 * L -1.2 * L -1.2 * L]; [+1.2 * L +1.2 * L +1.2 * L]])
traces = let traces = tbox
    for fes in fesa
        labl  = fes.label[1]
        tm = plot_solid(fens, fes; facecolor=colors[labl], opacity = 0.3);
        traces = cat(traces, tm; dims = 1)
    end
    traces
end
# Next we add to the "traces" the graphics representing all the nodes in the
# model as bright red dots.
using FinEtoolsFlexBeams.VisUtilModule: plot_nodes
traces = cat(traces, plot_nodes(fens; color = "rgb(255, 15, 5)"); dims = 1)

# Finally, the graphics is presented.
layout = default_layout_3d(;width=900, height=900)
pl = render(traces; layout = layout)

true


text:
function test_nlp_mi_cvx(s)
    global solver = s

    test_dir = dirname(@__FILE__)

    @testset "Convex MINLP Models" begin end
end


text:
function rrule(::typeof(partialsort), xs::AbstractVector, k::Union{Integer,OrdinalRange}; kwargs...)
    inds = partialsortperm(xs, k; kwargs...)
    ys = xs[inds]

    function partialsort_pullback(Δys)
        function partialsort_add!(Δxs)
            Δxs[inds] += Δys
            return Δxs
        end

        Δxs = InplaceableThunk(@thunk(partialsort_add!(zero(xs))), partialsort_add!)

        return NO_FIELDS, Δxs, DoesNotExist()
    end

    return ys, partialsort_pullback
end

function rrule(::typeof(sort), xs::AbstractVector; kwargs...)
    inds = sortperm(xs; kwargs...)
    ys = xs[inds]

    function sort_pullback(Δys)
        function sort_add!(Δxs)
            Δxs[inds] += Δys
            return Δxs
        end

        Δxs = InplaceableThunk(@thunk(sort_add!(zero(Δys))), sort_add!)

        return NO_FIELDS, Δxs
    end
    return ys, sort_pullback
end


text:
function read_dots_folds(input_file)
    coords = []
    folds = []
    for line in readlines(input_file)
        if ',' in line
            x,y = map(x->parse(Int,x), split(line, ','))
            push!(coords, [x,y])
        elseif '=' in line
            axis, value = split(split(line)[3], '=')
            value = parse(Int, value)
            push!(folds, [axis, value])
        end
    end
    return coords, folds
end

function fold(input_coords, fold_axis, fold_value)
    println("folding ",fold_axis, " ", fold_value)
    output_coords = []
    for coord in input_coords
        if fold_axis == "y" #fold up
            if coord[2] < fold_value
                push!(output_coords, coord)
            else
                folded_coord = [coord[1], 2*fold_value - coord[2]]
                push!(output_coords, folded_coord)
            end
        else # fold left
            if coord[1] < fold_value
                push!(output_coords, coord)
            else
                folded_coord = [2*fold_value - coord[1], coord[2]]
                push!(output_coords, folded_coord)
            end
        end
    end
    # deduplicate by putting into a set
    return Set(output_coords)
end

function test_folds()
    test_coords = [[3,4],[11,18]]
    @assert fold(test_coords, "x", 10) == Set([[3,4],[9,18]])
    @assert fold(test_coords, "y", 10) == Set([[3,4],[11,2]])
    new_dots = fold(test_coords, "y", 10)
    @assert fold(new_dots, "x", 6) == Set([[3,4],[1,2]])
end

test_folds()

function do_all_folds(coords, folds)
    for f in folds
        coords = fold(coords, f[1], f[2])
        println("number of coords ", length(coords))
    end
    return coords
end

function draw_grid(coords)
    xmax = maximum(x->x[1], coords)+1
    ymax = maximum(x->x[2], coords)+1
    grid = ones(String, ymax, xmax)

    # now draw on the coords
    for coord in coords
        grid[coord[2]+1,coord[1]+1] = "#"
    end
    return grid

end

#
input_dots, folds = read_dots_folds("input.txt")
#
output_dots = do_all_folds(input_dots, folds)
#
println("number of output dots ", length(output_dots))
#
draw_grid(output_dots)


text:
#! /usr/bin/env julia

# Copyright 2017 John Hanley.
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# The software is provided "AS IS", without warranty of any kind, express or
# implied, including but not limited to the warranties of merchantability,
# fitness for a particular purpose and noninfringement. In no event shall
# the authors or copyright holders be liable for any claim, damages or
# other liability, whether in an action of contract, tort or otherwise,
# arising from, out of or in connection with the software or the use or
# other dealings in the software.

# Answers the question: "what fraction of my files came from .deb packages?"
#
# usage:
#     src/backup_report.jl

using SQLite

function all_files(fspec)
    uid = ccall((:getuid, "libc"), Int32, ())
    pat = Regex("^([\\d-]+ \\d{2}:\\d{2})\\s+"
                * "([\\.\\d]+)\\s+"
                * "([rwx\\w-]+)\\s+"
                * "([\\w-]+)\\s+"
                * "(\\d+)\\s+"
                * "(\\d+)\\s+"
                * "(.*)")
    Task() do
        open(fspec) do fin
            for line in eachline(fin)
                m = match(pat, rstrip(line))
                if ! (m === nothing)
                    mtime, ctime, perm, owner, links, size, file = m.captures
                    if isfile(file)
                        st = stat(file)
                        # Mine, or world read, or (ordinary) group read.
                        isreadable = (st.uid == uid) ||
                            (stat(file).mode & 0o004) > 0 ||
                            ((stat(file).mode & 0o040) > 0 && !is_security_group(st.gid))
                        if isreadable
                            produce(file)
                        end
                    end
                end
            end
        end
    end
end

function is_security_group(gid)
    return gid == 42 ||  # shadow
        gid == 107       # mlocate
end

function report(db, all_files)
    sel = "select src, pkg" *
        "  from file_source" *
        "  where pathname = ?"

    for file in all_files

        result = SQLite.query(db, sel; values=[file])

        if size(result, 1) == 0
            println(file)
        end
    end
end

dir = "/tmp/file_source/"
db = SQLite.DB(dir * "pkg_contents.sqlite")
report(db, all_files(dir * "find.txt"))


text:
const true_params = MountainCarParameters(-0.0025, 3)
const range_of_values = 0:0.005:0.045
# const range_of_values = 0.03:0.001:0.045

const theta1_values = -0.0023:-0.00005:-0.0027
const theta2_values = 2.9:0.01:3.1
const gridsearch_data_path = "/Users/avemula1/Developer/TAML/Dalinar/data/gridsearch.jld"


text:
module ControlSimTools

using ComponentArrays
using DataStructures: CircularBuffer, DataStructures
using DiffEqBase: DEProblem, AbstractDiffEqFunction, AbstractTimeseriesSolution, isinplace
using DiffEqCallbacks: PeriodicCallback, SavedValues
using OffsetArrays: OffsetArray
using RecursiveArrayTools: RecursiveArrayTools, DiffEqArray, VectorOfArray
using UnPack

include("input_handling.jl")
include("ContinuousDiscreteSystems/hybrid_systems.jl")
include("ContinuousDiscreteSystems/discrete_buffer_variable.jl")

export apply_inputs
export ContinuousDiscreteProblem, DiscreteBufferVariable, DiscreteVariableCallback

end # module