text:
fundict = Dict()
catdict = Dict()
equivdict = Dict()

spacer(x) = replace(x, "\n" => " \n \n ")

function funadd(name::String, category)
  x -> (fundict[name] = x; catdict[name] = category; equivdict[name] = [name])
end

function funadd(name::Array{String}, category)
  x -> for n in name
      fundict[n] = x;
      catdict[n] = category;
      equivdict[n] = name
    end
end

## Repl
"Print a table of 10 feasible results (head and tail)." |> funadd(["show", "s"], "REPL")
"Print a table of all feasible results." |> funadd("showall", "REPL")
"Go back one command." |> funadd(["back", "b"], "REPL")
"Go forward one command. Only works if you have gone back first." |> funadd(["next", "n", "f"], "REPL")
"Empty the current variable space." |> funadd("clear", "REPL")
"Empty the current as well as the logicset space." |> funadd("clearall", "REPL")
"Save the current variable space for use with `restore`." |> funadd("preserve", "REPL")
"Restore the last saved variable space." |> funadd("restore", "REPL")
"Show command history with # feasible." |> funadd(["history", "h"], "REPL")
"List variables names." |> funadd(["keys","k"], "REPL")
"Toggles the dashboard printout - primarily used for debugging." |> funadd(["dashboard", "dash"], "REPL")

"""
Show all logical sets currently in working REPL memmory. Exiting the REPL:
`julia> returnactivelogicset()` returns the most recent set
`julia> returnlogicset()` returns a vector of all sets worked on int REPL
""" |> funadd(["logicset","ls"], "REPL")

"""
Returns documentation related to a function or command to the user.

**NOTE** `julia> help(...)` returns identical information.
###### Example
abstractlogic> ? >
  Operator: >

  a > b : a is greater than b
""" |> funadd(["?", "help"], "REPL")

"""
check feasibility of subsequent command
### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :2 3 3

`abstractlogic>` a > b ; b > c

a > b                    feasible outcomes 9 ✓           :3 2 3
b > c                    feasible outcomes 1 ✓✓          :3 2 1

`abstractlogic>` check a == 3
Check: a == 3 ... a == 3                         feasible outcomes 1 ✓✓          :3 2 1
true, 1 out of 1 possible combinations 'true'.
""" |> funadd("check", "REPL")

"""
search feasibility of {{i}} match
### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :2 3 3

`abstractlogic>` a > b ; b > c

a > b                    feasible outcomes 9 ✓           :3 2 3
b > c                    feasible outcomes 1 ✓✓          :3 2 1

`abstractlogic>` search == 3

Checking: a == 3
Checking: b == 3
Checking: c == 3

:a is a match with 1 feasible combinations out of 1.
:b is a not match with 0 feasible combinations out of 1.
:c is a not match with 0 feasible combinations out of 1.
""" |> funadd("search", "REPL")

"""
`a, b, c ∈ 1:3` | {{i}} != {{!i}}: add `a`, `b` and `c` to the feasible space while
constraining `a`, `b`, and `c` each to not equal each other.
""" |> funadd("||", "Generator")

## Generators
"""
`a, b ∈ 1:3`: add `a` and `b` to the feasible space with potential values 1, 2, 3
`a, b ∈ 1,2,3`: add `a` and `b` to the feasible space with potential values 1, 2, 3
`a, b ∈ up, down`: add `a` and `b` to the feasible space with potential values 'up', 'down'
""" |> funadd(["∈","in"], "Generator")

"""
`a, b, c ∈ 1:3` | {{i}} != {{!i}}: add `a`, `b` and `c` to the feasible space while
constraining `a`, `b`, and `c` each to not equal each other.
""" |> funadd("||", "Generator")

## Operators
"`a > b`  : a is greater than b" |> funadd(">", "Operator")
"`a >= b` : a is greater than or equal to b" |> funadd(">=", "Operator")
"`a == b` : a is equal to b" |> funadd(["==", "="], "Operator")
"`a < b`  : a is less than b" |> funadd("<", "Operator")
"`a <= b` : a is less than or equal to b" |> funadd("<=", "Operator")

"`a >> b`  : a is at least two greater than b. (`3 >> 1`) -> `true`" |>
  funadd(">>", "Operator")
"`a << b`  : a is at least two less than b. (`3 << 1`) -> `true`" |>
  funadd("<<", "Operator")

"`a & : a` is odd.\n`a & b, c`: a, b, and c is odd" |> funadd("&", "Operator")
"`a ! : a` is odd.\n`a ! b, c`: a, b, and c is even" |> funadd("!", "Operator")
"`a ^ b` checks if one side is odd and one side is even" |> funadd("^", "Operator")

"""
`a, b` on the left (or right) is treated as saying that anything on that side
must conform to whatever Operator follows.
""" |> funadd(",", "Operator")

"""
`a|b` on the left (or right) is treated as saying that at least one value on
that side must conform to whatever Operator follows.
""" |> funadd("|", "Operator")

"""
`{#,#}` or `{#}`: it is often helpful to qualify the number of matches that the `or` syntax can make.
Specifying `{#}` at the end of a `a = b|c` match will force `#` matches.
Specifying `{#1,#2}` will force a number of matches no less than `#1` and no more than `#2`.
#### Example
``abstractlogic>` a | b | c | d = 1 {2,3}`
Requires that at least two but no more than three of `a`, `b`, `c`, or `d` to equal 1.
""" |> funadd("{", "Operator")

## Superoperators

"""
`x &&& y`: both x and y must be `true`. Command end operator `;` is often preferred
to joint evaluators since it is more efficient. In `x ; y` efficiency is gained
by evaluating `x` first, reducing the feasible set, then evaluating `y`.
**Note** `x and y` is equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b &&& b > c
a > b &&& b > c                  feasible outcomes 1 ✓✓          :3 2 1
""" |> funadd(["&&&", "and"], "Superoperator")

"""
`x ||| y` either x or y must be `true`. If both are `false` then total expression will return false.
**Note** `x or y` is equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b ||| b > c
a > b ||| b > c                  feasible outcomes 17 ✓          :2 2 1
""" |> funadd(["|||", "or"], "Superoperator")

"""
`x !== y` either x or y must be `true` but not both. If both are `false` or
both are `true` then it will return `false`.
**Note** `x xor y` and `x ^^^ y` is equivalent (`x` and `y` can only be `true` or `false`)
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b ^^^ b > c
a > b ^^^ b > c                  feasible outcomes 16 ✓          :3 3 2
""" |> funadd(["^^^","xor","!=="], "Superoperator")

"""
`x === y` if `x = true` then `y = true` or if `x = false` then `y = false`.
**Note** `x iff y` and `x <=> y` are equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b === b > c
a > b === b > c                  feasible outcomes 11 ✓          :1 3 3
""" |> funadd(["===", "iff", "<=>"], "Superoperator")

"""
`x ==> y` if `x = true` then `y` must be `true`.
**Note** `if x then y` and `x then y` are equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b ==> b > c
a > b ==> b > c                  feasible outcomes 19 ✓          :3 3 3
""" |> funadd(["==>", "then", "if...then"], "Superoperator")

"""
`x <== y` if `y = true` then `x` must be `true`.
**Note** `x if y` is equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b <== b > c
a > b <== b > c                  feasible outcomes 19 ✓          :3 1 1
""" |> funadd(["<==", "if"], "Superoperator")

## Metaoperators

"""
`x &&&& y`: both x and y must be `true`. Command end operator `;` is often preferred
to joint evaluators since it is more efficient. In `x ; y` efficiency is gained
by evaluating `x` first, reducing the feasible set, then evaluating `y`.
**Note** `x AND y` is an equivalent expression.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b &&&& b > c
a > b &&&& b > c                  feasible outcomes 1 ✓✓          :3 2 1
""" |> funadd(["&&&&", "AND"], "Metaoperator")

"""
`x |||| y` either x or y must be `true`. If both are `false` then total expression will return false.
**Note** `x or y` is equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b |||| b > c
a > b ||| b > c                  feasible outcomes 17 ✓          :2 2 1
""" |> funadd(["||||", "OR"], "Metaoperator")

"""
`x !=== y` either x or y must be `true` but not both. If both are `false` or
both are `true` then it will return `false`.
**Note** `x xor y` and `x ^^^ y` is equivalent (`x` and `y` can only be `true` or `false`)
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b !=== b > c
a > b !=== b > c                  feasible outcomes 16 ✓          :3 3 2
""" |> funadd(["^^^^","XOR","!==="], "Metaoperator")

"""
`x ==== y` if `x = true` then `y = true` or if `x = false` then `y = false`.
**Note** `x IFF y` and `x <==> y` are equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b ==== b > c
a > b ==== b > c                  feasible outcomes 11 ✓          :1 3 3
""" |> funadd(["====", "IFF", "<==>"], "Metaoperator")

"""
`x ===> y` if `x = true` then `y` must be `true`.
**Note** `IF x THEN y` and `x THEN y` are equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b ===> b > c
a > b ===> b > c                  feasible outcomes 19 ✓          :3 3 3
""" |> funadd(["===>", "THEN", "IF...THEN"], "Metaoperator")

"""
`x <=== y` if `y = true` then `x` must be `true`.
**Note** `x IF y` is equivalent
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` a > b <=== b > c
a > b <=== b > c                  feasible outcomes 19 ✓          :3 1 1
""" |> funadd(["<===", "IF"], "Metaoperator")

## Wildcards
"""
* `{{i}}` the most basic wildcard reference is `{{i}}`. It can take any variable name from the set of variables so far defined.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} > 1
{{i}} > 1  >>> a > 1 >>> b > 1 >>> c > 1
                 Feasible Outcomes: 8    Perceived Outcomes: 8 ✓         :3 3 3
""" |> funadd(["{{i}}", "i"], "Wildcard")

"""
* `{{!i}}` any variable not `i`.
**NOTE** In general `!i` is not that helpful as there are few operations which `!i` does not return an empty set. And those cases in which it doesn't seem to be more efficiently expressed as `>i` or `<i`.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} == {{!i}}
{{i}} == {{!i}}
`>>>` a == b
`>>>` a == c
`>>>` b == a
`>>>` b == c
`>>>` c == a
`>>>` c == b
                 feasible outcomes 3 ✓           :3 3 3
""" |> funadd(["{{!i}}", "!i"], "Wildcard")

"""
`{{>i}}` any variable to the right of `i`.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} == {{>i}}
{{i}} == {{>i}}
`>>>` a == b
`>>>` a == c
`>>>` b == c
                 feasible outcomes 3 ✓           :3 3 3
""" |> funadd(["{{>i}}", ">i"], "Wildcard")

"""
`{{<i}}` any variable to the left of `i`.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} == {{<i}}
{{i}} == {{<i}}
`>>>` b == a
`>>>` c == a
`>>>` c == b
                 feasible outcomes 3 ✓           :3 3 3
""" |> funadd(["{{<i}}", "<i"], "Wildcard")

"""
`{{<<i}}` any variable at least two to the left of `i`.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} == {{<<i}}
{{i}} == {{<<i}}
`>>>` c == a
         feasible outcomes 9 ✓           :3 1 3
""" |> funadd(["{{<<i}}", "<<i"], "Wildcard")

"""
`{{>>i}}` any variable at least two to the right of `i`.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} == {{>>i}}
{{i}} == {{>>i}}
`>>>` a == c
         feasible outcomes 9 ✓           :2 1 2
""" |> funadd(["{{>>i}}", ">>i"], "Wildcard")

"""
`{{i+n}}`, `{{i-n}}` takes n steps right `+` or left `-` of `i`.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} != {{i+1}}
{{i}} != {{i+1}}
`>>>` a != b
`>>>` b != c
         feasible outcomes 12 ✓          :1 2 1
""" |> funadd(["{{i+n}}", "{{i-n}}", "i+n", "i-n"], "Wildcard")

"""
* `{{i+n!}}`, `{{i-n!}}` forces n steps right `+` or left `-` of `i`. If that
step cannot be taken the variable is replaced by `999` which with numeric
comparisons this will have different effects.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} = 1 ==> {{i+1!}} = 2
{{i}} = 1 ==> {{i+1!}} = 2
`>>>` a = 1 ==> b = 2
`>>>` b = 1 ==> c = 2
`>>>` c = 1 ==> 999 = 2
         feasible outcomes 12 ✓          :1 2 3
""" |> funadd(["{{i+n!}}", "{{i-n!}}", "i+n!", "i-n!"], "Wildcard")

"""
`{{N}}` evaluates every permutation given the wildcard set and evaluates the entire
expression as `true` only if the number of values across wildcards is equal to `N`.
**Note** `{{N}}` must appear at the end of the *superoperator* expression.
#### Example
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} = 2 {{2}}
{{i}} = 2 {{2}}
`>>>` a = 2
`>>>` b = 2
`>>>` c = 2
                 feasible outcomes 6 ✓           :2 2 3
Forces at least two variables from `a`, `b`, and `c` to equal `2`.
""" |> funadd(["{{N}}", "N"], "Wildcard")

"""
`{{n1,n2}}` evaluates every permutation given the wildcard set and evaluates the entire
expression as `true` only if the number of values across wildcards is at least `n1` and no
more than `n2`. If `n1` is empty then it is set to zero. If `n2` is empty it is set to a
large number.
**Note** `{{n1,n2}}` must appear at the end of the *superoperator* expression.
**Hint** `{{,n2}}` can be read as "at most" while `{{n1,}}` can be read "at least".
#### Example
```julia
`abstractlogic>` a,b,c ∈ 1:3
a,b,c ∈ 1:3              feasible outcomes 27 ✓          :1 1 3
`abstractlogic>` {{i}} = 3 {{1,2}}
{{i}} = 3 {{1,2}}
`>>>` a = 3
`>>>` b = 3
`>>>` c = 3
         feasible outcomes 18 ✓          :2 3 3
```
Forces at least one or two of the variables to equal `3`.
""" |> funadd(["{{n1,n2}}", "n1,n2", ",n2", "n1,"], "Wildcard")

"""
`{{j}}` `j` is exactly the same as `i` except that it ignores the parts of variable names that come after `.`s. So `a.1` just becomes `a` in the substitution mechanism.
#### Example
```julia
`abstractlogic>` a.1,a.2 ,b.1 ,b.2, c.1, c.2  ∈ 1:3
a.1,a.2 ,b.1 ,b.2, c.1, c.2  ∈ 1:3       feasible outcomes 729 ✓         :2 3 3 2 3 1
`abstractlogic>` {{j}}.1 < {{j}}.2
{{j}}.1 < {{j}}.2
`>>>` a.1 < a.2
`>>>` b.1 < b.2
`>>>` c.1 < c.2
         feasible outcomes 27 ✓          :1 2 1 2 1 3

`abstractlogic>` back
Last command: "a.1,a.2 ,b.1 ,b.2, c.1, c.2  ∈ 1:3" - Feasible Outcomes: 729     :1 3 2 1 2 3

`abstractlogic>` {{j}}.1 == {{j+1}}.2
{{j}}.1 == {{j+1}}.2
`>>>` a.1 == b.2
`>>>` b.1 == c.2
         feasible outcomes 81 ✓          :3 1 1 3 1 1
""" |> funadd(["{{j}}", "j"], "Wildcard")

## Process the Different files
repllist          = [v[1] for v in catdict if v[2]=="REPL"]     |> sort
generatorlist     = [v[1] for v in catdict if v[2]=="Generator"]     |> sort
operatorlist      = [v[1] for v in catdict if v[2]=="Operator"]      |> sort
superoperatorlist = [v[1] for v in catdict if v[2]=="Superoperator"] |> sort
metaoperatorlist  = [v[1] for v in catdict if v[2]=="Metaoperator"]  |> sort
wildcardlist      = [v[1] for v in catdict if v[2]=="Wildcard"]      |> sort

"""
    help(x)

A lookup related to `AbstractLogic` operators, commands, and wildcards.

Example
```
julia> help("")
`` not found. Search any of the following:

REPL Commands: ?, b, back, check, clear, clearall, dash, dashboard, f, h, help, history, k, keys, logicset, ls, n, next, preserve, restore, s, search, show, showall

Generators: in, ||, ∈

Operators: !, &, ,, <, <<, <=, =, ==, >, >=, >>, ^, {, |

Superoperators: !==, &&&, <==, <=>, ===, ==>, ^^^, and, if, if...then, iff, or, then, xor, |||

Metaoperators: !===, &&&&, <===, <==>, ====, ===>, AND, IF, IF...THEN, IFF, OR, THEN, XOR, ^^^^, ||||

Wildcards: !i, ,n2, <<i, <i, >>i, >i, N, i, i+n, i+n!, i-n, i-n!, j, n1,, n1,n2, {{!i}}, {{<<i}}, {{<i}}, {{>>i}}, {{>i}}, {{N}}, {{i+n!}}, {{i+n}}, {{i-n!}}, {{i-n}}, {{i}}, {{j}}, {{n1,n2}}

julia> help("i")

Wildcard: i

variantes: {{i}}, i

•    {{i}} the most basic wildcard reference is {{i}}. It can take any variable name from the set of variables so far defined.

Example
---------

abstractlogic> a,b,c ∈ 1:3

a,b,c ∈ 1:3 feasible outcomes 27 ✓ :1 1 3

abstractlogic> {{i}} > 1

{{i}} > 1 >>> a > 1 >>> b > 1 >>> c > 1

            Feasible Outcomes: 8    Perceived Outcomes: 8 ✓         :3 3 3

```
"""
function help(x::String)
  x = replace(x, r"^(\?|help)"=>"") |> strip
  if (haskey(fundict, x))
    also = (length(equivdict[x]) == 1 ? "" : "**variantes:** `" * join(equivdict[x],", ") * "` \n")
    "**$(catdict[x])**: `$x` \n " * also * " \n " * fundict[x] |>
    spacer |>  printmarkdown
    return
  end
  """
  *`$x` not found*. Search any of the following:
  REPL Commands: `$(join(repllist, "`, `"))`
  Generators: `$(join(generatorlist, "`, `"))`
  Operators: `$(join(operatorlist, "`, `"))`
  Superoperators: `$(join(superoperatorlist, "`, `"))`
  Metaoperators: `$(join(metaoperatorlist, "`, `"))`
  Wildcards: `$(join(wildcardlist, "`, `"))`
 """ |> spacer |> printmarkdown
end

# help("=")
# help(":")
# help("in")
# help("&&&")
# help("&&&")

help() = help("")


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb02id(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    KMAX  = 20
    LMAX  = 20
    MMAX = 20
    NMAX = 20
    RBMAX = 20
    RCMAX = 20
    LDB  = KMAX*MMAX
    LDC  = KMAX*MMAX
    LDTC = MMAX*KMAX
    LDTR = KMAX
    LDWORK = 2*NMAX*LMAX*( LMAX + KMAX ) + ( 6 + NMAX )*LMAX + MMAX*KMAX*( LMAX + 1 ) + RBMAX + RCMAX
    TC = Array{Float64,2}(undef, LDTC,LMAX)
    TR = Array{Float64,2}(undef, LDTR,NMAX*LMAX)
    B = Array{Float64,2}(undef, LDB,RBMAX)
    C = Array{Float64,2}(undef, LDC,RCMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    K = parse(BlasInt, vs[1])
    L = parse(BlasInt, vs[2])
    M = parse(BlasInt, vs[3])
    N = parse(BlasInt, vs[4])
    RB = parse(BlasInt, vs[5])
    RC = parse(BlasInt, vs[6])
    JOB =  vs[7][1]
if ( K<=0 || K>KMAX )
elseif ( L<=0 || L>LMAX )
elseif ( M<=0 || M>MMAX )
elseif ( N<=0 || N>NMAX )
elseif ( ( LSAME( JOB, 'O' ) || LSAME( JOB, 'A' ) ) && ( ( RB<=0 ) || ( RB>RBMAX ) ) )
elseif ( ( LSAME( JOB, 'U' ) || LSAME( JOB, 'A' ) ) && ( ( RC<=0 ) || ( RC>RCMAX ) ) )
else

    vs = String[]
    _isz,_jsz = (M*K,L)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       TC[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (K,N*L-L)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       TR[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
if ( LSAME( JOB, 'O' ) || LSAME( JOB, 'A' ) )

    vs = String[]
    _isz,_jsz = (M*K,RB)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
end # if
if ( LSAME( JOB, 'U' ) || LSAME( JOB, 'A' ) )

    vs = String[]
    _isz,_jsz = (N*L,RC)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       C[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
end # if
# interp call 1

    INFO = SLICOT.mb02id!(JOB, K, L, M, N, RB, RC, TC, TR, B, C)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
if ( LSAME( JOB, 'O' ) || LSAME( JOB, 'A' ) )
# interp output 1
    println(io,"B:")
    _nc = RB
    _nr = N*L
    show(io,"text/plain",B[1:_nr,1:_nc])
    println(io,)

end # if
if ( LSAME( JOB, 'U' ) || LSAME( JOB, 'A' ) )
# interp output 2
    println(io,"C:")
    _nc = RC
    _nr = M*K
    show(io,"text/plain",C[1:_nr,1:_nc])
    println(io,)

end # if
end # if
end # if
    close(f)
end # run_X()

ext: jl


text:
################################################################################
#####  result info

type PostgresResultInfo <: DatabaseResultInfo
    msg::UTF8String
    severity::UTF8String
    state::Symbol
    code::UTF8String
    primary::UTF8String
    detail::UTF8String
    hint::UTF8String
    pos::UTF8String
end

function PostgresResultInfo(ptr::Ptr{PGresult})

        s = p -> p == C_NULL ? "" : utf8(p)
        getstate = sym -> s(Libpq.PQresultErrorField(ptr, Libpq.error_field[sym]))

        msg = utf8(Libpq.PQresultErrorMessage(ptr))
        severity =  getstate(:severity)
        code =      getstate(:sqlstate)
        state =     Libpq.error_state[code[1:2]]
        primary =   getstate(:message_primary)
        detail =    getstate(:message_detail)
        hint =      getstate(:message_hint)
        pos =       getstate(:statement_position)
        PostgresResultInfo(msg, severity, state, code, primary, detail, hint, pos)
end

function Base.show(io::IO, info::PostgresResultInfo)
    println("""PostgresResultInfo(
        \tmsg:$(strip(info.msg))
        \tseverity:$(info.severity)
        \tstate:$(info.state)
        \tcode:$(info.code)
        \tprimary:$(info.primary)
        \tdetail:$(info.detail)
        \thint:$(info.hint)
        \tpos:$(info.pos)
)""")
end

function notice_callback(id::Ptr{Void},  ptr::Ptr{PGresult})
    i = PostgresResultInfo(ptr)
    if i.state in (:warning, :invalid_transaction_state)
        warn(i.msg)
    elseif i.state == :successful_completion
        info(i.msg)
    end
    C_NULL
end

const notice_callback_ptr = cfunction(
                                notice_callback,
                                Ptr{Void},
                                (Ptr{Void}, Ptr{PGresult},))

"""Errors that propigate from the Postgres server."""
type PostgresServerError <: DatabaseServerError
    info::PostgresResultInfo
end

function Base.showerror(io::IO, err::PostgresServerError)
    for line in split(strip(err.info.msg), "\n")
        println(io, line)
    end
end

################################################################################
#####  Result

"""Interface type to retrieve values from the server."""
type PostgresResult
    types::Vector{AbstractPostgresType}
    colnames::Vector{UTF8String}
    nrows::Integer
    ncols::Integer
    ptr::Nullable{Ptr}
end

function check_status(p::Ptr{PGresult})

    s = Libpq.PQresultStatus(p)
    code = get(Libpq.exec_status, s, nothing)

    if code == nothing
        error("unknown status $s, $msg")

    elseif code == :fatal_error
        info = PostgresResultInfo(p)
        Libpq.PQclear(p)
        throw(PostgresServerError(info))
    end
    code
end

function check_command(p::Ptr{PGresult})
    code = check_status(p)
    if code == :command_ok
        cmd = utf8(Libpq.PQcmdStatus(p))
        num = utf8(Libpq.PQcmdTuples(p))
        info("$cmd $num")
    end
end

function PostgresResult(p::Ptr{PGresult}, types::Dict)

    r = nothing
    code = check_status(p)
    if code in (:tuples_ok, :empty_query, :command_ok)
        #code == :command_ok ? nonfatal_error(p) : nothing
        oids = [Int(Libpq.PQftype(p, col)) for col in 0:(Libpq.PQnfields(p)-1)]
        types = [get(types, oid, types[0]) for oid in oids]
        colnames = [utf8(Libpq.PQfname(p, col)) for col in 0:(Libpq.PQnfields(p)-1)]
        nrows = Libpq.PQntuples(p)
        ncols = Libpq.PQnfields(p)
        r = PostgresResult(types, colnames, nrows, ncols, Nullable(p))
    else
        error("unhandled server code: $code")
    end
    check_command(p)
    r
end

function Base.show(io::IO, r::PostgresResult)
    t = join(["$(typeof(t.naval)), " for t in r.types])[1:end-2]
    print(io, "$(r.nrows)x$(r.ncols){$t} PostgresResult")
end

"""Frees resources held by server."""
function close(r::PostgresResult)
    if !isnull(r.ptr)
        Libpq.PQclear(get(r.ptr))
    end
    r.ptr = Nullable{ResultPtr}()
end

function unsafe_column{T}(
    ptr::Ptr{Libpq.PGresult},
    col::Int,
    nrows::Int,
    t::AbstractPostgresType{T}
    )

    mask = Vector{Bool}(nrows)
    vals = Vector{T}(nrows)
    for row in 1:nrows
        isnull = Libpq.PQgetisnull(ptr, row-1, col-1)
        @inbounds mask[row] = isnull
        if isnull==1
            @inbounds vals[row] = t.naval
        else
            v = Libpq.PQgetvalue(ptr, row-1, col-1)
            @inbounds vals[row] = unsafe_parse(t, v)
        end
    end
    if isa(t, PostgresEnumType)
        PooledDataArray(vals, collect(t.enumvals), mask)
    else
        DataArray(vals, mask)
    end
end

"""`column(result, col)`
Returns a DataArray of the given column number."""
function column(result::PostgresResult, col::Int)
    if !(1 <= col <= result.ncols)
        throw(BoundsError(result, (row, col)))
    end
    unsafe_column(get(result.ptr), col, result.nrows, result.types[col])
end

"""`row(result, row)`
Returns a Vector{Any} of the given row number."""
function row(result::PostgresResult, row::Int)
    #tuple([result[row,col] for col in (1:result.ncols)] ...)
    v = Vector{Any}(result.ncols)
    for col in 1:result.ncols
        v[col] = result[row, col]
    end
    v
end

function Base.getindex(result::PostgresResult, row::Int, col::Int)
    if !(1 <= row <= result.nrows && 1 <= col <= result.ncols)
        throw(BoundsError(result, (row, col)))
    end
    t = result.types[col]
    col -= 1
    row -= 1
    ptr = get(result.ptr)
    if Libpq.PQgetisnull(ptr, row, col)==1
        return Nullable{typeof(naval(t))}()
    else
        p = Libpq.PQgetvalue(ptr, row, col)
        return Nullable{typeof(naval(t))}(unsafe_parse(t, p))
    end
end

Base.getindex(result::PostgresResult, ::Colon, col::Int) =
    column(result, col)

Base.getindex(result::PostgresResult, row_::Int, ::Colon) =
    row(result, row_)

Base.length(r::PostgresResult) = r.nrows
Base.size(r::PostgresResult) = (r.nrows,r.ncols)
Base.isempty(r::PostgresResult) = (r.nrows==0 && r.ncols==0)

Base.start(r::PostgresResult) = 1
Base.endof(r::PostgresResult) = r.nrows
Base.done(r::PostgresResult, i::Int) = i >= r.nrows
Base.next(r::PostgresResult, i::Int) = (row(r, i), i+1)
#eltype

max_stars_repo_name: NCarson/Postgres

max_issues_repo_name: NCarson/Postgres

max_forks_repo_name: NCarson/Postgres


text:
module DirectSearch

include("./Types.jl")
include("./Cache.jl")
include("./Constraints.jl")
include("./Core.jl")
include("./StoppingCondition.jl")
include("./Mesh.jl")
include("./Search.jl")
include("./Poll.jl")
include("./Report.jl")

#Include preset methods here
include("./LTMADS.jl")
include("./OrthoMADS.jl")
#include("./GPS.jl")

include("./BiMADS.jl")

include("./test_utils.jl")

end # module


text:
## in mpmath module
## These are not right!!!
## see hyper and meijerg to indicate what needs to be done for these special function
## they really need to be coordinated with `Julia`'s as well.
mpmath_fns = (:hyp0f1,
           :hyp1f1, :hyp1f2,
           :hyp2f0, :hyp2f1, :hyp2f2, :hyp2f3,
           :hyp3f2,
           :hypercomb,
#           :meijerg,
           :bihyper,
           :hyper2d,
           :appellf1, :appellf2, :appellf3, :appellf4,
           :ber,:bei,:ker,:kei,
           :struveh,:struvel,
           :angerj,
           :webere,
           :coulombc,
           :legenp, :legenq,
           :chebyt, :chebyu,
           :pcfd, :pcfu, :pcfv, :pcfw,
           :lommels1, :lommels2,
           :coulombf, :coulombg,
           :hyperu,
           :whitm, :whitw,
           :scorergi, :scorerhi,
           :spherharm,
           :airyaizero, :airybizero,
           :besseljzero, :besselyzero
           )
for fn in mpmath_fns
    meth = string(fn)
    @eval ($fn)(xs::SymOrNumber...;kwargs...) = mpmath_meth($meth, xs...; kwargs...)
    eval(Expr(:export, fn))
end

## Initialize mpmath
## includes trying to find the module!
## automatic mappings may throw warning about strings, though it is expected these
## will be addressed by PyCall
function init_mpmath()
    try
        PyCall.mpmath_init()
    catch err
        return()
    end
    ## try to load mpmath module
    copy!(mpmath, PyCall.pyimport_conda("mpmath", "mpmath"))
    if mpmath != PyCall.PyNULL()
        ## ignore warnings for now...
        mpftype = mpmath["mpf"]
        pytype_mapping(mpftype, BigFloat) ## Raises warning!!!
        mpctype = mpmath["mpc"]
        pytype_mapping(mpctype, Complex{BigFloat})
    end

    ## Call a function in the mpmath module, giving warning and returning NaN if module is not found
    global mpmath_meth(meth, args...; kwargs...) = begin
        if isa(mpmath, Void)
            warn("The mpmath module of Python is not installed. http://docs.sympy.org/dev/modules/mpmath/setup.html#download-and-installation")
            return(Sym(NaN))
        end

        fn = mpmath[@compat(Symbol(meth))]
        ans = call_sympy_fun(fn, args...; kwargs...)
        ## make nicer...
        if isa(ans, Vector)
            ans = convert(Vector{Sym}, ans)
        end
        ans
    end
end


text:
# SPDX-License-Identifier: X11
# 2020-09-01
# Making Triangle (200pt)

function main()
  N = parse(Int, readline())
  L = parse.(Int, split(readline()))
  sort!(L)

  c = 0
  for i ∈ 1:N - 2
    for j ∈ i + 1:N - 1
      for k ∈ j + 1:N
        L[i] != L[j] != L[k] && L[i] + L[j] > L[k] && (c += 1)
      end
    end
  end

  println(c)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
include( "testingUtils.jl" )
include( "../src/bootstrapPosteriorEstimation.jl" )
include( "../src/directSolution.jl" )
using DataFrames

println( "running direct solution unit tests" );

predMatTraining, predMatEval, tTraining, tEval, errMatTraining = makeupPredictions();

predMat = convert( Matrix{Float64}, predMatTraining );
t       = Float64.( tTraining );

pUni    = zeros( Float64, 16 );
fill!( pUni, 1 / 16 );

pBayes       = bootstrapPosteriorEstimation( errMatTraining, 100, 10000 );
pDirectMSE   = directOptimNaiveMSE( predMat, t, pBayes );
pDirectHinge = directOptimHinge( predMat, t, pBayes );
yBayes       = predMatEval * pBayes;
yDirectMSE   = predMatEval * pDirectMSE;
yDirectHinge = predMatEval * pDirectHinge;
yUni         = predMatEval * pUni;

performanceBayes       = mean( lossFunctions.MSE( yBayes, tEval ) );
performanceDirectMSE   = mean( lossFunctions.MSE( yDirectMSE, tEval ) );
performanceDirectHinge = mean( lossFunctions.MSE( yDirectHinge, tEval ) );
performanceUni         = mean( lossFunctions.MSE( yUni, tEval ) );

#== test MSE optimization ==#
@test performanceBayes > performanceDirectMSE
@test performanceDirectHinge > performanceDirectMSE
@test performanceUni > performanceDirectMSE

performanceBayes       = mean( lossFunctions.hingeLoss( yBayes, tEval ) );
performanceDirectMSE   = mean( lossFunctions.hingeLoss( yDirectMSE, tEval ) );
performanceDirectHinge = mean( lossFunctions.hingeLoss( yDirectHinge, tEval ) );
performanceUni         = mean( lossFunctions.hingeLoss( yUni, tEval ) );

#== test hingeLoss optimization ==#
@test performanceBayes > performanceDirectHinge
@test performanceDirectMSE > performanceDirectHinge
@test performanceUni > performanceDirectHinge


text:
############################### ODE Integrators ################################
# restrict to ODE integrators? could use for gradient flow.
function simulate(U::Update, x, ti::Real, tf::Real)
    t = ti
    while t < tf
        x, t = U(x, t)
    end
    x, t
end

################################# Runge Kutta ##################################
# TODO: adaptive step size, have to make them mutable?
struct Euler{T, F} <: Update{T}
    f::F
    Δ::T
end
const RK1 = Euler
(E::Euler)(x, t) = update(E, x, t)
function update(R::Euler, x, t::Real)
    x += R.Δ * R.f(t, x) # TODO: make t second argument
    t += R.Δ
    return x, t
end

struct Midpoint{T, F} <: Update{T}
    f::F
    Δ::T
end
const RK2 = Midpoint
(R::RK2)(x, t) = update(R, x, t)
function update(R::Midpoint, x, t::Real)
    x += R.Δ * R.f(t + R.Δ/2, x + R.Δ/2 * R.f(t, x))
    t += R.Δ
    return x, t
end

struct RungeKutta4{T, F} <: Update{T}
    f::F
    Δ::T
end
const RK4 = RungeKutta4
(R::RK4)(x, t) = update(R, x, t)
function update(R::RK4, x, t::Real)
    k1 = R.Δ * R.f(t, x) # explitic euler
    k2 = R.Δ * R.f(t + R.Δ/2, x + k1/2) # midpoint
    k3 = R.Δ * R.f(t + R.Δ/2, x + k2/2) # midpoint
    k4 = R.Δ * R.f(t + R.Δ, x + k3) # midpoint
    x += (k1 + 2k2 + 2k3 + k4)/6 # round-off lower if we factor out R.Δ?
    t += R.Δ
    return x, t
end

# generalized Runge-Kutta method where A, b, c constitute Butcher tableau
# make N part of type?
struct RungeKutta{T, F, AT<:AbstractMatrix{T}, B, C} <: Update{T}
    f::F
    A::AT # coefficient matrix
    b::B # vector of weightings for intermediate steps, for adaptive, can be matrix (with lower order method embedded)
    c::C # vector of fractional time-steps
    Δ::T
end
const RK = RungeKutta
# general implicit scheme
function update(R::RK, x::Real, t::Real)
    g(k) = sum(abs2, k - @. R.f(t + R.Δ * c, x + R.Δ * A*k)) # equation for k
    k = zero(R.b)
    fixedpoint!(ArmijoStep(Newton(g, k)), k) # if f is linear in k, this is one step
    x += R.Δ * dot(R.b, k)
    t += R.Δ
    return x, t
end
# general explicit scheme (i.e. Butcher tableau is lower triangular)
function update(R::RK{T, F, <:LowerTriangular}, x::Real, t::Real) where {T, F}
    k = zero(R.b)
    for i in eachindex(k)
        k[i] = R.f(t + R.Δ * c[i], x + R.Δ * dot(view(A, i, :), k))
    end
    x += R.Δ * dot(R.b, k)
    t += R.Δ
    return x, t
end

using StaticArrays
# TODO: add second b vector to embed lower order method for error estimation
# TODO: step sizes
function RK(f, Δ::Real, n::Int, explicit = true)
    if explicit
        A = zeros(n, n)
        if n == 1
            b = SVector(1)
            c = SVector(0)
        elseif n == 2
            A[2, 2] = 1/2
            b = SVector(0, 1)
            c = SVector(0, 1/2)
        elseif n == 4
            A[2, 1] = 1/2
            A[3, 2] = 1/2
            A[4, 3] = 1
            A = SMatrix{n,n}(0)
            b = SVector(1/6, 1/3, 1/3, 1/6)
            c = SVector(0, 1/2, 1/2, 1)
        else
            error("order $n RK method not implemented")
        end
        A = LowerTriangular(SMatrix{n,n}(A))
        RK(f, A, b, c, Δ)
    else
        error("implicit RK not implemented")
    end
end


text:
function remove_lln!(ex::Expr)
    for i in length(ex.args):-1:1
        if ex.args[i] isa LineNumberNode
            deleteat!(ex.args, i)
        elseif ex.args[i] isa Expr
            remove_lln!(ex.args[i])
        end
    end
end

macro enter(command)
    remove_lln!(command)
    :(JSONRPC.send_notification(conn_endpoint[], "debugger/enter", $(string(command))))
end

macro run(command)
    remove_lln!(command)
    :(JSONRPC.send_notification(conn_endpoint[], "debugger/run", $(string(command))))
end

max_stars_repo_name: fonsp/julia-vscode

max_issues_repo_name: fonsp/julia-vscode

max_forks_repo_name: fonsp/julia-vscode


text:
# Code for Figure 2.2
module Sutton2017

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

    include( "./fig_2_2.jl" )
    include( "./fig_2_4.jl" )

export
    run
end


text:
import Test: @test, @testset, @test_broken
import Luna: Maths, Nonlinear, PhysData
import FFTW

# test full / nothg / envelope kerr effect
Nt = collect(range(0, length=2^16))
t = @. (Nt - 2^16/2)*3.430944979182369e-16/4
E = exp.(-0.5.*(t./10e-15).^2).*cos.(2π*PhysData.c/800e-9.*t)
kerrfield = Nonlinear.Kerr_field(1.0)
outf = similar(E)
fill!(outf, 0.0)
kerrfield(outf, E, 1.0)
kerrfieldn = Nonlinear.Kerr_field_nothg(1.0, 2^16)
outn = similar(E)
fill!(outn, 0.0)
kerrfieldn(outn, E, 1.0)
kerrfieldenv = Nonlinear.Kerr_env(1.0)
Eenv = Maths.hilbert(E)
oute = similar(Eenv)
fill!(oute, 0.0)
kerrfieldenv(oute, Eenv, 1.0)

outfω = FFTW.rfft(outf)
outnω = FFTW.rfft(outn)
outeω = FFTW.rfft(real.(oute))
# we compare only low (non THG) frequencies
# note that these are not expected to be exact, because we have dropped not just the THG term
# but also cross terms between positive and negative frequencies
@test isapprox(abs.(outnω[1800:2400]), abs.(outfω[1800:2400]), rtol=1e-15)
@test isapprox(abs.(outeω[1800:2400]), abs.(outfω[1800:2400]), rtol=1e-15)


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module DataFrames

##############################################################################
##
## Dependencies
##
##############################################################################

using Compat
using Reexport
@reexport using StatsBase
@reexport using DataArrays
using GZip
using SortingAlgorithms
using Base: Sort, Order
using Docile

@document

##############################################################################
##
## Exported methods and types (in addition to everything reexported above)
##
##############################################################################

export @~,

       AbstractDataFrame,
       DataFrame,
       DataFrameRow,
       Formula,
       GroupApplied,
       GroupedDataFrame,
       ModelFrame,
       ModelMatrix,
       SubDataFrame,

       aggregate,
       by,
       coefnames,
       colwise,
       combine,
       complete_cases,
       complete_cases!,
       deleterows!,
       describe,
       eachcol,
       eachrow,
       eltypes,
       groupby,
       melt,
       meltdf,
       names!,
       ncol,
       nonunique,
       nrow,
       nullable!,
       order,
       pool,
       pool!,
       printtable,
       read_rda,
       readtable,
       rename!,
       rename,
       showcols,
       stack,
       stackdf,
       unique!,
       unstack,
       writetable

##############################################################################
##
## Load files
##
##############################################################################

for (dir, filename) in [
        ("other", "utils.jl"),
        ("other", "index.jl"),

        ("abstractdataframe", "abstractdataframe.jl"),
        ("dataframe", "dataframe.jl"),
        ("subdataframe", "subdataframe.jl"),
        ("groupeddataframe", "grouping.jl"),
        ("dataframerow", "dataframerow.jl"),

        ("abstractdataframe", "iteration.jl"),
        ("abstractdataframe", "join.jl"),
        ("abstractdataframe", "reshape.jl"),

        ("abstractdataframe", "io.jl"),
        ("dataframe", "io.jl"),

        ("abstractdataframe", "show.jl"),
        ("groupeddataframe", "show.jl"),
        ("dataframerow", "show.jl"),

        ("abstractdataframe", "sort.jl"),
        ("dataframe", "sort.jl"),

        ("statsmodels", "formula.jl"),
        ("statsmodels", "statsmodel.jl"),

        ("", "RDA.jl"),
        ("", "deprecated.jl")
    ]

    include(joinpath(dir, filename))
end

end # module DataFrames


text:
struct Product{A,B}
    a::A
    b::B
end
@inline function Base.size(p::Product)
    M = size(p.a, 1)
    (M, Base.tail(size(p.b))...)
end
@inline function Base.size(p::Product, i::Integer)
    i == 1 && return size(p.a, 1)
    size(p.b, i)
end
@inline Base.length(p::Product) = prod(size(p))
@inline Base.broadcastable(p::Product) = p
@inline Base.ndims(p::Type{Product{A,B}}) where {A,B} = ndims(B)

Base.Broadcast._broadcast_getindex_eltype(::Product{A,B}) where {T, A <: AbstractVecOrMat{T}, B <: AbstractVecOrMat{T}} = T
function Base.Broadcast._broadcast_getindex_eltype(p::Product)
    promote_type(
        Base.Broadcast._broadcast_getindex_eltype(p.a),
        Base.Broadcast._broadcast_getindex_eltype(p.b)
    )
end

# recursive_eltype(::Type{A}) where {T, A <: AbstractArray{T}} = T
# recursive_eltype(::Type{NTuple{N,T}}) where {N,T<:Union{Float32,Float64}} = T
# recursive_eltype(::Type{Float32}) = Float32
# recursive_eltype(::Type{Float64}) = Float64
# recursive_eltype(::Type{Tuple{T}}) where {T} = T
# recursive_eltype(::Type{Tuple{T1,T2}}) where {T1,T2} = promote_type(recursive_eltype(T1), recursive_eltype(T2))
# recursive_eltype(::Type{Tuple{T1,T2,T3}}) where {T1,T2,T3} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3))
# recursive_eltype(::Type{Tuple{T1,T2,T3,T4}}) where {T1,T2,T3,T4} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3), recursive_eltype(T4))
# recursive_eltype(::Type{Tuple{T1,T2,T3,T4,T5}}) where {T1,T2,T3,T4,T5} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3), recursive_eltype(T4), recursive_eltype(T5))

# function recursive_eltype(::Type{Broadcasted{S,A,F,ARGS}}) where {S,A,F,ARGS}
#     recursive_eltype(ARGS)
# end

@inline *ˡ(a::A, b::B) where {A,B} = Product{A,B}(a, b)
@inline Base.Broadcast.broadcasted(::typeof(*ˡ), a::A, b::B) where {A, B} = Product{A,B}(a, b)
const ∗ = *ˡ
# TODO: Need to make this handle A or B being (1 or 2)-D broadcast objects.
function add_broadcast!(
    ls::LoopSet, mC::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Product{A,B}}, elementbytes::Int = 8
) where {A, B}
    K = gensym(:K)
    mA = gensym(:Aₘₖ)
    mB = gensym(:Bₖₙ)
    pushpreamble!(ls, Expr(:(=), mA, Expr(:(.), bcname, QuoteNode(:a))))
    pushpreamble!(ls, Expr(:(=), mB, Expr(:(.), bcname, QuoteNode(:b))))
    pushpreamble!(ls, Expr(:(=), K, Expr(:call, :size, mB, 1)))

    k = gensym(:k)
    ls.loops[k] = Loop(k, K)
    m = loopsyms[1];
    if ndims(B) == 1
        bloopsyms = Symbol[k]
        cloopsyms = Symbol[m]
        reductdeps = Symbol[m, k]
    elseif ndims(B) == 2
        n = loopsyms[2];
        bloopsyms = Symbol[k,n]
        cloopsyms = Symbol[m,n]
        reductdeps = Symbol[m, k, n]
    else
        throw("B must be a vector or matrix.")
    end
    # load A
    # loadA = add_load!(ls, gensym(:A), productref(A, mA, m, k), elementbytes)
    loadA = add_broadcast!(ls, gensym(:A), mA, Symbol[m,k], A, elementbytes)
    # load B
    loadB = add_broadcast!(ls, gensym(:B), mB, bloopsyms, B, elementbytes)
    # set Cₘₙ = 0
    # setC = add_constant!(ls, zero(promote_type(recursive_eltype(A), recursive_eltype(B))), cloopsyms, mC, elementbytes)
    setC = if elementbytes == 4
        add_constant!(ls, 0f0, cloopsyms, mC, elementbytes)
    else#if elementbytes == 4
        add_constant!(ls, 0.0, cloopsyms, mC, elementbytes)
    end
    # compute Cₘₙ += Aₘₖ * Bₖₙ
    reductop = Operation(
        ls, mC, elementbytes, :vmuladd, compute, reductdeps, Symbol[k], Operation[loadA, loadB, setC]
    )
    pushop!(ls, reductop, mC)
end

struct LowDimArray{D,T,N,A<:DenseArray{T,N}} <: DenseArray{T,N}
    data::A
end
@inline Base.pointer(A::LowDimArray) = pointer(A.data)
Base.size(A::LowDimArray) = Base.size(A.data)
@generated function VectorizationBase.stridedpointer(A::LowDimArray{D,T,N}) where {D,T,N}
    s = Expr(:tuple, [Expr(:ref, :strideA, n) for n ∈ 1+D[1]:N if D[n]]...)
    f = D[1] ? :PackedStridedPointer : :SparseStridedPointer
    Expr(:block, Expr(:meta,:inline), Expr(:(=), :strideA, Expr(:call, :strides, Expr(:(.), :A, QuoteNode(:data)))),
         Expr(:call, Expr(:(.), :VectorizationBase, QuoteNode(f)), Expr(:call, :pointer, Expr(:(.), :A, QuoteNode(:data))), s))
end
function LowDimArray{D}(data::A) where {D,T,N,A <: AbstractArray{T,N}}
    LowDimArray{D,T,N,A}(data)
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{<:LowDimArray{D,T,N}}, elementbytes::Int = 8
) where {D,T,N}
    fulldims = Union{Symbol,Int}[loopsyms[n] for n ∈ 1:N if D[n]]
    ref = ArrayReference(bcname, fulldims, Ref{Bool}(false))
    add_load!(ls, destname, ref, elementbytes)::Operation
end
function add_broadcast_adjoint_array!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{A}, elementbytes::Int = 8
) where {T,N,A<:AbstractArray{T,N}}
    parent = gensym(:parent)
    pushpreamble!(ls, Expr(:(=), parent, Expr(:call, :parent, bcname)))
    ref = ArrayReference(parent, Union{Symbol,Int}[loopsyms[N + 1 - n] for n ∈ 1:N], Ref{Bool}(false))
    add_load!( ls, destname, ref, elementbytes )::Operation
end
function add_broadcast_adjoint_array!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{<:AbstractVector}, elementbytes::Int = 8
)
    ref = ArrayReference(bcname, Union{Symbol,Int}[loopsyms[2]], Ref{Bool}(false))
    add_load!( ls, destname, ref, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Adjoint{T,A}}, elementbytes::Int = 8
) where {T, A <: AbstractArray{T}}
    add_broadcast_adjoint_array!( ls, destname, bcname, loopsyms, A, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Transpose{T,A}}, elementbytes::Int = 8
) where {T, A <: AbstractArray{T}}
    add_broadcast_adjoint_array!( ls, destname, bcname, loopsyms, A, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{<:AbstractArray{T,N}}, elementbytes::Int = 8
) where {T,N}
    add_load!(ls, destname, ArrayReference(bcname, @view(loopsyms[1:N]), Ref{Bool}(false)), elementbytes)
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{T}, elementbytes::Int = 8
) where {T<:Union{Integer,Float32,Float64}}
    op = add_constant!(ls, destname, elementbytes) # or replace elementbytes with sizeof(T) ? u
    pushpreamble!(ls, Expr(:(=), mangledvar(op), bcname))
    op
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{SubArray{T,N,A,S,B}}, elementbytes::Int = 8
) where {T,N,N2,A<:AbstractArray{T,N2},B,N3,S <: Tuple{Int,Vararg{Any,N3}}}
    inds = Vector{Union{Int,Symbol}}(undef, N+1)
    inds[1] = Symbol("##DISCONTIGUOUSSUBARRAY##")
    inds[2:end] .= @view(loopsyms[1:N])
    add_load!(ls, destname, ArrayReference(bcname, inds, Ref{Bool}(false)), elementbytes)
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Broadcasted{S,Nothing,F,A}},
    elementbytes::Int = 8
) where {N,S<:Base.Broadcast.AbstractArrayStyle{N},F,A}
    instr = get(FUNCTIONSYMBOLS, F) do
        Instruction(bcname, :f)
    end
    args = A.parameters
    Nargs = length(args)
    bcargs = Expr(:(.), bcname, QuoteNode(:args))
    # this is the var name in the loop
    parents = Operation[]
    deps = Symbol[]
    reduceddeps = Symbol[]
    for (i,arg) ∈ enumerate(args)
        argname = gensym(:arg)
        pushpreamble!(ls, Expr(:(=), argname, Expr(:ref, bcargs, i)))
        # dynamic dispatch
        parent = add_broadcast!(ls, gensym(:temp), argname, loopsyms, arg, elementbytes)::Operation
        pushparent!(parents, deps, reduceddeps, parent)
    end
    op = Operation(
        length(operations(ls)), destname, elementbytes, instr, compute, deps, reduceddeps, parents
    )
    pushop!(ls, op, destname)
end

# size of dest determines loops
# function vmaterialize!(
@generated function vmaterialize!(
    dest::AbstractArray{T,N}, bc::BC
) where {T <: Union{Float32,Float64}, N, BC <: Broadcasted}
    # we have an N dimensional loop.
    # need to construct the LoopSet
    loopsyms = [gensym(:n) for n ∈ 1:N]
    ls = LoopSet()
    sizes = Expr(:tuple)
    for (n,itersym) ∈ enumerate(loopsyms)
        Nsym = gensym(:N)
        ls.loops[itersym] = Loop(itersym, Nsym)
        push!(sizes.args, Nsym)
    end
    pushpreamble!(ls, Expr(:(=), sizes, Expr(:call, :size, :dest)))
    elementbytes = sizeof(T)
    add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
    add_store!(ls, :dest, ArrayReference(:dest, loopsyms, Ref{Bool}(false)), elementbytes)
    resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
    q = lower(ls)
    push!(q.args, :dest)
    pushfirst!(q.args, Expr(:meta,:inline))
    q
    # ls
end
@generated function vmaterialize!(
    dest′::Union{Adjoint{T,A},Transpose{T,A}}, bc::BC
) where {T <: Union{Float32,Float64}, N, A <: AbstractArray{T,N}, BC <: Broadcasted}
    # we have an N dimensional loop.
    # need to construct the LoopSet
    loopsyms = [gensym(:n) for n ∈ 1:N]
    ls = LoopSet()
    pushpreamble!(ls, Expr(:(=), :dest, Expr(:call, :parent, :dest′)))
    sizes = Expr(:tuple)
    for (n,itersym) ∈ enumerate(loopsyms)
        Nsym = gensym(:N)
        ls.loops[itersym] = Loop(itersym, Nsym)
        push!(sizes.args, Nsym)
    end
    pushpreamble!(ls, Expr(:(=), sizes, Expr(:call, :size, :dest′)))
    elementbytes = sizeof(T)
    add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
    add_store!(ls, :dest, ArrayReference(:dest, reverse(loopsyms), Ref{Bool}(false)), elementbytes)
    resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
    q = lower(ls)
    push!(q.args, :dest′)
    pushfirst!(q.args, Expr(:meta,:inline))
    q
    # ls
end

@inline function vmaterialize(bc::Broadcasted)
    ElType = Base.Broadcast.combine_eltypes(bc.f, bc.args)
    vmaterialize!(similar(bc, ElType), bc)
end


text:
struct Gisette <: Name end

# mandatory methods
format(::Type{Gisette}) = TabularData
task(::Type{Gisette}) = TwoClass
source(::Type{Gisette}) = "http://archive.ics.uci.edu/ml/datasets/gisette"

function downloadlink(::Type{Gisette})
    return "https://archive.ics.uci.edu/ml/machine-learning-databases/gisette/" .* [
        "GISETTE/gisette_train.data",
        "GISETTE/gisette_train.labels",
        "GISETTE/gisette_valid.data",
        "gisette_valid.labels"
    ]
end

function preprocess(N::Type{Gisette})
    return [
        path -> preprocess_csv(N, path, :train; toremove = 5001)
        path -> preprocess_targets(N, path, :train)
        path -> preprocess_csv(N, path, :valid; toremove = 5001)
        path -> preprocess_targets(N, path, :valid)
    ]
end

nattributes(::Type{Gisette}) = (5000, )
ninstances(::Type{Gisette}) = (6000, 1000, 0)
positivelabel(::Type{Gisette}) = 1

# optional methods
function checksum(::Type{Gisette})
    return "d5a3e85f7db845564b066f01818a86b70bc2beefa2af3737abfc8a156b30a96d"
end

name(::Type{Gisette}) = "Gisette"
author(::Type{Gisette}) = ["Isabelle Guyon"]

function licence(::Type{Gisette})
    return """
    This data set is licensed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license. This allows for the sharing and adaptation of the datasets for any purpose, provided that the appropriate credit is given.

    URL: https://creativecommons.org/licenses/by/4.0/legalcode
    """
end

function citation(::Type{Gisette})
    return """
    @misc{misc_gisette_170,
        author       = {Guyon, Isabelle, Gunn, Steve, Ben-Hur, Asa & Dror, Gideon},
        title        = {{Gisette}},
        year         = {2008},
        howpublished = {UCI Machine Learning Repository}
    }
    """
end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
module TopOptTools
using Plots, LinearAlgebra, ForwardDiff, DataFrames
pyplot(framestyle = :box, legend = false, color = :black)

include("fem2d.jl")
include("opt.jl")

end


text:
  # wrap function -  using periodic boundary conditions to resize i and j of cells
  function wrap_cell(dms,id,jd)
    i = dms
    j = dms

    if id < 1
      while id < 1
        id = id + i
      end
    elseif id > i
      while id > i
        id = id - i
      end
    end

    if jd < 1
      while jd < 1
        jd = jd + j
      end
    elseif jd > j
      while jd > j
        jd = jd - j
      end
    end
    return id, jd
  end

  function wrapcell(nc,i,j)
    if i < 1
      i = nc + i
    elseif i > nc
      i = i - nc
    end
    if j < 1
      j = nc + j
    elseif j > nc
     j = j - nc
    end
    return i,j
  end


text:
using XYZFiles
using Test

@testset "read" begin
    allAtoms = XYZFiles.read("test.xyz", ["energy",])

    io = open("test.xyz", "r")
    allAtoms1 = XYZFiles.read(io, ["energy",])
    close(io)

    io = open("test.xyz", "r")
    allAtoms2 = XYZFiles.read(io)
    close(io)

    allAtoms3 = XYZFiles.read("test.xyz")

    XYZFiles.write("test_output.xyz", allAtoms1)
    @test length(allAtoms) == 3
    @test length(allAtoms1[1]) == 8
    @test length(allAtoms2[2]) == 6
    @test length(allAtoms3[3]) == 4
end

max_stars_repo_name: yaozhenghangma/XYZ

max_issues_repo_name: yaozhenghangma/XYZ

max_forks_repo_name: yaozhenghangma/XYZ


text:
import MPI
using Serialization

export MPIDistributor

"""
    MPIDistributor{GID, PID, LID}(comm::MPIComm{GID, PID, LID})
Creates an Distributor to work with MPIComm.  Created by
createDistributor(::MPIComm{GID, PID, LID})
"""
mutable struct MPIDistributor{GID <: Integer, PID <: Integer, LID <: Integer} <: Distributor{GID, PID, LID}
    comm::MPIComm{GID, PID, LID}

    lengths_to::Vector{GID}
    procs_to::Vector{PID}
    indices_to::Vector{GID}

    lengths_from::Vector{GID}
    procs_from::Vector{PID}
    indices_from::Vector{GID}

    resized::Bool
    sizes::Vector{GID}

    sizes_to::Vector{GID}
    starts_to::Vector{GID}
    #starts_to_ptr::Array{Integer}
    #indices_to_ptr::Array{Integer}

    sizes_from::Vector{GID}
    starts_from::Vector{GID}
    #sizes_from_ptr::Array{Integer}
    #starts_from_ptr::Array{Integer}

    numRecvs::GID
    numSends::GID
    numExports::GID

    selfMsg::GID

    maxSendLength::GID
    totalRecvLength::GID

    request::Vector{MPI.Request}
    status::Vector{MPI.Status}

    #sendArray::Array{UInt8}

    planReverse::Union{MPIDistributor{GID, PID, LID}, Nothing}

    importObjs::Union{Vector{<:Vector}, Nothing}

    datatype::DataType

    #never seem to be used
    #lastRoundBytesSend::Integer
    #lastRoundBytesRecv::Integer

    function MPIDistributor(comm::MPIComm{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
        new{GID, PID, LID}(comm, [], [], [], [], [], [], false, [], [], [], [], [],
            0, 0, 0, 0, 0, 0, [], [],  nothing,
            nothing, Nothing)
    end
end

#### internal methods ####
function createSendStructure(dist::MPIDistributor{GID, PID, LID}, pid::PID,
        nProcs::PID, exportPIDs::AbstractArray{PID}
        ) where {GID <: Integer, PID <: Integer, LID <: Integer}

    numExports = length(exportPIDs)
    dist.numExports = numExports

    starts = zeros(GID, nProcs+1)

    nactive = 0
    noSendBuff = true
    numDeadIndices::GID = 0  #for GIDs not owned by any processors

    for i = 1:numExports
        if noSendBuff && i > 1 && exportPIDs[i] < exportPIDs[i-1]
            noSendBuff = false
        end
        if exportPIDs[i] >= 1
            starts[exportPIDs[i]] += 1
            nactive += 1
        else
            numDeadIndices += 1
        end
    end

    dist.selfMsg = starts[pid] != 0
    dist.numSends = 0

    if noSendBuff #grouped by processor, no send buffer or indices_to needed
        for i = 1:nProcs
            if starts[i] > 0
                dist.numSends += 1
            end
        end

        dist.procs_to = Vector{PID}(undef, dist.numSends)
        dist.starts_to = Vector{GID}(undef, dist.numSends)
        dist.lengths_to = Vector{GID}(undef, dist.numSends)

        index = numDeadIndices+1
        for i = 1:dist.numSends
            dist.starts_to[i] = index
            proc = exportPIDs[index]
            dist.procs_to[i] = proc
            index += starts[proc]
        end

        perm = sortperm(dist.procs_to)
        dist.procs_to = dist.procs_to[perm]
        dist.starts_to = dist.starts_to[perm]

        # line 430

        dist.maxSendLength = 0

        for i = 1:dist.numSends
            proc = dist.procs_to[i]
            dist.lengths_to[i] = starts[proc]
            if (proc != pid) && (dist.lengths_to[i] > dist.maxSendLength)
                maxSendLength = dist.lengths_to[i]
            end
        end
    else #not grouped by processor, need send buffer and indices_to
        if starts[1] != 0
            dist.numSends = 1
        end

        for i = 2:nProcs
            if starts[i] != 0
                dist.numSends += 1
            end
            starts[i] += starts[i-1]
        end

        for i = nProcs:-1:2
            starts[i] = starts[i-1] + 1
        end
        starts[1] = 1

        if nactive > 0
            dist.indices_to = Array{GID, 1}(undef, nactive)
        end

        for i = 1:numExports
            if exportPIDs[i] >= 1
                dist.indices_to[starts[exportPIDs[i]]] = i
                starts[exportPIDs[i]] += 1
            end
        end

        #reconstruct starts array to index into indices_to

        for i = nProcs:-1:2
            starts[i] = starts[i-1]
        end
        starts[1] = 1
        starts[nProcs+1] = nactive+1

        if dist.numSends > 0
            dist.lengths_to = Vector{GID}(undef, dist.numSends)
            dist.procs_to = Vector{PID}(undef, dist.numSends)
            dist.starts_to = Vector{GID}(undef, dist.numSends)
        end

        j::GID = 1
        dist.maxSendLength = 0
        for i = 1:nProcs
            if starts[i+1] != starts[i]
                dist.lengths_to[j] = starts[i+1] - starts[i]
                dist.starts_to[j] = starts[i]
                if (i != pid) && (dist.lengths_to[j] > dist.maxSendLength)
                    dist.maxSendLength = dist.lengths_to[j]
                end
                dist.procs_to[j] = i
                j += 1
            end
        end
    end

    dist.numSends -= dist.selfMsg

    dist
end

function computeRecvs(dist::MPIDistributor{GID, PID, LID}, myProc::PID, nProcs::PID) where GID <: Integer where PID <: Integer where LID <: Integer

    msgCount = zeros(Int, nProcs)

    for proc in dist.procs_to
        msgCount[proc] += 1
    end

    #bug fix for reduce-scatter bug applied since no reduce_scatter is present in julia's MPI
    rawCounts = MPI.Reduce(msgCount, MPI.SUM, 0, dist.comm.mpiComm)
    if rawCounts isa Nothing
        counts = Int[]
    else
        counts = rawCounts
    end
    totalRecvs = MPI.Scatter(counts, 1, 0, dist.comm.mpiComm)[1]

    dist.lengths_from = zeros(Int, totalRecvs)
    dist.procs_from = zeros(PID, totalRecvs)

    #using NEW_COMM_PATTERN (see line 590)

    if dist.request == []
        dist.request = Vector{MPI.Request}(undef, totalRecvs - dist.selfMsg)
    end

    #line 616

    lengthWrappers = [Array{Int, 1}(undef, 1) for i in 1:(totalRecvs - dist.selfMsg)]
    for i = 1:(totalRecvs - dist.selfMsg)
        dist.request[i] = MPI.Irecv!(lengthWrappers[i], MPI.ANY_SOURCE, MPI.ANY_TAG, dist.comm.mpiComm)
    end

    barrier(dist.comm)

    for i = 1:(dist.numSends+dist.selfMsg)
        if dist.procs_to[i] != myProc
            #have to use Rsend in MPIUtil
            MPI_Rsend(dist.lengths_to[i], dist.procs_to[i]-1, 1, dist.comm.mpiComm)
        else
            dist.lengths_from[totalRecvs] = dist.lengths_to[i]
            dist.procs_from[totalRecvs] = myProc
        end
    end

    if totalRecvs > dist.selfMsg
        dist.status = MPI.Waitall!(dist.request)
    end

    for i = 1:(totalRecvs - dist.selfMsg)
        dist.lengths_from[i] = lengthWrappers[i][1]
    end

    for i = 1:(totalRecvs - dist.selfMsg)
        dist.procs_from[i] = MPI.Get_source(dist.status[i])+1
    end

    perm = sortperm(dist.procs_from)
    dist.procs_from = dist.procs_from[perm]
    dist.lengths_from = dist.lengths_from[perm]

    dist.starts_from = Vector{GID}(undef, totalRecvs)
    j = GID(1)
    for i = 1:totalRecvs
        dist.starts_from[i] = j
        j += dist.lengths_from[i]
    end

    dist.totalRecvLength = 0
    for i = 1:totalRecvs
        dist.totalRecvLength += dist.lengths_from[i]
    end

    dist.numRecvs = totalRecvs - dist.selfMsg

    dist
end

function computeSends(dist::MPIDistributor{GID, PID, LID},
        remoteGIDs::AbstractArray{GID, 1}, remotePIDs::AbstractArray{PID, 1}
        )::Tuple{AbstractArray{GID, 1}, AbstractArray{PID, 1}
        } where {GID <:Integer, PID <:Integer, LID <:Integer}
    numImports = length(remoteGIDs)

    tmpPlan = MPIDistributor(dist.comm)

    importObjs = Vector{Tuple{GID, PID}}(undef, numImports)
    for i = 1:numImports
        importObjs[i] = (remoteGIDs[i], myPid(dist.comm))#remotePIDs[i])
    end

    numExports = createFromSends(tmpPlan, copy(remotePIDs))

    exportIDs = Vector{GID}(undef, numExports)
    exportProcs = Vector{PID}(undef, numExports)

    exportObjs = resolve(tmpPlan, importObjs)
    for i = 1:numExports
        exportIDs[i] = exportObjs[i][1]
        exportProcs[i] = exportObjs[i][2]
    end
    (exportIDs, exportProcs)
end

"""
Creates a reverse distributor for the given MPIDistributor
"""
function createReverseDistributor(dist::MPIDistributor{GID, PID, LID}
        ) where {GID <: Integer, PID <: Integer, LID <: Integer}
    myProc = myPid(dist.comm)

    if dist.planReverse == nothing
        totalSendLength = reduce(+, dist.lengths_to)

        maxRecvLength::GID = 0
        for i = 1:dist.numRecvs
            if dist.procs_from[i] != myProc
                maxRecvLength = max(maxRecvLength, dist.lengths_from[i])
            end
        end

        reverse = MPIDistributor(dist.comm)
        dist.planReverse = reverse

        reverse.lengths_to = dist.lengths_from
        reverse.procs_to = dist.procs_from
        reverse.indices_to = dist.indices_from
        reverse.starts_to = dist.starts_from

        reverse.lengths_from = dist.lengths_to
        reverse.procs_from = dist.procs_to
        reverse.indices_from = dist.indices_to
        reverse.starts_from = dist.starts_to

        reverse.numSends = dist.numRecvs
        reverse.numRecvs = dist.numSends
        reverse.selfMsg = dist.selfMsg

        reverse.maxSendLength = maxRecvLength
        reverse.totalRecvLength = totalSendLength

        reverse.request = Vector{MPI.Request}(undef, reverse.numRecvs)
        reverse.status  = Vector{MPI.Status}(undef, reverse.numRecvs)
    end

    nothing
end

#### Internal Utilities ####
"""
    arraytobytes(value::Vector{T})::Vector{UInt8}

Converts the given vector to a list of bytes.
"""
function arraytobytes(value::AbstractVector{T}) where T
    if isbitstype(T)
        # if T is a bits type, then just use the buffer as it
        sz = sizeof(value)
        bytes_ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, value))
        unsafe_wrap(Array, bytes_ptr, sz)
    else
        # otherwise write using an IOBuffer
        buffer = IOBuffer(;sizehint=sizeof(value))
        write(buffer, value)
        take!(buffer)
    end
end

function arraytobytes(value::AbstractVector{<:AbstractVector})
    buffer = IOBuffer()
    for v in value
        write(buffer, UInt32(length(v)))
        for elt in v
            write(buffer, elt)
        end
    end
    take!(buffer)
end

"""
    bytestoarray(bytes::Vector{UInt8}, ::Type{T})::Vector{T} where T

Converts a list of bytes to a vector with contents of the given type
"""
function bytestoarray(bytes::AbstractVector{UInt8}, ::Type{T}) where T
    sz = fld(length(bytes), sizeof(T))
    value = Vector{T}(undef, sz)
    if isbitstype(T)
        bytes_ptr = convert(Ptr{T}, Base.unsafe_convert(Ptr{UInt8}, bytes))
        value_ptr = Base.unsafe_convert(Ptr{T}, value)
        Base.unsafe_copyto!(value_ptr, bytes_ptr, sz)
    else
        buffer = IOBuffer(bytes)
        i = 1
        while !eof(buffer)
            value[i] = read(buffer, T)
            i+=1
        end
    end
    value
end

function bytestoarray(bytes::AbstractVector{UInt8}, ::Type{Vector{T}}) where T
    buffer = IOBuffer(bytes)
    value = Vector{Vector{T}}()
    while !eof(buffer)
        len = read(buffer, UInt32)
        elt = Vector{T}(undef, len)
        for i in 1:len
            elt[i] = read(buffer, T)
        end
        push!(value, elt)
    end
    value
end

#### General Interface ####

getComm(dist::MPIDistributor) = dist.comm

#### Distributor interface ####

function createFromSends(dist::MPIDistributor{GID, PID, LID}, exportPIDs::AbstractArray{PID, 1})::Integer where GID <:Integer where PID <:Integer where LID <:Integer
    pid = myPid(dist.comm)
    nProcs = numProc(dist.comm)
    createSendStructure(dist, pid, nProcs, exportPIDs)
    computeRecvs(dist, pid, nProcs)
    if dist.numRecvs > 0
        if dist.request == []
            dist.request = Vector{MPI.Request}(undef, dist.numRecvs)
            dist.status = Vector{MPI.Status}(undef, dist.numRecvs)
        end
    end
    dist.totalRecvLength
end

function createFromRecvs(dist::MPIDistributor{GID, PID, LID},
        remoteGIDs::AbstractArray{GID, 1}, remotePIDs::AbstractArray{PID, 1}
        )::Tuple{AbstractArray{GID, 1}, AbstractArray{PID, 1}
        } where {GID <: Integer, PID <: Integer, LID <: Integer}
    if length(remoteGIDs) != length(remotePIDs)
        throw(InvalidArgumentError("remote lists must be the same length"))
    end
    (exportGIDs, exportPIDs) = computeSends(dist, remoteGIDs, remotePIDs)
    createFromSends(dist, exportPIDs)

    exportGIDs, exportPIDs
end

function resolvePosts(dist::MPIDistributor{GID, PID, LID}, exportObjs::AbstractArray{T, 1}) where {T, GID<:Integer, PID<:Integer, LID<:Integer}
    myProc = myPid(dist.comm)

    dist.datatype = T

    nBlocks::GID = dist.numSends + dist.selfMsg
    procIndex::GID = 1
    while procIndex <= nBlocks && dist.procs_to[procIndex] < myProc
        procIndex += 1
    end
    if procIndex == nBlocks
        procIndex = 1
    end

    exportBytes = Vector{Vector{UInt8}}(undef, dist.numRecvs + dist.selfMsg)

    j::GID = 1

    if length(dist.indices_to) == 0 #data already grouped by processor
        for i = 1:nBlocks
            # note that the range needs to be signed ints because Julia doesn't treat unsigened ints as first class
            exportBytes[i] = arraytobytes(view(exportObjs, Int64(j):Int64(j+dist.lengths_to[i]-1)))
            j += dist.lengths_to[i]
        end
    else #data not grouped by proc, must be grouped first
        for i = 1:nBlocks
            j = dist.starts_to[i]
            sendArray = Array{T}(undef, dist.lengths_to[i])
            for k = 1:dist.lengths_to[i]
                sendArray[k] = exportObjs[dist.indices_to[j+k-1]]
            end
            exportBytes[i] = arraytobytes(sendArray)
        end
    end

    ## get sizes of data begin received ##
    lengthRequests = Vector{MPI.Request}(undef, dist.numRecvs)
    lengths = Vector{Vector{Int}}(undef, dist.numRecvs + dist.selfMsg)
    for i = 1:dist.numRecvs + dist.selfMsg
        lengths[i] = Vector{Int}(undef, 1)
    end

    j = 1

    for i = 1:dist.numRecvs + dist.selfMsg
        if dist.procs_from[i] != myProc
            lengthRequests[j] = MPI.Irecv!(lengths[i], dist.procs_from[i]-1, MPI.ANY_TAG, dist.comm.mpiComm)
            j += 1
        end
    end

    barrier(dist.comm)

    for i = 1:dist.numSends + dist.selfMsg
        p = procIndex + 1
        if p > nBlocks
            p -= nBlocks
        end
        if dist.procs_to[i] != myProc
            MPI_Rsend(length(exportBytes[i]), dist.procs_to[i]-1, 2, dist.comm.mpiComm)
        else
            lengths[i][1] = length(exportBytes[i])
        end
    end

    MPI.Waitall!(lengthRequests)
    #at this point `lengths` should contain the sizes of incoming data

    importObjs = Vector{Vector{UInt8}}(undef, dist.numRecvs+dist.selfMsg)
    for i = 1:length(importObjs)
        importObjs[i] = Vector{UInt8}(undef, lengths[i][1])
    end

    dist.importObjs = importObjs

    ## back to the regularly scheduled program ##

    k::GID = 0
    j = 0
    selfRecvAddress::GID = 0
    for i = 1:dist.numRecvs + dist.selfMsg
        if dist.procs_from[i] != myProc
            MPI.Irecv!(importObjs[i], dist.procs_from[i]-1, MPI.ANY_TAG, dist.comm.mpiComm)
        else
            selfRecvAddress = i
        end
    end

    barrier(dist.comm)

    #line 844

    selfNum::GID = 0

#    if dist.indices_to == [] #data already grouped by processor
    for i = 1:nBlocks
        p = i + procIndex - 1
        if p > nBlocks
            p -= nBlocks
        end
        if dist.procs_to[p] != myProc
            MPI_Rsend(exportBytes[p], dist.procs_to[p]-1, 3, dist.comm.mpiComm)
        else
            selfNum = p
        end
    end

    if dist.selfMsg != 0
        importObjs[selfRecvAddress] = exportBytes[selfNum]
    end

    nothing
end

function resolveWaits(dist::MPIDistributor)::Array
    barrier(dist.comm)#run into issues deserializing otherwise
    if dist.numRecvs> 0
        dist.status = MPI.Waitall!(dist.request)
    end

    if dist.importObjs == nothing
        throw(InvalidStateError("Cannot resolve waits when no posts have been made"))
    end

    T = dist.datatype
    barrier(dist.comm)#run into issues deserializing otherwise
    importObjs = dist.importObjs::Vector{Vector{UInt8}}
    deserializedObjs = Vector{Vector{T}}(undef, length(importObjs))
    for i = 1:length(importObjs)
        deserializedObjs[i] = bytestoarray(importObjs[i], T)
    end
    results = deserializedObjs

    dist.importObjs = nothing

    reduce(vcat, results; init = T[])
end

function resolveReversePosts(dist::MPIDistributor{GID, PID, LID},
        exportObjs::AbstractVector{T}
        ) where {GID <: Integer, PID <: Integer, LID <: Integer, T}
    if dist.indices_to != []
        throw(InvalidStateError("Cannot do reverse comm when data is not blocked by processor"))
    end

    if dist.planReverse == nothing
        createReverseDistributor(dist)
    end

    resolvePosts(dist.planReverse, exportObjs)
end

function resolveReverseWaits(dist::MPIDistributor{GID, PID, LID}
    )::Vector where {GID <: Integer, PID <: Integer, LID <: Integer}
    if dist.planReverse == nothing
        throw(InvalidStateError("Cannot resolve reverse waits if there is no reverse plan"))
    end

    resolveWaits(dist.planReverse)
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   SOLFSMY.txt
#
#   This file contains the following indices:
#
#       F10, F81c, S10, S81c, M10, M81c, Y10, Y81c
#
#   in which 81c means the 81-day averaged centered value.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

################################################################################
#                       Private Structures and Variables
################################################################################

"""
    _SOLFSMY_Structure

Structure to store the interpolations of the data in `SOLFSMY.TXT` file.

# Fields

* `F10`: 10.7-cm solar flux [10⁻²² W/(m² Hz)].
* `F81a`: 10.7-cm averaged solar flux, 81-day centered on input time.
* `S10`: EUV index.
* `S81a`: EUV 81-day averaged centered index.
* `M10`: MG2 index scaled to F10.
* `M81a`: MG2 81-day averaged centered index.
* `Y81a`: Solar X-ray & Lya 81-day averaged centered index.
* `Y81a`: Solar X-ray & Lya 81-day averaged centered index.

"""
struct _SOLFSMY_Structure
    F10::_space_indices_itp_constant{Float64}
    F81a::_space_indices_itp_constant{Float64}
    S10::_space_indices_itp_constant{Float64}
    S81a::_space_indices_itp_constant{Float64}
    M10::_space_indices_itp_constant{Float64}
    M81a::_space_indices_itp_constant{Float64}
    Y10::_space_indices_itp_constant{Float64}
    Y81a::_space_indices_itp_constant{Float64}
end

# Remote file: SOLFSMY.TXT
_solfsmy = @RemoteFile(
    "http://sol.spacenvironment.net/jb2008/indices/SOLFSMY.TXT",
    file="SOLFSMY.TXT",
    updates=:daily
   )

# Optional variable that will store the `SOLFSMY.TXT` data.
@OptionalData(_solfsmy_data, _SOLFSMY_Structure,
              "Run `init_space_indices()` with `:solfsmy` in `enabled_files` array to initialize required data.")

################################################################################
#                               Public Functions
################################################################################

#                                   Getters
# ==============================================================================

# Create a function `get_S` for each `S` field in `_SOLFSMY_Structure`.
for field in fieldnames(_SOLFSMY_Structure)
    field_name = string(field)
    func_name  = "get_" * field_name
    func       = Symbol("get_" * string(field))
    qfield     = Meta.quot(field)

    @eval begin
        export $func
"""
    $($func_name)(JD::Number)

Get the value of the index `$($field_name)` at Julian Day `JD`.

This function requires the initialization of the variable `_solfsmy_data`.
Otherwise, an exception will be raised. To initialize it, run
`init_space_indices()`.

"""
        ($func)(JD::Number) = getfield(get(_solfsmy_data), $qfield)(JD)
    end
end

################################################################################
#                              Private Functions
################################################################################

"""
    _init_solfsmy(;force_download = false, local_path = nothing)

Initialize the data in the file `SOLFSMY.TXT` by creating `_solfsmy_data`. The
initialization process is composed of:

1. Download the file, if it is necessary;
2. Parse the file;
3. Create the interpolations and the structures.

If the keyword `force_download` is `true`, then the file will always be
downloaded.

The user can also specify a location for the file using the keyword
`local_path`. If it is `nothing`, which is the default, then the file will be
downloaded.

"""
function _init_solfsmy(;force_download = false, local_path = nothing)
    # Update the remote files if no path is given.
    if local_path == nothing
        download(_solfsmy; force = force_download)
        local_path = path(_solfsmy)
    end

    push!(_solfsmy_data,   _parse_solfsmy(local_path))

    nothing
end

"""
    _parse_solfsmy(path::AbstractString)

Parse the `SOLFSMY.TXT` file in `path` and retur an instance of the structure
`_SOLFSMY_Structure` with the initialized interpolations.

The format of the file `SOLFSMY.TXT` must be:

    YYYY DDD   JulianDay  F10   F81c  S10   S81c  M10   M81c  Y10   Y81c  Ssrc

"""
function _parse_solfsmy(path::AbstractString)
    # Allocate the raw data.
    JD   = Float64[]
    F10  = Float64[]
    F81a = Float64[]
    S10  = Float64[]
    S81a = Float64[]
    M10  = Float64[]
    M81a = Float64[]
    Y10  = Float64[]
    Y81a = Float64[]

    # Read the raw data in the file.
    open(path) do file
        line_num = 0

        for ln in eachline(file)
            line_num += 1

            # Ignore comments.
            (ln[1] == '#') && continue

            tokens   = split(ln)
            num_toks = length(tokens)

            # Ignore blank lines.
            (num_toks == 0) && continue

            # Ignore and print warning about lines with bad format.
            if num_toks != 12
                @warn "Line $line_num of file SOLFSMY.TXT has invalid format!"
                continue
            end

            # Parse data.
            push!(JD  , parse(Float64, tokens[ 3]))
            push!(F10 , parse(Float64, tokens[ 4]))
            push!(F81a, parse(Float64, tokens[ 5]))
            push!(S10 , parse(Float64, tokens[ 6]))
            push!(S81a, parse(Float64, tokens[ 7]))
            push!(M10 , parse(Float64, tokens[ 8]))
            push!(M81a, parse(Float64, tokens[ 9]))
            push!(Y10 , parse(Float64, tokens[10]))
            push!(Y81a, parse(Float64, tokens[11]))
        end
    end

    # Create the interpolations for each parameter.
    #
    # The parameters in `SOLFSMY.TXT` are computed at 12 UT. Hence, if we
    # interpolate by the nearest-neighbor, then it will always return the data
    # related to that day.
    knots    = (JD,)
    itp_F10  = interpolate(knots, F10 , Gridded(Constant()))
    itp_F81a = interpolate(knots, F81a, Gridded(Constant()))
    itp_S10  = interpolate(knots, S10 , Gridded(Constant()))
    itp_S81a = interpolate(knots, S81a, Gridded(Constant()))
    itp_M10  = interpolate(knots, M10 , Gridded(Constant()))
    itp_M81a = interpolate(knots, M81a, Gridded(Constant()))
    itp_Y10  = interpolate(knots, Y10 , Gridded(Constant()))
    itp_Y81a = interpolate(knots, Y81a, Gridded(Constant()))

    _SOLFSMY_Structure(itp_F10, itp_F81a,
                       itp_S10, itp_S81a,
                       itp_M10, itp_M81a,
                       itp_Y10, itp_Y81a)
end


text:
module GeometricEquations

    using GeometricBase

    import Base: Callable
    import GeometricBase: invariants, parameters, periodicity

    export NullInvariants, NullParameters, NullPeriodicity

    export OptionalAbstractArray, OptionalArray,
           OptionalFunction, OptionalNamedTuple,
           OptionalInvariants, OptionalParameters,
           OptionalPeriodicity

    export State, StateVector

    export GeometricEquation
    export AbstractEquationODE, AbstractEquationPODE,
           AbstractEquationDAE, AbstractEquationPDAE,
           AbstractEquationSDE, AbstractEquationPSDE

    export AbstractProblem, GeometricProblem

    export ODE, IODE, PODE, HODE, LODE, SODE
    export DAE, IDAE, PDAE, HDAE, LDAE, SPDAE
    export SDE, PSDE, SPSDE

    export ODEProblem,  IODEProblem, PODEProblem,
           HODEProblem, LODEProblem, SODEProblem
    export DAEProblem,  IDAEProblem, PDAEProblem,
           HDAEProblem, LDAEProblem, SPDAEProblem
    export SDEProblem,  PSDEProblem, SPSDEProblem

    export ODEEnsemble,  IODEEnsemble, PODEEnsemble,
           HODEEnsemble, LODEEnsemble, SODEEnsemble
    export DAEEnsemble,  IDAEEnsemble, PDAEEnsemble,
           HDAEEnsemble, LDAEEnsemble, SPDAEEnsemble
    export SDEEnsemble,  PSDEEnsemble, SPSDEEnsemble

    export nsamples, nconstraints
    export equation, functions, solutions
    export invariants, parameters, periodicity
    export initial_conditions

    export hassolution, hasvectorfield, hasprimary, hassecondary,
           hasinvariants, hasparameters, hasperiodicity,
           hashamiltonian, haslagrangian

    export tspan, tstep, tbegin, tend

    include("utils.jl")

    include("geometric_equation.jl")

    include("ode.jl")
    include("hode.jl")
    include("iode.jl")
    include("lode.jl")
    include("pode.jl")

    include("dae.jl")
    include("hdae.jl")
    include("idae.jl")
    include("ldae.jl")
    include("pdae.jl")

    include("sode.jl")
    include("spdae.jl")

    include("sde.jl")
    include("psde.jl")
    include("spsde.jl")

    include("conversion.jl")

    include("abstract_problem.jl")
    include("geometric_problem.jl")
end