text:
using Test
using LazyJSON

@testset "SplicedString" begin

sv = LazyJSON.SplicedStrings.splice_vector
SS = LazyJSON.SplicedStrings.SplicedString

FI(n, i) = n << LazyJSON.SplicedStrings.fragment_bits | i + LazyJSON.SplicedStrings.index_offset

@test SS("Foo").v == ["Foo"]
@test SS(SubString("Food", 1, 3)).v == ["Foo"]
@test SS(SubString(SS("Food"), 1, FI(0,3))).v == ["Foo"]
@test SS(SubString(SS("Food"), 1, FI(0,3))).v[1] isa SubString{String}
@test SS("Foo", "Bar").v == ["Foo", "Bar"]
@test SS(SS("Foo", "Bar")).v == ["Foo", "Bar"]

@test sv(         "Foo"                              ) == ["Foo"]
@test sv(   [     "Foo",       "Bar"                ]) == ["Foo", "Bar"]
@test sv(   (     "Foo",       "Bar"                )) == ["Foo", "Bar"]
@test sv(     SS(["Foo",       "Bar"])               ) == ["Foo", "Bar"]
@test sv(   [ SS(["Foo"]), SS(["Bar"])              ]) == ["Foo", "Bar"]
@test sv(   [ SS( "Foo",       "Bar" ), SS( "Boo" ) ]) == ["Foo", "Bar", "Boo"]
@test sv(   [ SS(["Foo",       "Bar"]), SS( "Boo" ) ]) == ["Foo", "Bar", "Boo"]
@test sv(   [ SS(["Foo",       "Bar"]), SS(["Boo"]) ]) == ["Foo", "Bar", "Boo"]
@test sv(   [ SS(["Foo",       "Bar"]),     "Boo"   ]) == ["Foo", "Bar", "Boo"]
@test sv(   [     "Foo",   SS(["Bar",       "Boo"]) ]) == ["Foo", "Bar", "Boo"]

@test SS(         "Foo"                              ) == "Foo"
@test SS(   [     "Foo"                             ]) == "Foo"
@test SS(         "Foo",       "Bar"                 ) == "FooBar"
@test SS(   [     "Foo",       "Bar"                ]) == "FooBar"
@test SS(     SS(["Foo",       "Bar"])               ) == "FooBar"
@test SS(   [ SS(["Foo",       "Bar"])              ]) == "FooBar"
@test SS(     SS(["Foo"]), SS(["Bar"])               ) == "FooBar"
@test SS(   [ SS(["Foo"]), SS(["Bar"])              ]) == "FooBar"
@test SS(     SS(["Foo",       "Bar"]), SS(["Boo"])  ) == "FooBarBoo"
@test SS(         "Foo",       "Bar",       "Boo"    ) == "FooBarBoo"
@test SS(     SS(["Foo",       "Bar"]),     "Boo"    ) == "FooBarBoo"
@test SS(   [ SS(["Foo",       "Bar"]),     "Boo"   ]) == "FooBarBoo"
@test SS(         "Foo",   SS(["Bar",       "Boo"])  ) == "FooBarBoo"
@test SS(   [     "Foo",   SS(["Bar",       "Boo"]) ]) == "FooBarBoo"
@test SS(   [     "Foo",   SS(["Bar",       "Boo"]) ]) == "FooBarBoo"

# Simple insert
ss = SS("Foo", "Bar")
splice!(ss, FI(1,1), FI(1,0), "d")
@test ss == "FoodBar"
@test ss.v == ["Foo", "d", "Bar"]

# Replace fragment prefix
ss = SS("Foo", "Bar")
splice!(ss, FI(1,1), FI(1,1), "Be")
@test ss == "FooBear"
@test ss.v == ["Foo", "Be", "ar"]

# Replace whole fragment
splice!(ss, FI(1,1), FI(1,2), "C")
@test ss == "FooCar"
@test ss.v == ["Foo", "C", "ar"]

# Replace two whole fragments
ss = SS("Foo", "XXX", "YYY", "Bar")
splice!(ss, FI(1,1), FI(2,3), "_")
@test ss == "Foo_Bar"
@test ss.v == ["Foo", "_", "Bar"]

# Replace three whole fragments
ss = SS("Foo", "XXX", "YYY", "Bar")
splice!(ss, FI(1,1), FI(3,3), "_")
@test ss == "Foo_"
@test ss.v == ["Foo", "_"]

# Replace all fragments
ss = SS("Foo", "XXX", "YYY", "Bar")
splice!(ss, 1 : FI(3,3), "_")
@test ss == "_"
@test ss.v == ["_"]

# Replace all fragments with nothing
ss = SS("Foo", "XXX", "YYY", "Bar")
splice!(ss, firstindex(ss), lastindex(ss), "")
@test ss == ""
@test ss.v == []

ss = SS("AAA", "BBB", "CCC")
i = findfirst(isequal('B'), ss)
splice!(ss, i , i+2, "")
@test ss == "AAACCC"
@test ss.v == ["AAA", "CCC"]

# Replace two whole fragments and part of another
ss = SS("Foo", "XXX", "YYY", "Bar")
splice!(ss, FI(1,1), FI(3,2), "u")
@test ss == "Foour"
@test ss.v == ["Foo", "u", "r"]

# Replace fragment suffix
ss = SS("Foo", "Bar")
splice!(ss, FI(0,3), FI(0,3), "g")
@test ss == "FogBar"
@test ss.v == ["Fo", "g", "Bar"]

# Replace fragment suffix and prefix
ss = SS("Foo", "Bar")
splice!(ss, FI(0,3), FI(1,2), "lde")
@test ss == "Folder"
@test ss.v == ["Fo", "lde", "r"]

# Replace fragment suffix and whole fragment and prefix
ss = SS("Foo", "XXX", "Bar")
splice!(ss, FI(0,3), FI(2,2), "lde")
@test ss == "Folder"
@test ss.v == ["Fo", "lde", "r"]

# Replace fragment suffix and whole fragment and prefix with multiple others
ss = SS("Foo", "XXX", "Bar")
splice!(ss, FI(0,2), FI(2,2), ["w", "oooo", "aaaaaa", "rr"])
@test ss == "Fwooooaaaaaarrr"
@test ss.v == ["F", "w", "oooo", "aaaaaa", "rr", "r"]

ss = SS("Foo", "XXX", "Bar")
splice!(ss, FI(0,2), FI(2,2), SS(SS("w", "oooo"), "aaaaaa", "rr"))
@test ss == "Fwooooaaaaaarrr"
@test ss.v == ["F", "w", "oooo", "aaaaaa", "rr", "r"]

ss = SS("Foo", "XXX", "Bar")
splice!(ss, FI(0,2), FI(2,2), [SS("w", "oooo"), "aaaaaa", SS("rr")])
@test ss == "Fwooooaaaaaarrr"
@test String(ss) == "Fwooooaaaaaarrr"
@test ss.v == ["F", "w", "oooo", "aaaaaa", "rr", "r"]

ss = SS(["Hello", " ", "world", "!"])
@test ss == "Hello world!"
@test String(ss) == "Hello world!"

ss = SS("xxxx", "Hello", " ", "world", "!", "yyyy")
@test SubString(ss, findfirst(isequal('l'), ss),
                    findfirst(isequal('r'), ss)) == "llo wor"

sss = SubString(ss, findfirst(isequal('l'), ss),
                    findfirst(isequal('r'), ss))
@test sss isa SubString{SS}

ss = SS("AAA", sss, "BBB")
@test ss.v == ["AAA", "llo", " ", "wor", "BBB"]

ss = SS(["Hello", " ", "world", "!"])
@test ss == "Hello world!"

i =      0  ; @test thisind(ss, i) == i     ; @test !isvalid(ss, i) ; @test nextind(ss, i) == 1
i =      1  ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == FI(0,2) ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('H')
i = FI(0,2) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1   ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('e')
i = FI(0,3) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1   ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('l')
i = FI(0,4) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1   ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('l')
i = FI(0,5) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == FI(1,1) ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('o')
i = FI(1,1) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == FI(2,1) ; @test prevind(ss, i) == FI(0,5) ; @test codeunit(ss, i) == UInt8(' ')
i = FI(2,1) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1   ; @test prevind(ss, i) == FI(1,1) ; @test codeunit(ss, i) == UInt8('w')
i = FI(2,2) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1   ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('o')
i = FI(2,3) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1   ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('r')
i = FI(2,4) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1   ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('l')
i = FI(2,5) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == FI(3,1) ; @test prevind(ss, i) == i - 1   ; @test codeunit(ss, i) == UInt8('d')
i = FI(3,1) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1   ; @test prevind(ss, i) == FI(2,5) ; @test codeunit(ss, i) == UInt8('!')
i = FI(3,2) ; @test thisind(ss, i) == i     ; @test !isvalid(ss, i)                                   ; @test prevind(ss, i) == i - 1
              @test ncodeunits(ss) == FI(3,1)
              @test length(ss) == 12
              @test length(ss,      1 , FI(3,1)) == 12
              @test length(ss, FI(0,2), FI(2,5)) == 10
              @test length(ss, FI(0,3), FI(2,4)) == 8
              @test length(ss, FI(0,4), FI(2,3)) == 6
              @test length(ss, FI(0,5), FI(2,2)) == 4
              @test length(ss, FI(1,1), FI(2,1)) == 2
              @test length(ss, FI(1,1), FI(1,1)) == 1
              @test length(ss, FI(1,1), FI(0,5)) == 0

ss = SS(["\u1234x", "x\u1234"])
@test ss == "ሴxxሴ"

i =      0  ; @test thisind(ss, i) == i     ; @test !isvalid(ss, i) ; @test nextind(ss, i) == 1
i =      1  ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == FI(0,2) + 2
i = FI(0,2) ; @test thisind(ss, i) == i - 1 ; @test !isvalid(ss, i) ; @test nextind(ss, i) == i + 2
i = FI(0,3) ; @test thisind(ss, i) == i - 2 ; @test !isvalid(ss, i) ; @test nextind(ss, i) == i + 1
i = FI(0,4) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == FI(1,1)
i = FI(1,1) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1
i = FI(1,2) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 3
i = FI(1,3) ; @test thisind(ss, i) == i - 1 ; @test !isvalid(ss, i) ; @test nextind(ss, i) == i + 2
i = FI(1,4) ; @test thisind(ss, i) == i - 2 ; @test !isvalid(ss, i) ; @test nextind(ss, i) == i + 1
i = FI(1,5) ; @test thisind(ss, i) == i     ; @test !isvalid(ss, i)
              @test ncodeunits(ss) == FI(1,4)
              @test length(ss) == 4

ss = SS(["x\u1234", "\u1234x"])
@test ss == "xሴሴx"

i =      0  ; @test thisind(ss, i) == i     ; @test !isvalid(ss, i) ; @test nextind(ss, i) == 1
i =      1  ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == FI(0,2)
i = FI(0,2) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == FI(1,1)
i = FI(0,3) ; @test thisind(ss, i) == i - 1 ; @test !isvalid(ss, i) ; @test nextind(ss, i) == FI(1,1)
i = FI(0,4) ; @test thisind(ss, i) == i - 2 ; @test !isvalid(ss, i) ; @test nextind(ss, i) == FI(1,1)
i = FI(1,1) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 3
i = FI(1,2) ; @test thisind(ss, i) == i - 1 ; @test !isvalid(ss, i) ; @test nextind(ss, i) == i + 2
i = FI(1,3) ; @test thisind(ss, i) == i - 2 ; @test !isvalid(ss, i) ; @test nextind(ss, i) == i + 1
i = FI(1,4) ; @test thisind(ss, i) == i     ; @test  isvalid(ss, i) ; @test nextind(ss, i) == i + 1
i = FI(1,5) ; @test thisind(ss, i) == i     ; @test !isvalid(ss, i)
              @test ncodeunits(ss) == FI(1,4)
              @test length(ss) == 4

ss = SS("one", "two", "three")
@test 'o' in ss
@test 'n' in ss
@test 'e' in ss
@test 't' in ss
@test 'w' in ss
@test 'o' in ss
@test 'h' in ss
@test 'r' in ss
i = findfirst(isequal('w'), ss)
@test ss[i] == 'w'

#=
ss = SS("one", "two", "three")
cu = LazyJSON.SplicedStrings.densecodeunits(ss)

@test  cu[1] == UInt8('o')
@test  cu[2] == UInt8('n')
@test  cu[3] == UInt8('e')
@test  cu[4] == UInt8('t')
@test  cu[5] == UInt8('w')
@test  cu[6] == UInt8('o')
@test  cu[7] == UInt8('t')
@test  cu[8] == UInt8('h')
@test  cu[9] == UInt8('r')
@test cu[10] == UInt8('e')
@test cu[11] == UInt8('e')
=#

ncu = LazyJSON.SplicedStrings.nextcodeunit

@test ncu(ss,         0) == (     1,  UInt8('o'))
@test ncu(ss,        1 ) == (FI(0,2), UInt8('n'))
@test ncu(ss,   FI(0,2)) == (FI(0,3), UInt8('e'))
@test ncu(ss,   FI(0,3)) == (FI(1,1), UInt8('t'))

end # @testset "SplicedString"


text:
function velocity(H,Br,m,y,lat,long,bx,by,pd,q,t)

# computes the velocity field on a lat,long grid
# at time t using a mode from pWaves().
#
# input
#
# H    - layer thickness (km)
# Np   - dimensionless buoyancy frequency (N/Omega)
# Br   - radial magnetic field  (mT)
# m    - angular order
# y    - meridional coordinate
# lat  - grid in latitude   (radian)
# long - grid in longitude  (radian)
# bx,by - x and y components of magnetic perturbation
# pd    - period  (years)
# q     - quality factor
# t     - time (in years)
#
# output
#
# vy,vx     -  theta and phi components of velocity at CMB
# ay,ax    -  theta and phi components of acceleration

# reconstruct frequency of mode (in seconds)
sectoyear = 365.25 * 24 * 60 * 60;
tsec = t * sectoyear;
omega0 = 2 * pi / (pd*sectoyear);
omega = omega0 - abs(0.5*omega0/q)im;

# set model parameters (compare with ModelParameters.jl)
B = Br * 1.0e-3;
mu = 4 * pi * 1.0e-7;      # permeability
sigma = 1.0e6;              # conductivity (S/m)
eta = 1.0 / (mu * sigma);   # diffusivity (m^2/s)
kz = pi / (H * 1000.0);      # radial wavenumber
chi = 1.0 + (eta * kz^2/omega)im;

# add southern hemisphere
nlat = length(lat);
nsol = length(y);
yw = [-reverse(y);y];
latw = asin.(yw);
bxw = [zeros(nsol);bx];
byw = [zeros(nsol);by];

#  interpolate wave on to colat,long grid
nodes = (latw,);
itp = Interpolations.interpolate(nodes,bxw,Gridded(Linear()));
bxgrid = itp(lat);
itp= Interpolations.interpolate(nodes,byw,Gridded(Linear()));
bygrid = itp(lat);

# evaluate velocity and acceleration on grid
vx = Array{Float64,2}(undef,length(lat),length(long));
vy = Array{Float64,2}(undef,length(lat),length(long));
ax = Array{Float64,2}(undef,length(lat),length(long));
ay = Array{Float64,2}(undef,length(lat),length(long));

# store in reverse order in latitude
for j = 1 : length(lat)
    for k = 1 : length(long)

       phase = exp((m*long[k] - omega*tsec)im);
       vx[nlat-j+1,k] = real(((omega*chi*bxgrid[j]*phase)/(kz*B))im);
       vy[nlat-j+1,k] = real(((omega*chi*bygrid[j]*phase)/(kz*B))im);
       ax[nlat-j+1,k] = real((omega^2*chi*bxgrid[j]*phase)/(kz*B));
       ay[nlat-j+1,k] = real((omega^2*chi*bygrid[j]*phase)/(kz*B));

    end
end

# normalize and convert to km/yr and km/yr^2
sectoyear = 365.25 * 24 * 60 * 60;
vymax = maximum(vy) * sectoyear;
aymax = maximum(ay) * sectoyear^2;

vx *= sectoyear/vymax;
vy *= sectoyear/vymax;
ax *= sectoyear^2/vymax;
ay *= sectoyear^2/vymax;

return vx,vy,ax,ay

end

ext: jl

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby


text:
const PSA = Pair{Symbol,Any}

compound_declarations = [
"package" => CommandDeclaration[
PSA[:name => "test",
    :api => API.test,
    :should_splat => false,
    :arg_count => 0 => Inf,
    :arg_parser => parse_package,
    :option_spec => [
        PSA[:name => "coverage", :api => :coverage => true],
    ],
    :completions => complete_installed_packages,
    :description => "run tests for packages",
    :help => md"""
    test [--coverage] pkg[=uuid] ...

Run the tests for package `pkg`. This is done by running the file `test/runtests.jl`
in the package directory. The option `--coverage` can be used to run the tests with
coverage enabled. The `startup.jl` file is disabled during testing unless
julia is started with `--startup-file=yes`.
""",
],
PSA[:name => "help",
    :short_name => "?",
    :api => identity, # dummy API function
    :arg_count => 0 => Inf,
    :arg_parser => ((x,y) -> x),
    :completions => complete_help,
    :description => "show this message",
    :help => md"""
    [?|help]

List available commands along with short descriptions.

    [?|help] cmd

If `cmd` is a partial command, display help for all subcommands.
If `cmd` is a full command, display help for `cmd`.
""",
],
PSA[:name => "instantiate",
    :api => API.instantiate,
    :option_spec => [
        PSA[:name => "project", :short_name => "p", :api => :manifest => false],
        PSA[:name => "manifest", :short_name => "m", :api => :manifest => true],
        PSA[:name => "verbose", :short_name => "v", :api => :verbose => true],
    ],
    :description => "downloads all the dependencies for the project",
    :help => md"""
    instantiate [-v|--verbose]
    instantiate [-v|--verbose] [-m|--manifest]
    instantiate [-v|--verbose] [-p|--project]

Download all the dependencies for the current project at the version given by the project's manifest.
If no manifest exists or the `--project` option is given, resolve and download the dependencies compatible with the project.
""",
],
PSA[:name => "remove",
    :short_name => "rm",
    :api => API.rm,
    :should_splat => false,
    :arg_count => 1 => Inf,
    :arg_parser => parse_package,
    :option_spec => [
        PSA[:name => "project", :short_name => "p", :api => :mode => PKGMODE_PROJECT],
        PSA[:name => "manifest", :short_name => "m", :api => :mode => PKGMODE_MANIFEST],
    ],
    :completions => complete_installed_packages,
    :description => "remove packages from project or manifest",
    :help => md"""
    [rm|remove] [-p|--project] pkg[=uuid] ...

Remove package `pkg` from the project file. Since the name `pkg` can only
refer to one package in a project this is unambiguous, but you can specify
a `uuid` anyway, and the command is ignored, with a warning, if package name
and UUID do not match. When a package is removed from the project file, it
may still remain in the manifest if it is required by some other package in
the project. Project mode operation is the default, so passing `-p` or
`--project` is optional unless it is preceded by the `-m` or `--manifest`
options at some earlier point.

    [rm|remove] [-m|--manifest] pkg[=uuid] ...

Remove package `pkg` from the manifest file. If the name `pkg` refers to
multiple packages in the manifest, `uuid` disambiguates it. Removing a package
from the manifest forces the removal of all packages that depend on it, as well
as any no-longer-necessary manifest packages due to project package removals.
""",
],
PSA[:name => "add",
    :api => API.add,
    :should_splat => false,
    :arg_count => 1 => Inf,
    :arg_parser => ((x,y) -> parse_package(x,y; add_or_dev=true)),
    :option_spec => [
        PSA[:name => "preserve", :takes_arg => true, :api => :preserve => do_preserve],
    ],
    :completions => complete_add_dev,
    :description => "add packages to project",
    :help => md"""
    add [--preserve=<opt>] pkg[=uuid] [@version] [#rev] ...

Add package `pkg` to the current project file. If `pkg` could refer to
multiple different packages, specifying `uuid` allows you to disambiguate.
`@version` optionally allows specifying which versions of packages to add. Version specifications
are of the form `@1`, `@1.2` or `@1.2.3`, allowing any version with a prefix
that matches, or ranges thereof, such as `@1.2-3.4.5`. A git revision can be
specified by `#branch` or `#commit`.

If a local path is used as an argument to `add`, the path needs to be a git repository.
The project will then track that git repository just like it would track a remote repository online.
If the package is not located at the top of the git repository, a subdirectory can be specified with
`path:subdir/path`.

`Pkg` resolves the set of packages in your environment using a tiered approach.
The `--preserve` command line option allows you to key into a specific tier in the resolve algorithm.
The following table describes the command line arguments to `--preserve` (in order of strictness).

| Argument | Description                                                                         |
|:---------|:------------------------------------------------------------------------------------|
| `all`    | Preserve the state of all existing dependencies (including recursive dependencies)  |
| `direct` | Preserve the state of all existing direct dependencies                              |
| `semver` | Preserve semver-compatible versions of direct dependencies                          |
| `none`   | Do not attempt to preserve any version information                                  |
| `tiered` | Use the tier which will preserve the most version information (this is the default) |

!!! compat "Julia 1.5"
    Subdirectory specification requires at least Julia 1.5.

**Examples**
```
pkg> add Example
pkg> add --preserve=all Example
pkg> add Example@0.5
pkg> add Example#master
pkg> add Example#c37b675
pkg> add https://github.com/JuliaLang/Example.jl#master
pkg> add git@github.com:JuliaLang/Example.jl.git
pkg> add Example=7876af07-990d-54b4-ab0e-23690620f79a
```
""",
],
PSA[:name => "develop",
    :short_name => "dev",
    :api => API.develop,
    :should_splat => false,
    :arg_count => 1 => Inf,
    :arg_parser => ((x,y) -> parse_package(x,y; add_or_dev=true)),
    :option_spec => [
        PSA[:name => "strict", :api => :strict => true],
        PSA[:name => "local", :api => :shared => false],
        PSA[:name => "shared", :api => :shared => true],
    ],
    :completions => complete_add_dev,
    :description => "clone the full package repo locally for development",
    :help => md"""
    [dev|develop] [--shared|--local] pkg[=uuid] ...

Make a package available for development. If `pkg` is an existing local path, that path will be recorded in
the manifest and used. Otherwise, a full git clone of `pkg` is made. The location of the clone is
controlled by the `--shared` (default) and `--local` arguments. The `--shared` location defaults to
`~/.julia/dev`, but can be controlled with the `JULIA_PKG_DEVDIR` environment variable. When `--local` is given,
the clone is placed in a `dev` folder in the current project.
This operation is undone by `free`.

**Examples**
```jl
pkg> develop Example
pkg> develop https://github.com/JuliaLang/Example.jl
pkg> develop ~/mypackages/Example
pkg> develop --local Example
```
""",
],
PSA[:name => "free",
    :api => API.free,
    :should_splat => false,
    :arg_count => 1 => Inf,
    :arg_parser => parse_package,
    :completions => complete_installed_packages,
    :description => "undoes a `pin`, `develop`, or stops tracking a repo",
    :help => md"""
    free pkg[=uuid] ...

Free a pinned package `pkg`, which allows it to be upgraded or downgraded again. If the package is checked out (see `help develop`) then this command
makes the package no longer being checked out.
""",
],
PSA[:name => "pin",
    :api => API.pin,
    :should_splat => false,
    :arg_count => 1 => Inf,
    :arg_parser => parse_package,
    :completions => complete_installed_packages,
    :description => "pins the version of packages",
    :help => md"""
    pin pkg[=uuid] ...

Pin packages to given versions, or the current version if no version is specified. A pinned package has its version fixed and will not be upgraded or downgraded.
A pinned package has the symbol `⚲` next to its version in the status list.

**Examples**
```
pkg> pin Example
pkg> pin Example@0.5.0
pkg> pin Example=7876af07-990d-54b4-ab0e-23690620f79a@0.5.0
```
""",
],
PSA[:name => "build",
    :api => API.build,
    :should_splat => false,
    :arg_count => 0 => Inf,
    :arg_parser => parse_package,
    :option_spec => [
        PSA[:name => "verbose", :short_name => "v", :api => :verbose => true],
    ],
    :completions => complete_installed_packages,
    :description => "run the build script for packages",
    :help => md"""
    build [-v|--verbose] pkg[=uuid] ...

Run the build script in `deps/build.jl` for `pkg` and all of its dependencies in depth-first recursive order.
If no packages are given, run the build scripts for all packages in the manifest.
The `-v`/`--verbose` option redirects build output to `stdout`/`stderr` instead of the `build.log` file.
The `startup.jl` file is disabled during building unless julia is started with `--startup-file=yes`.
""",
],
PSA[:name => "resolve",
    :api => API.resolve,
    :description => "resolves to update the manifest from changes in dependencies of developed packages",
    :help => md"""
    resolve

Resolve the project i.e. run package resolution and update the Manifest. This is useful in case the dependencies of developed
packages have changed causing the current Manifest to be out of sync.
""",
],
PSA[:name => "activate",
    :api => API.activate,
    :arg_count => 0 => 1,
    :arg_parser => parse_activate,
    :option_spec => [
        PSA[:name => "shared", :api => :shared => true],
        PSA[:name => "temp", :api => :temp => true],
    ],
    :completions => complete_activate,
    :description => "set the primary environment the package manager manipulates",
    :help => md"""
    activate
    activate [--shared|--temp] [path]

Activate the environment at the given `path`, or the home project environment if no `path` is specified.
The active environment is the environment that is modified by executing package commands.
When the option `--shared` is given, `path` will be assumed to be a directory name and searched for in the
`environments` folders of the depots in the depot stack. In case no such environment exists in any of the depots,
it will be placed in the first depot of the stack.
Use the `--temp` option to create temporary environments which are removed when the julia
process is exited.
""" ,
],
PSA[:name => "update",
    :short_name => "up",
    :api => API.up,
    :should_splat => false,
    :arg_count => 0 => Inf,
    :arg_parser => parse_package,
    :option_spec => [
        PSA[:name => "project",  :short_name => "p", :api => :mode => PKGMODE_PROJECT],
        PSA[:name => "manifest", :short_name => "m", :api => :mode => PKGMODE_MANIFEST],
        PSA[:name => "major", :api => :level => UPLEVEL_MAJOR],
        PSA[:name => "minor", :api => :level => UPLEVEL_MINOR],
        PSA[:name => "patch", :api => :level => UPLEVEL_PATCH],
        PSA[:name => "fixed", :api => :level => UPLEVEL_FIXED],
    ],
    :completions => complete_installed_packages,
    :description => "update packages in manifest",
    :help => md"""
    [up|update] [-p|--project]  [opts] pkg[=uuid] [@version] ...
    [up|update] [-m|--manifest] [opts] pkg[=uuid] [@version] ...

    opts: --major | --minor | --patch | --fixed

Update `pkg` within the constraints of the indicated version
specifications. These specifications are of the form `@1`, `@1.2` or `@1.2.3`, allowing
any version with a prefix that matches, or ranges thereof, such as `@1.2-3.4.5`.
In `--project` mode, package specifications only match project packages, while
in `manifest` mode they match any manifest package. Bound level options force
the following packages to be upgraded only within the current major, minor,
patch version; if the `--fixed` upgrade level is given, then the following
packages will not be upgraded at all.
""",
],
PSA[:name => "generate",
    :api => API.generate_deprecated,
    :arg_count => 1 => 1,
    :arg_parser => ((x,y) -> map(expanduser, unwrap(x))),
    :description => "generate files for a new project",
    :help => md"""
    generate pkgname

Create a project called `pkgname` in the current folder.
""",
],
PSA[:name => "precompile",
    :api => API.precompile,
    :description => "precompile all the project dependencies",
    :help => md"""
    precompile

Precompile all the dependencies of the project in parallel.
The `startup.jl` file is disabled during precompilation unless julia is started with `--startup-file=yes`.

Errors will only throw when precompiling the top-level dependencies, given that
not all manifest dependencies may be loaded by the top-level dependencies on the given system.

This method is called automatically after any Pkg action that changes the manifest.
Any packages that have previously errored during precompilation won't be retried in auto mode
until they have changed. To disable automatic precompilation set the environment variable `JULIA_PKG_PRECOMPILE_AUTO=0`
""",
],
PSA[:name => "status",
    :short_name => "st",
    :api => API.status,
    :should_splat => false,
    :arg_count => 0 => Inf,
    :arg_parser => parse_package,
    :option_spec => [
        PSA[:name => "project",  :short_name => "p", :api => :mode => PKGMODE_PROJECT],
        PSA[:name => "manifest", :short_name => "m", :api => :mode => PKGMODE_MANIFEST],
        PSA[:name => "diff", :short_name => "d", :api => :diff => true],
    ],
    :completions => complete_installed_packages,
    :description => "summarize contents of and changes to environment",
    :help => md"""
    [st|status] [-d|--diff] [pkgs...]
    [st|status] [-d|--diff] [-p|--project] [pkgs...]
    [st|status] [-d|--diff] [-m|--manifest] [pkgs...]

Show the status of the current environment. In `--project` mode (default), the
status of the project file is summarized. In `--manifest` mode the output also
includes the recursive dependencies of added packages given in the manifest.
If there are any packages listed as arguments the output will be limited to those packages.
The `--diff` option will, if the environment is in a git repository, limit
the output to the difference as compared to the last git commit.

!!! compat "Julia 1.1"
    `pkg> status` with package arguments requires at least Julia 1.1.

!!! compat "Julia 1.3"
    The `--diff` option requires Julia 1.3. In earlier versions `--diff`
    is the default for environments in git repositories.
""",
],
PSA[:name => "gc",
    :api => API.gc,
    :option_spec => [
        PSA[:name => "all", :api => :collect_delay => Hour(0)],
        PSA[:name => "verbose", :short_name => "v", :api => :verbose => true],
    ],
    :description => "garbage collect packages not used for a significant time",
    :help => md"""
    gc [-v|--verbose] [--all]

Free disk space by garbage collecting packages not used for a significant time.
The `--all` option will garbage collect all packages which can not be immediately
reached from any existing project.
Use verbose mode for detailed output.
""",
],
PSA[:name => "undo",
    :api => API.undo,
    :description => "undo the latest change to the active project",
    :help => md"""
    undo

Undoes the latest change to the active project.
""",
],
PSA[:name => "redo",
    :api => API.redo,
    :description => "redo the latest change to the active project",
    :help => md"""
    redo

Redoes the changes from the latest [`undo`](@ref).
""",
],
], #package
"registry" => CommandDeclaration[
PSA[:name => "add",
    :api => Registry.add,
    :should_splat => false,
    :arg_count => 0 => Inf,
    :arg_parser => ((x,y) -> parse_registry(x,y; add = true)),
    :description => "add package registries",
    :help => md"""
    registry add reg...

Add package registries `reg...` to the user depot. Without arguments
it adds known registries, i.e. the General registry and registries
served by the configured package server.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

!!! compat "Julia 1.5"
    `registry add` without arguments requires at least Julia 1.5.

**Examples**
```
pkg> registry add General
pkg> registry add https://www.my-custom-registry.com
pkg> registry add
```
""",
],
PSA[:name => "remove",
    :short_name => "rm",
    :api => Registry.rm,
    :should_splat => false,
    :arg_count => 1 => Inf,
    :arg_parser => parse_registry,
    :description => "remove package registries",
    :help => md"""
    registry [rm|remove] reg...

Remove package registries `reg...`.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

**Examples**
```
pkg> registry [rm|remove] General
```
""",
],
PSA[:name => "update",
    :short_name => "up",
    :api => Registry.update,
    :should_splat => false,
    :arg_count => 0 => Inf,
    :arg_parser => parse_registry,
    :description => "update package registries",
    :help => md"""
    registry [up|update]
    registry [up|update] reg...

Update package registries `reg...`. If no registries are specified
all registries will be updated.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

**Examples**
```
pkg> registry up
pkg> registry up General
```
""",
],
PSA[:name => "status",
    :short_name => "st",
    :api => Registry.status,
    :description => "information about installed registries",
    :help => md"""
    registry [st|status]

Display information about installed registries.

!!! compat "Julia 1.1"
    Pkg's registry handling requires at least Julia 1.1.

**Examples**
```
pkg> registry status
```
""",
]
], #registry
] #command_declarations


text:
abstract type Objective end

mutable struct TrackingObjective{Q,U,C,B} <: Objective
    q::Vector{Q}
    u::Vector{U}
    γ::Vector{C}
    b::Vector{B}
end

function TrackingObjective(model, env, H::Int;
    q = [Diagonal(zeros(SizedVector{model.dim.q})) for t = 1:H],
    u = [Diagonal(zeros(SizedVector{model.dim.u})) for t = 1:H],
    γ = [Diagonal(zeros(SizedVector{model.dim.c})) for t = 1:H],
    b = [Diagonal(zeros(SizedVector{model.dim.c * friction_dim(env)})) for t = 1:H])
    return TrackingObjective(q, u, γ, b)
end

mutable struct TrackingVelocityObjective{Q,V,U,C,B} <: Objective
    q::Vector{Q}
    v::Vector{V}
    u::Vector{U}
    γ::Vector{C}
    b::Vector{B}
end

function TrackingVelocityObjective(model, env, H::Int;
    q = [Diagonal(zeros(SizedVector{model.dim.q})) for t = 1:H],
    v = [Diagonal(zeros(SizedVector{model.dim.q})) for t = 1:H],
    u = [Diagonal(zeros(SizedVector{model.dim.u})) for t = 1:H],
    γ = [Diagonal(zeros(SizedVector{model.dim.c})) for t = 1:H],
    b = [Diagonal(zeros(SizedVector{model.dim.c * friction_dim(env)})) for t = 1:H])
    return TrackingVelocityObjective(q, v, u, γ, b)
end


text:
# # Horizontal convection example
#
# In "horizontal convection", a non-uniform buoyancy is imposed on top of an initially resting fluid.
#
# This example demonstrates:
#
#   * How to use computed `Field`s for output.
#   * How to post-process saved output using `FieldTimeSeries`.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, CairoMakie"
# ```

# ## Horizontal convection
#
# We consider here two-dimensional horizontal convection of an incompressible flow ``\boldsymbol{u} = (u, w)``
# on the ``(x, z)``-plane (``-L_x/2 \le x \le L_x/2`` and ``-H \le z \le 0``). The flow evolves
# under the effect of gravity. The only forcing on the fluid comes from a prescribed, non-uniform
# buoyancy at the top-surface of the domain.
#
# We start by importing `Oceananigans` and `Printf`.

using Oceananigans
using Printf

# ### The grid

H = 1.0          # vertical domain extent
Lx = 2H          # horizontal domain extent
Nx, Nz = 128, 64 # horizontal, vertical resolution

grid = RectilinearGrid(size = (Nx, Nz),
                          x = (-Lx/2, Lx/2),
                          z = (-H, 0),
                       halo = (3, 3),
                   topology = (Bounded, Flat, Bounded))

# ### Boundary conditions
#
# At the surface, the imposed buoyancy is
# ```math
# b(x, z = 0, t) = - b_* \cos (2 \pi x / L_x) \, ,
# ```
# while zero-flux boundary conditions are imposed on all other boundaries. We use free-slip
# boundary conditions on ``u`` and ``w`` everywhere.

b★ = 1.0

@inline bˢ(x, y, t, p) = - p.b★ * cos(2π * x / p.Lx)

b_bcs = FieldBoundaryConditions(top = ValueBoundaryCondition(bˢ, parameters=(; b★, Lx)))

# ### Non-dimensional control parameters and Turbulence closure
#
# The problem is characterized by three non-dimensional parameters. The first is the domain's
# aspect ratio, ``L_x / H`` and the other two are the Rayleigh (``Ra``) and Prandtl (``Pr``)
# numbers:
#
# ```math
# Ra = \frac{b_* L_x^3}{\nu \kappa} \, , \quad \text{and}\, \quad Pr = \frac{\nu}{\kappa} \, .
# ```
#
# The Prandtl number expresses the ratio of momentum over heat diffusion; the Rayleigh number
# is a measure of the relative importance of gravity over viscosity in the momentum equation.
#
# For a domain with a given extent, the nondimensional values of ``Ra`` and ``Pr`` uniquely
# determine the viscosity and diffusivity, i.e.,
#
# ```math
# \nu = \sqrt{\frac{Pr b_* L_x^3}{Ra}} \quad \text{and} \quad \kappa = \sqrt{\frac{b_* L_x^3}{Pr Ra}} \, .
# ```
#
# We use isotropic viscosity and diffusivities, `ν` and `κ` whose values are obtain from the
# prescribed ``Ra`` and ``Pr`` numbers. Here, we use ``Pr = 1`` and ``Ra = 10^8``:

Pr = 1.0    # Prandtl number
Ra = 1e8    # Rayleigh number

ν = sqrt(Pr * b★ * Lx^3 / Ra)  # Laplacian viscosity
κ = ν * Pr                     # Laplacian diffusivity
nothing # hide

# ## Model instantiation
#
# We instantiate the model with the fifth-order WENO advection scheme, a 3rd order
# Runge-Kutta time-stepping scheme, and a `BuoyancyTracer`.

model = NonhydrostaticModel(; grid,
                            advection = WENO5(),
                            timestepper = :RungeKutta3,
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            closure = ScalarDiffusivity(; ν, κ),
                            boundary_conditions = (; b=b_bcs))

# ## Simulation set-up
#
# We set up a simulation that runs up to ``t = 40`` with a `JLD2OutputWriter` that saves the flow
# speed, ``\sqrt{u^2 + w^2}``, the buoyancy, ``b``, andthe vorticity, ``\partial_z u - \partial_x w``.

simulation = Simulation(model, Δt=1e-2, stop_time=40.0)

# ### The `TimeStepWizard`
#
# The TimeStepWizard manages the time-step adaptively, keeping the Courant-Freidrichs-Lewy
# (CFL) number close to `0.75` while ensuring the time-step does not increase beyond the
# maximum allowable value for numerical stability.

wizard = TimeStepWizard(cfl=0.75, max_change=1.2, max_Δt=1e-1)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(50))

# ### A progress messenger
#
# We write a function that prints out a helpful progress message while the simulation runs.

progress(sim) = @printf("i: % 6d, sim time: % 1.3f, wall time: % 10s, Δt: % 1.4f, CFL: %.2e\n",
                        iteration(sim), time(sim), prettytime(sim.run_wall_time),
                        sim.Δt, AdvectiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(50))

# ### Output
#
# We use computed `Field`s to diagnose and output the total flow speed, the vorticity, ``\zeta``,
# and the buoyancy, ``b``. Note that computed `Field`s take "AbstractOperations" on `Field`s as
# input:

u, v, w = model.velocities # unpack velocity `Field`s
b = model.tracers.b        # unpack buoyancy `Field`

## total flow speed
s = sqrt(u^2 + w^2)

## y-component of vorticity
ζ = ∂z(u) - ∂x(w)
nothing # hide

# We create a `JLD2OutputWriter` that saves the speed, and the vorticity.
# We then add the `JLD2OutputWriter` to the `simulation`.

saved_output_filename = "horizontal_convection.jld2"

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; s, b, ζ),
                                                      schedule = TimeInterval(0.5),
                                                      filename = saved_output_filename,
                                                      overwrite_existing = true)
nothing # hide

# Ready to press the big red button:

run!(simulation)

# ## Load saved output, process, visualize
#
# We animate the results by opening the JLD2 file, extracting data for the iterations we ended
# up saving at, and ploting the saved fields. From the saved buoyancy field we compute the
# buoyancy dissipation, ``\chi = \kappa |\boldsymbol{\nabla} b|^2``, and plot that also.
#
# To start we load the saved fields are `FieldTimeSeries` and prepare for animating the flow by
# creating coordinate arrays that each field lives on.

using CairoMakie
using Oceananigans
using Oceananigans.Fields
using Oceananigans.AbstractOperations: volume

saved_output_filename = "horizontal_convection.jld2"

## Open the file with our data
s_timeseries = FieldTimeSeries(saved_output_filename, "s")
b_timeseries = FieldTimeSeries(saved_output_filename, "b")
ζ_timeseries = FieldTimeSeries(saved_output_filename, "ζ")

times = b_timeseries.times

## Coordinate arrays
xs, ys, zs = nodes(s_timeseries[1])
xb, yb, zb = nodes(b_timeseries[1])
xζ, yζ, zζ = nodes(ζ_timeseries[1])
nothing # hide

χ_timeseries = deepcopy(b_timeseries)

for i in 1:length(times)
  bᵢ = b_timeseries[i]
  χ_timeseries[i] .= κ * (∂x(bᵢ)^2 + ∂z(bᵢ)^2)
end

# Now we're ready to animate using Makie.

@info "Making an animation from saved data..."

n = Observable(1)

title = @lift @sprintf("t=%1.2f", times[$n])

sₙ = @lift interior(s_timeseries[$n], :, 1, :)
ζₙ = @lift interior(ζ_timeseries[$n], :, 1, :)
bₙ = @lift interior(b_timeseries[$n], :, 1, :)
χₙ = @lift interior(χ_timeseries[$n], :, 1, :)

slim = 0.6
blim = 0.6
ζlim = 9
χlim = 0.025

axis_kwargs = (xlabel = "x / H",
               ylabel = "z / H",
               limits = ((-Lx/2, Lx/2), (-H, 0)),
               aspect = Lx/H,
               titlesize = 20)

fig = Figure(resolution = (600, 1100))

ax_s = Axis(fig[2, 1];
            title = "speed √[(u²+w²)/(b⋆H)]", axis_kwargs...)

ax_b = Axis(fig[3, 1];
            title = "buoyancy, b/b⋆", axis_kwargs...)

ax_ζ = Axis(fig[4, 1];
            title = "vorticity, (∂u/∂z - ∂w/∂x) √(H/b⋆)", axis_kwargs...)

ax_χ = Axis(fig[5, 1];
            title = "buoyancy dissipation, κ|∇b|² √(H/b⋆⁵)", axis_kwargs...)

fig[1, :] = Label(fig, title, textsize=24, tellwidth=false)

hm_s = heatmap!(ax_s, xs, zs, sₙ;
                colorrange = (0, slim),
                colormap = :speed)
Colorbar(fig[2, 2], hm_s)

hm_b = heatmap!(ax_b, xb, zb, bₙ;
                colorrange = (-blim, blim),
                colormap = :thermal)
Colorbar(fig[3, 2], hm_b)

hm_ζ = heatmap!(ax_ζ, xζ, zζ, ζₙ;
                colorrange = (-ζlim, ζlim),
                colormap = :balance)
Colorbar(fig[4, 2], hm_ζ)

hm_χ = heatmap!(ax_χ, xs, zs, χₙ;
                colorrange = (0, χlim),
                colormap = :dense)
Colorbar(fig[5, 2], hm_χ)

# And, finally, we record a movie.

frames = 1:length(times)

record(fig, "horizontal_convection.mp4", frames, framerate=8) do i
    msg = string("Plotting frame ", i, " of ", frames[end])
    print(msg * " \r")
    n[] = i
end
nothing #hide

# ![](horizontal_convection.mp4)

# At higher Rayleigh numbers the flow becomes much more vigorous. See, for example, an animation
# of the voricity of the fluid at ``Ra = 10^{12}`` on [vimeo](https://vimeo.com/573730711).

# ### The Nusselt number
#
# Often we are interested on how much the flow enhances mixing. This is quantified by the
# Nusselt number, which measures how much the flow enhances mixing compared if only diffusion
# was in operation. The Nusselt number is given by
#
# ```math
# Nu = \frac{\langle \chi \rangle}{\langle \chi_{\rm diff} \rangle} \, ,
# ```
#
# where angle brackets above denote both a volume and time average and  ``\chi_{\rm diff}`` is
# the buoyancy dissipation that we get without any flow, i.e., the buoyancy dissipation associated
# with the buoyancy distribution that satisfies
#
# ```math
# \kappa \nabla^2 b_{\rm diff} = 0 \, ,
# ```
#
# with the same boundary conditions same as our setup. In this case we can readily find that
#
# ```math
# b_{\rm diff}(x, z) = b_s(x) \frac{\cosh \left [2 \pi (H + z) / L_x \right ]}{\cosh(2 \pi H / L_x)} \, ,
# ```
# which implies ``\langle \chi_{\rm diff} \rangle = \kappa b_*^2 \pi \tanh(2 \pi Η /Lx)``.
#
# We use the loaded `FieldTimeSeries` to compute the Nusselt number from buoyancy and the volume
# average kinetic energy of the fluid.
#
# First we compute the diffusive buoyancy dissipation, ``\chi_{\rm diff}`` (which is just a
# scalar):

χ_diff = κ * b★^2 * π * tanh(2π * H / Lx)
nothing # hide

# We then create two `ReducedField`s to store the volume integrals of the kinetic energy density
# and the buoyancy dissipation. We need the `grid` to do so; the `grid` can be recoverd from
# any `FieldTimeSeries`, e.g.,

grid = b_timeseries.grid

∫ⱽ_s² = Field{Nothing, Nothing, Nothing}(grid)
∫ⱽ_mod²_∇b = Field{Nothing, Nothing, Nothing}(grid)

# We recover the time from the saved `FieldTimeSeries` and construct two empty arrays to store
# the volume-averaged kinetic energy and the instantaneous Nusselt number,

t = b_timeseries.times

kinetic_energy, Nu = zeros(length(t)), zeros(length(t))
nothing # hide

# Now we can loop over the fields in the `FieldTimeSeries`, compute `KineticEnergy` and `Nu`,
# and plot.

for i = 1:length(t)
    s = s_timeseries[i]
    sum!(∫ⱽ_s², s^2 * volume)
    kinetic_energy[i] = 0.5 * ∫ⱽ_s²[1, 1, 1]  / (Lx * H)

    b = b_timeseries[i]
    sum!(∫ⱽ_mod²_∇b, (∂x(b)^2 + ∂z(b)^2) * volume)
    Nu[i] = (κ *  ∫ⱽ_mod²_∇b[1, 1, 1]) / χ_diff
end

fig = Figure(resolution = (850, 450))

ax_KE = Axis(fig[1, 1], xlabel = "time", ylabel = "KE / (b⋆H)")
lines!(ax_KE, t, kinetic_energy; linewidth = 3)

ax_Nu = Axis(fig[2, 1], xlabel = "time", ylabel = "Nu")
lines!(ax_Nu, t, Nu; linewidth = 3)

current_figure() # hide


text:
#ENV["JULIA_NUM_THREADS"] = "4"
#ENV["UV_THREADPOOL_SIZE"] = "4"
#ENV["JULIA_REVISE"] = "manual" # manually call revise() to update code
#ENV["JULIA_REVISE_INCLUDE"] = "1" # automatic tracking included files
#ENV["JULIA_REVISE_POLL"] = "1"
ENV["AI_THREADS"] = "1"
julia = joinpath(Sys.BINDIR, "julia")
start = abspath(@__DIR__,"startup.jl")
run(`$julia $start`)

max_stars_repo_name: Gilga/AI-Cywwtaip

max_issues_repo_name: Gilga/AI-Cywwtaip

max_forks_repo_name: Gilga/AI-Cywwtaip


text:
@testset "Map" begin

    using Tables, DataFrames
    df = Tables.namedtupleiterator(DataFrame(Chrom = repeat(["chr1"], 3), Pos = [1,2,3]))

    m = df |> GOR.map(x -> Base.merge(x, (Value = x.Pos * 2,)))

    @test length(m) == 3
    @test m |> collect |> length == 3
    @test m |> DataFrame |> size == (3, 3)
    @test eltype(m) == NamedTuple{(:Chrom, :Pos, :Value),Tuple{String,Int64,Int64}}

end


text:
## Sparse matrix performance
include("../perfutil.jl")

## create some sparse matrices (need to be square):
seed = 1
srand(seed)

small = 10^3
med = 10^4
large = 10^5
huge = 10^6
# # 1 entry per line
# ss = {}
# push!(ss, sprand(small, small, 1e-3))
# #push!(ss, sprand(med, med, 1e-4))
# push!(ss, sprand(large, large, 1e-5))

# 10 entries per line
ts = {}
push!(ts, sprand(small, small, 1e-2))
#push!(ts, sprand(med, med, 1e-3))
push!(ts, sprand(large, large, 1e-4))

# 100 entries per line
us = {}
push!(us, sprand(small, small, 1e-1))
#push!(us, sprand(med, med, 1e-2))
push!(us, sprand(large, large, 1e-3))
#push!(us, sprand(huge, huge, 1e-4))

# # 1000 entries per line
# vs = {}
# push!(vs, sprand(small, small, 1.0))
# #push!(vs, sprand(med, med, 1e-1))
# push!(vs, sprand(large, large, 1e-2))
# #push!(vs, sprand(huge, huge, 1e-3))

## using uint32 (works up to 10^9)
uus = {}
for u in us
    push!(uus, SparseMatrixCSC(u.m, u.n, uint32(u.colptr), uint32(u.rowval), u.nzval))
end

## getindex
rep = 20
reps = rep^2
function integer_indexing(A)
    # index with two random integers
    nI, nJ = size(A)
    rI = 1:nI
    rJ = 1:nJ
    tmp = zero(eltype(A))
    for i in rand(rI, reps)
        for j in rand(rJ, rep)
            tmp += A[i,j]
        end
    end
    tmp
end

function row_indexing(A, rowinds)
    # index rows with rowinds and columns with a random integer
    nI, nJ = size(A)
    rI = 1:nI
    rJ = 1:nJ
    tmp = zero(eltype(A))
    for j in rand(rJ, reps)
        tmp += sum(A[rowinds,j])
    end
    tmp
end

function col_indexing(A, colinds)
    # index rows with a random integer and columns with colinds
    nI, nJ = size(A)
    rI = 1:nI
    rJ = 1:nJ
    tmp = zero(eltype(A))
    for i in rand(rI, div(reps,10) )
        tmp += sum(A[i,colinds])
    end
    tmp
end

function row_col_indexing(A, rowinds, colinds)
    # index rows with rowinds and columns with colinds
    # we need:
    (maximum(rowinds)+rep < size(A,1) && maximum(colinds)+rep < size(A, 2)) || error("bad rowinds or colinds")
    nI, nJ = size(A)
    rI = 1:nI
    rJ = 1:nJ
    for i in 1:10
        for j in 1:10
            tmp2 = A[rowinds.+i, colinds.+j]
        end
    end
end

function one_arg_indexing(A, lininds)
    # This is for 1d-indexing and indexing with one array of logicals.
    # Both return a nx1 sparse matrix.
    tmp = zero(eltype(A))
    if isa(eltype(A), Bool)
        tmp = sum(A[lininds])
    else
        for i in 1:rep
            tmp += sum(A[lininds])
        end
    end
    tmp
end

# test performance with matrices in us, uus and ts for
# - integer indexing
# - range indexing
# - ordered array indexing
# - disordered array indexing
# - 1d indexing with integers and booleans

# setup:
# - indices
intinds = nothing
logicalinds = nothing # needs to be generated for a specific matrix size.
rangeinds = 121:237
orderedinds = [rangeinds]
disorderedinds = orderedinds[randperm(length(orderedinds))]

inds = [(intinds, "integers"), (logicalinds, "logical array"), (rangeinds, "a range"),
        (orderedinds, "a ordered array"), (disorderedinds, "a disordered array")]

# - matrix sizes
sizes = [(1, "small", "Small sparse matrix"), (2, "medium", "Medium sparse matrix")]

# - matrix types
mattyp = [(ts, "10 entries/column"), (us, "100 entries/column"), (uus, "100 entries/column uint32")]

# - functions
funs = [(integer_indexing, 1, "indexing"), (one_arg_indexing, 1, "1d indexing"),
        (row_indexing, 2, "indexing rows"), (col_indexing, 2, "indexing columns"),
        (row_col_indexing, 3, "indexing rows & columns")]

# performance tests:
counters = [1,1] # for small and medium matrix
# integer indexing
for (sz,s1,s2) in sizes # size of the matrix
    for (mt, ms) in mattyp  # type of the matrix
        m = mt[sz]
        c = counters[sz]
        @timeit integer_indexing(m) "sparse_getindex_$s1$c" "$s2 with $ms, indexing with integers"
        counters[sz] += 1
    end
end

# range & array indexing
for (sz,s1,s2) in sizes # size of the matrix
    for (mt, ms) in mattyp  # type of the matrix
        m = mt[sz]
        for (fun, nargs, funstr) in funs[3:5] # indexing test function
            for (ind,indstr) in inds[3:5] # type of indices
                c = counters[sz]
                if indstr=="logical array"
                    # make a logical array of the right size
                    ind = sprandbool(size(m,1)..., 1e-5)
                end
                if nargs==2
                    @timeit fun(m, ind)  "sparse_getindex_$s1$c" "Sparse matrix with $ms, $funstr with $indstr"
                elseif nargs==3
                    @timeit fun(m, ind, ind) "sparse_getindex_$s1$c" "Sparse matrix with $ms, $funstr with $indstr"
                else
                    error("Something is amiss here.")
                end
                counters[sz] += 1
            end
        end
    end
end

# linear indexing
for (sz,s1,s2) in sizes # size of the matrix
    for (mt, ms) in mattyp  # type of the matrix
        m = mt[sz]
        for (ind,indstr) in inds[2:5] # type of indices
            if indstr=="logical array"
                if sz==2
                    continue # logical indexing with medium size sparse matrix takes too long
                end
                # make a logical array of the right size
                ind = sprandbool(size(m)..., 1e-5)
                c = counters[sz]
                @timeit one_arg_indexing(m, ind) "sparse_getindex_$s1$c" "$s2 with $ms, linear indexing with $indstr"
                counters[sz] += 1
            else
                c = counters[sz]
                @timeit one_arg_indexing(m, ind) "sparse_getindex_$s1$c" "$s2 with $ms, linear indexing with $indstr"
                counters[sz] += 1
            end
        end
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#######################################################################
import JSON, Conda
using Compat
import Compat.String

# remove deps.jl at exit if it exists, in case build.jl fails
try
#######################################################################

# print to stderr, since that is where Pkg prints its messages
eprintln(x...) = println(STDERR, x...)

# Make sure Python uses UTF-8 output for Unicode paths
ENV["PYTHONIOENCODING"] = "UTF-8"

function prog_version(prog)
    try
       return convert(VersionNumber, chomp(readstring(`$prog --version`)))
    catch
       return v"0.0"
    end
end

jupyter = ""
jupyter_vers = v"0.0"
for p in (haskey(ENV, "JUPYTER") ? (ENV["JUPYTER"],) : (isfile("JUPYTER") ? readchomp("JUPYTER") : "jupyter", "jupyter", "ipython", "ipython2", "ipython3", "ipython.bat"))
    v = prog_version(p)
    if v >= v"3.0"
       jupyter = p
       jupyter_vers = v
       break
    end
end
if jupyter_vers < v"3.0" || dirname(jupyter) == abspath(Conda.SCRIPTDIR)
    info("Installing Jupyter via the Conda package.")
    Conda.add("jupyter")
    jupyter = abspath(Conda.SCRIPTDIR,"jupyter")
    jupyter_vers = prog_version(jupyter)
    jupyter_vers < v"3.0" && error("failed to find $jupyter 3.0 or later")
end
info("Found Jupyter version $jupyter_vers: $jupyter")

#######################################################################
# Warn people upgrading from older IJulia versions:
try
    juliaprof = chomp(readstring(pipeline(`$ipython locate profile julia`,
                                          stderr=DevNull)))
    warn("""You should now run IJulia just via `$jupyter notebook`, without
            the `--profile julia` flag.  IJulia no longer maintains the profile.
            Consider deleting $juliaprof""")
end

#######################################################################
# Install Jupyter kernel-spec file.

# Is IJulia being built from a debug build? If so, add "debug" to the description.
debugdesc = ccall(:jl_is_debugbuild,Cint,())==1 ? "-debug" : ""

spec_name = "julia-$(VERSION.major).$(VERSION.minor)"*debugdesc
juliakspec = abspath(spec_name)

binary_name = is_windows() ? "julia.exe" : "julia"
kernelcmd_array = String[joinpath(JULIA_HOME,("$binary_name")), "-i"]
ijulia_dir = get(ENV, "IJULIA_DIR", Pkg.dir("IJulia")) # support non-Pkg IJulia installs
append!(kernelcmd_array, ["--startup-file=yes", "--color=yes", joinpath(ijulia_dir,"src","kernel.jl"), "{connection_file}"])

ks = Dict(
    "argv" => kernelcmd_array,
    "display_name" => "Julia " * Base.VERSION_STRING * debugdesc,
    "language" => "julia",
)

destname = "kernel.json"
mkpath(juliakspec)
dest = joinpath(juliakspec, destname)

eprintln("Writing IJulia kernelspec to $dest ...")

open(dest, "w") do f
    # indent by 2 for readability of file
    write(f, JSON.json(ks, 2))
end

copy_config(src, dest) = cp(src, joinpath(dest, src), remove_destination=true)

copy_config("logo-32x32.png", juliakspec)
copy_config("logo-64x64.png", juliakspec)

eprintln("Installing julia kernelspec $spec_name")

# remove these hacks when
# https://github.com/jupyter/notebook/issues/448 is closed and the fix
# is widely available -- just run `$jupyter kernelspec ...` then.
notebook = String[]
try
    run(`$jupyter kernelspec install --replace --user $juliakspec`)
    push!(notebook, jupyter, "notebook")
catch
    @static if is_unix()
        run(`$jupyter-kernelspec install --replace --user $juliakspec`)
        push!(notebook, jupyter * "-notebook")
    end

    # issue #363:
    @static if is_windows()
        jupyter_dir = dirname(jupyter)
        jks_exe = ""
        if jupyter_dir == abspath(Conda.SCRIPTDIR)
            jk_path = "$jupyter-kernelspec"
            if isfile(jk_path * "-script.py")
                jk_path *= "-script.py"
            end
            jn_path = "$jupyter-notebook"
            if isfile(jn_path * "-script.py")
                jn_path *= "-script.py"
            end
            python = abspath(Conda.PYTHONDIR, "python.exe")
        else
            jks_exe = joinpath(jupyter_dir, "jupyter-kernelspec.exe")
            if !isfile(jks_exe)
                jk_path = readchomp(`where.exe $jupyter-kernelspec`)
                jn_path = readchomp(`where.exe $jupyter-notebook`)
                # jupyter-kernelspec should start with "#!/path/to/python":
                python = strip(chomp(open(readline, jk_path, "r"))[3:end])
                # strip quotes, if any
                if python[1] == python[end] == '"'
                    python = python[2:end-1]
                end
            else
                jn_path = joinpath(jupyter_dir, "jupyter-notebook.exe")
                isfile(jn_path) || error("$jn_path not found")
            end
        end
        if isfile(jks_exe)
            run(`$jks_exe install --replace --user $juliakspec`)
        else
            run(`$python $jk_path install --replace --user $juliakspec`)
        end
        if endswith(jn_path, ".exe")
            push!(notebook, jn_path)
        else
            push!(notebook, python, jn_path)
        end
    end
end
deps = """
    const jupyter = "$(escape_string(jupyter))"
    const notebook_cmd = ["$(join(map(escape_string, notebook), "\", \""))"]
    const jupyter_vers = $(repr(jupyter_vers))
    """
if !isfile("deps.jl") || readstring("deps.jl") != deps
    write("deps.jl", deps)
end
write("JUPYTER", jupyter)

#######################################################################
catch
isfile("deps.jl") && rm("deps.jl") # remove deps.jl file on build error
end


text:
# Tests for the parzen density estimator and classifier
function t_parzen()

Ac = j4pr.DataGenerator.fish(20)   	# classification dataset
Ade = j4pr.datacell(collect(-1:0.1:1)) 	# density estimation dataset

Wclass = [j4pr.parzen(h, window=Φ, metric=m) for h in [0.1,1,3],
	  						  Φ in [:hat, :linear, :gaussian, :exponential, :cosine],
							  m in [Distances.Euclidean(), Distances.Cityblock()]
	]

Wde = [j4pr.parzen(h, window=Φ, metric=m) for h in [0.1,1,3],
	  						  Φ in [:hat, :linear, :gaussian, :exponential, :cosine],
							  m in [Distances.Euclidean(), Distances.Cityblock()]
	]

# Test classification
for w in Wclass
	Test.@test try
		wt1 = Ac |> w
		wt2 = j4pr.strip(Ac) |> w

		result = Ac |> wt1
		result2 = +Ac |> wt2
		true
	catch
		false
	end
end

# Test density estimation
for w in Wde
	Test.@test try
		wt1 = Ade |> w
		wt2 = +Ade |> w

		result = Ade |> wt1
		result2 = +Ade |> wt2
		true
	catch
		false
	end
end

# Test classification on a 1-D dataset
A1d = j4pr.datacell([rand(100);2*rand(100)],[zeros(100);ones(100)])
Test.@test try
	wt1 = A1d |> Wclass[1]
	result = A1d |> wt1
	true
	catch
		false
	end

# Test Printer
Test.@test try
	buf = IOBuffer()

	# test show for classifier
	wt1 = A1d |> Wclass[1]
	Base.show(buf,wt1.x.data)

	# test show for estimator
	wt1 = +A1d |> Wclass[1]
	Base.show(buf,wt1.x.data)

	true #works
	catch
		false
	end
end


text:
using Cleaner
using Documenter

DocMeta.setdocmeta!(Cleaner, :DocTestSetup, :(using Cleaner); recursive=true)

makedocs(;
    modules=[Cleaner],
    authors="Ronald Gamez",
    repo="https://github.com/TheRoniOne/Cleaner.jl/blob/{commit}{path}#{line}",
    sitename="Cleaner.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://TheRoniOne.github.io/Cleaner.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Guide" => Any[
            "First steps" => "man/first_steps.md",
            "A name-changing help" => "man/name_changing.md",
            "Getting the dirt out" => "man/dirt_removal.md",
            "Workflow tips" => "man/workflow_tips.md"
        ],
        "API" => Any[
            "Types" => "lib/types.md",
            "Functions" => "lib/functions.md"
        ]
    ],
    strict = true
)

deploydocs(;
    repo="github.com/TheRoniOne/Cleaner.jl",
)

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
using Plots

function wheeler(x, a=1.5)
    return exp(-(x[1]*x[2] - a)^2 - (x[2] - a)^2);
end

# plot

heatmap(-10:0.01:25, -3:0.01:6, (x, y)->wheeler([x, y]))


text:
function grid_search(f::Function, grid::Vector)
    min_value = Inf
    arg_min_value = Inf

    for el in grid
        if f(el) < min_value
            arg_min_value = el
        end
    end

    return arg_min_value
end


text:
using Mimi
using MimiFUND
using DataFrames

# Run FUND for baseline scenario
m = MimiFUND.get_model()

# set_dimension!(m, :time, 1950:2100) # set the timeline to be 1950 to 2100

run(m)

explore(m)

p = Mimi.plot(m, :impactagriculture, :aglevel)
p = Mimi.plot(m, :impactagriculture, :agrate)
p = Mimi.plot(m, :impactagriculture, :agcost)

scc_year = 2020

update_param!(m, :climatedynamics, :climatesensitivity, 4.5)

mm = MimiFUND.get_marginal_model(m, year = scc_year)   # The additional emissions pulse will be added in the specified year

run(mm)

MimiFUND.compute_scc(m, year = scc_year, eta = 1., prtp = 0.) * 1.68

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
"""
    fakedata_gaussian(pts_per_clust::Int = 10, nclust::Int = 2, xmin = 0.0, xmax = nclust * 7.0, nx = round(Int, (xmax - xmin) * 10), σ = 1.0, amin = 1.0, amax = 2.0)

DOCSTRING

# Arguments:
- `pts_per_clust`: DESCRIPTION
- `nclust`: DESCRIPTION
- `xmin`: DESCRIPTION
- `xmax`: DESCRIPTION
- `nx`: DESCRIPTION
- `σ`: DESCRIPTION
- `amin`: DESCRIPTION
- `amax`: DESCRIPTION
"""
function fakedata_gaussian(
    pts_per_clust::Int = 10,
    nclust::Int        = 2,
    xmin               = 0.0,
    xmax               = nclust * 7.0,
    nx                 = round(Int, (xmax - xmin) * 10),
    σ                  = 1.0,
    amin               = 1.0,
    amax               = 2.0,
)
    x = LinRange(xmin, xmax, nx)
    μ = LinRange(xmin, xmax, nclust + 2)[2:end-1]
    a = LinRange(xmin, xmax, pts_per_clust + 2)[2:end-1]

    npts = pts_per_clust * nclust
    data = zeros(nx, npts)
    labels = zeros(Int, npts)

    i = 1
    for c = 1:nclust
        for n = 1:pts_per_clust
            data[:, i] = gauss_func.(a[n], μ[c], σ, x)
            labels[i] = c
            i += 1
        end
    end

    return data, labels
end

gauss_func(a, μ, σ, x) = a * exp((-(x - μ) .^ 2) / (2 * σ^2))


text:
using JuMP
using DisjunctiveProgramming

m = Model()
@variable(m, -1<=x<=10)

@constraint(m, con1, x<=3)
@constraint(m, con2, 0<=x)
@constraint(m, con3, x<=9)
@constraint(m, con4, 5<=x)

@disjunction(m,(con1,con2),con3,con4,reformulation=:CHR,name=:y)


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "Matrix.jl" begin

   # function flatten(T::Matrix3)
   @testset "flatten" begin
      @test
      @test
      @test
      @test
   end

	# function flatten(T::Matrix4)
   @testset "flatten" begin
      @test
      @test
      @test
      @test
   end

	 # function dropW(T::Matrix4)
   @testset "dropW" begin
      @test
      @test
      @test
      @test
   end

	 # function translateMatrix(vt::Point3d)
   @testset "translateMatrix" begin
      @test
      @test
      @test
      @test
   end

	 # function scaleMatrix(vs::Point3d)
   @testset "scaleMatrix" begin
      @test
      @test
      @test
      @test
   end

	 # function lookAtMatrix(eye::Point3d, center::Point3d, up::Point3d)
   @testset "lookAtMatrix" begin
      @test
      @test
      @test
      @test
   end

	 # function perspectiveMatrix(fovy::Float64, aspect::Float64, zNear::Float64, zFar::Float64)
   @testset "perspectiveMatrix" begin
      @test
      @test
      @test
      @test
   end

	 # function getLookAt(T::Matrix4)
   @testset "getLookAt" begin
      @test
      @test
      @test
      @test
   end

end


text:
using OrdinaryDiffEq
import OrdinaryDiffEq: ODEProblem
using Petri

export linreg, peakgap, paramsweep, fluxes, u₀, makeplots_seir, makeplots_seird, savedata

#this function was removed in Julia 0.7
"""    linreg(x,y)

perform a simple linear regression for scalar case. Coefficients are returned
as (intercept, scalar)
"""
linreg(x, y) = hcat(fill!(similar(x), 1), x) \ y

"""    peakgap(s1, s2)

identify the time delay between the peak in two time series. Time series should
be stored as (x(t), t) for consistency with the tuples(sol) from OrdinaryDiffEq.
"""
function peakgap(s1, s2)
    @show p1, i1 = findmax(first.(s1))
    @show p2, i2 = findmax(first.(s2))
    return s2[i2][end] - s1[i1][end]
end

"""    peakgap(sol, i::Int, j::Int)

identify the time delay between the peak in two dimensions of an ODE solution structure.
"""
function peakgap(sol, i::Int, j::Int)
    s1 = [(u[i],t) for (u,t) in tuples(sol)]
    s2 = [(u[j],t) for (u,t) in tuples(sol)]
    return peakgap(s1,s2)
end

"""    paramsweep(f::Function, m::Petri.Model, initials, tspan, params)

solve an ODEProblem for each parameter setting in params and apply the function f to the solution.

Usage:
    paramsweep(m,initials,tspan,[p1,p2,p3]) do sol
        return f(sol)
    end
"""
function paramsweep(f::Function, m::Petri.Model, initials, tspan, params)
    map(params) do p
        prob = ODEProblem(fluxes(m), initials, tspan, p)
        sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
        return f(sol)
    end
end

"""    fluxes(m::Petri.Model)

a PetriNet interpreter that computes the mass action kinetics of a petri net.

Usage: pass this to ODEProblem to set up an ODE for a given model.
"""
function fluxes(m::Petri.Model)
    S = m.S
    T = m.Δ
    nS = length(S)
    nT = length(T)
    ϕ = zeros(Float64, nT)
    f(du, u, p, t) = begin
        for (i, t) in enumerate(T)
            ins = t[1]
            # TODO: accomodate multiplicites here
            ϕ[i] = p[i]*prod(u[ins])
        end
        for i in S
            du[i] = 0
        end
        for (i, t) in enumerate(T)
            ins = t[1]
            out = t[2]
            for s in ins
                # TODO: accomodate multiplicites here
                du[s] -= ϕ[i]
            end
            for s in out
                # TODO: accomodate multiplicites here
                du[s] += ϕ[i]
            end
        end
        return du
    end
    return f
end

u₀(m::Petri.Model) = begin
    zeros(Float64, length(m.S))
end

u₀(m::Petri.Model, initialS, initialI=1) = begin
    u0=zeros(Float64, length(m.S))
    u0[1] = initialS
    u0[2] = initialI
    return u0
end

function savedata(sol::ODESolution, colnames, fname::String)
    open(fname, "w") do fp
        println(fp, "time, $colnames")
        map(tuples(sol)) do (u,t)
            print(fp, "$t")
            map(u) do x
                print(fp, ",$x")
            end
            println(fp, "")
        end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Xdmf.jl/blob/master/LICENSE

using JuliaFEM
using JuliaFEM.Testing
using LightXML

@testset "create new Xdmf object" begin
    r = Xdmf()
    expected = "<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\" Version=\"3.0\"/>"
    @test string(r.xml) == expected
end

@testset "put and get to Xdmf, low level" begin
    io = Xdmf()
    # h5
    write(io.hdf, "/Xdmf/Domain/Geometry", [1 2 3])
    @test isapprox(read(io.hdf, "/Xdmf/Domain/Geometry"), [1 2 3])
    # xml
    obj = new_child(io.xml, "Domain")
    set_attribute(obj, "Name", "Test Domain")
    obj2 = find_element(io.xml, "Domain")
    @test attribute(obj2, "Name") == "Test Domain"
end

@testset "write data to HDF, automatically generate path" begin
    xdmf = Xdmf()
    di1 = new_dataitem(xdmf, [1 2 3])
    di2 = new_dataitem(xdmf, [4 5 6])
    @test contains(content(di1), "DataItem_1")
    @test contains(content(di2), "DataItem_2")
end

@testset "Xdmf filtering" begin
    grid1 = new_element("Grid")
    add_text(grid1, "I am first grid")
    grid2 = new_element("Grid")
    add_text(grid2, "I am second grid")
    set_attribute(grid2, "Name", "Frame 2")
    grid3 = new_element("Grid")
    add_text(grid3, "I am third grid")
    grids = [grid1, grid2, grid3]

    @test content(xdmf_filter(grids, "Grid")) == "I am first grid"
    @test content(xdmf_filter(grids, "Grid[1]")) == "I am first grid"
    @test content(xdmf_filter(grids, "Grid[2]")) == "I am second grid"
    @test content(xdmf_filter(grids, "Grid[3]")) == "I am third grid"
    @test content(xdmf_filter(grids, "Grid[end]")) == "I am third grid"
    @test content(xdmf_filter(grids, "Grid[@Name=Frame 2]")) == "I am second grid"
    @test xdmf_filter(grids, "Grid[0]") == nothing
    @test xdmf_filter(grids, "Grid[4]") == nothing
    @test xdmf_filter(grids, "Grid[@Name=Frame 3]") == nothing
    @test xdmf_filter(grids, "Domain/Grid[@Name=Frame 3]") == nothing
    @test xdmf_filter(grids, "Domain") == nothing
end

@testset "XML traverse" begin
    xdmf = Xdmf()
    domain = new_child(xdmf.xml, "Domain")
    grid = new_child(domain, "Grid")
    set_attribute(grid, "CollectionType", "Temporal")
    set_attribute(grid, "GridType", "Collection")

    frame1 = new_child(grid, "Grid")
    time1 = new_child(frame1, "Time")
    set_attribute(time1, "Value", 0.0)
    X1 = new_child(frame1, "Geometry")
    set_attribute(X1, "Type", "XY")

    frame2 = new_child(grid, "Grid")
    set_attribute(frame2, "Name", "Frame 2")
    time2 = new_child(frame2, "Time")
    set_attribute(time2, "Value", 1.0)
    X2 = new_child(frame2, "Geometry")
    set_attribute(X2, "Type", "XY")

    add_child(grid, frame1)
    add_child(grid, frame2)

    dataitem = new_dataitem(xdmf, "/Domain/Grid/Grid/2/Geometry", [1.0, 2.0])
    add_child(X2, dataitem)

    println(xdmf.xml)
    @test read(xdmf, "/Domain/Grid/Grid/Time/Value") == "0.0"
    @test read(xdmf, "/Domain/Grid/Grid[2]/Time/Value") == "1.0"
    @test read(xdmf, "/Domain/Grid/Grid[end]/Time/Value") == "1.0"
    @test read(xdmf, "/Domain/Grid/Grid[@Name=Frame 2]/Time/Value") == "1.0"
    @test isapprox(read(xdmf, "/Domain/Grid/Grid[2]/Geometry/DataItem"), [1.0, 2.0])
end

@testset "write fields from different problems to Xdmf file" begin

    X = Dict(
             1 => [0.0, 0.0, 0.0],
             2 => [1.0, 0.0, 0.0],
             3 => [1.0, 1.0, 0.0],
             4 => [0.0, 1.0, 0.0],
             5 => [0.0, 0.0, 0.5],
             6 => [1.0, 0.0, 0.5],
             7 => [1.0, 1.0, 0.5],
             8 => [0.0, 1.0, 0.5],
             9 => [0.0, 0.0, 1.0],
             10 => [1.0, 0.0, 1.0],
             11 => [1.0, 1.0, 1.0],
             12 => [0.0, 1.0, 1.0])

    u = Dict()
    u[0] = Dict(
                1 => [0.0, 0.0, 0.0],
                2 => [0.0, 0.0, 0.0],
                3 => [0.0, 0.0, 0.0],
                4 => [0.0, 0.0, 0.0],
                5 => [0.0, 0.0, 0.0],
                6 => [0.0, 0.0, 0.0],
                7 => [0.0, 0.0, 0.0],
                8 => [0.0, 0.0, 0.0],
                9 => [0.0, 0.0, 0.0],
                10 => [0.0, 0.0, 0.0],
                11 => [0.0, 0.0, 0.0],
                12 => [0.0, 0.0, 0.0])
    u[1] = Dict(
                1 => [0.0, 0.0, 0.0],
                2 => [0.0, 0.0, 0.0],
                3 => [0.0, 0.0, 0.0],
                4 => [0.0, 0.0, 0.0],
                5 => [0.0, 0.0, -0.1],
                6 => [0.0, 0.0, -0.1],
                7 => [0.0, 0.0, -0.1],
                8 => [0.0, 0.0, -0.1],
                9 => [0.0, 0.0, -0.2],
                10 => [0.0, 0.0, -0.2],
                11 => [0.0, 0.0, -0.2],
                12 => [0.0, 0.0, -0.2])

    T = Dict(
             1 => 10.0,
             2 => 10.0,
             3 => 10.0,
             4 => 10.0,
             5 => 20.0,
             6 => 20.0,
             7 => 20.0,
             8 => 20.0,
             9 => 30.0,
             10 => 30.0,
             11 => 30.0,
             12 => 30.0)

    rf = Dict()
    rf[0] = Dict(
                 1 => [0.0, 0.0, 0.0],
                 2 => [0.0, 0.0, 0.0],
                 3 => [0.0, 0.0, 0.0],
                 4 => [0.0, 0.0, 0.0])
    rf[1] = Dict(
                 1 => [0.0, 0.0, 1.0],
                 2 => [0.0, 0.0, 1.0],
                 3 => [0.0, 0.0, 1.0],
                 4 => [0.0, 0.0, 1.0])

    e1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
    e2 = Element(Hex8, [5, 6, 7, 8, 9, 10, 11, 12])
    e3 = Element(Quad4, [1, 2, 3, 4])
    update!([e1, e2, e3], "geometry", X)
    update!([e1, e2, e3], "displacement", 0.0 => u[0])
    update!([e1, e2, e3], "displacement", 1.0 => u[1])
    update!([e1, e2, e3], "temperature", T)
    update!(e3, "reaction force", 0.0 => rf[0])
    update!(e3, "reaction force", 1.0 => rf[1])

    p1 = Problem(Elasticity, "lower", 3)
    p1.elements = [e1]
    p2 = Problem(Elasticity, "upper", 3)
    p2.elements = [e2]
    p3 = Problem(Dirichlet, "bc", 3, "displacement")
    p3.elements = [e3]

    xdmf = Xdmf()
    xdmf.format = "XML"
    update_xdmf!(xdmf, p1, 0.0, ["displacement", "temperature"])
    update_xdmf!(xdmf, p2, 0.0, ["displacement"])
    update_xdmf!(xdmf, p3, 0.0, ["reaction force"])
    update_xdmf!(xdmf, p1, 1.0, ["displacement", "temperature"])
    update_xdmf!(xdmf, p2, 1.0, ["displacement"])
    update_xdmf!(xdmf, p3, 1.0, ["reaction force"])
    @test read(xdmf, "/Domain/Grid/Grid/Time/Value") == "0.0"
    @test read(xdmf, "/Domain/Grid/Grid[2]/Time/Value") == "1.0"
end


text:
# Initializing Python (surprisingly complicated; see also deps/build.jl)

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

# Base.with_env in Julia 0.3, withenv in Julia 0.4 (#10914)
if VERSION >= v"0.4.0-dev+4392"
    with_env(f::Function, key::AbstractString, val) = withenv(f, key=>val)
else
    const with_env = Base.with_env
end

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

function __init__()
    already_inited = 0 != ccall((@pysym :Py_IsInitialized), Cint, ())

    if !already_inited
        if !isempty(PYTHONHOME)
            if pyversion_build.major < 3
                ccall((@pysym :Py_SetPythonHome), Void, (Cstring,), PYTHONHOME)
            else
                ccall((@pysym :Py_SetPythonHome), Void, (Cwstring,), PYTHONHOME)
            end
        end
        if !isempty(pyprogramname)
            if pyversion_build.major < 3
                ccall((@pysym :Py_SetProgramName), Void, (Cstring,), pyprogramname)
            else
                ccall((@pysym :Py_SetProgramName), Void, (Cwstring,), pyprogramname)
            end
        end
        ccall((@pysym :Py_InitializeEx), Void, (Cint,), 0)
    end

    # cache the Python version as a Julia VersionNumber
    global const pyversion = convert(VersionNumber,
                                     split(bytestring(ccall(@pysym(:Py_GetVersion),
                                                            Ptr{Uint8}, ())))[1])
    if pyversion_build.major != pyversion.major
        error("PyCall built with Python $pyversion_build, but now using Python $pyversion; ",
              "you need to relaunch Julia and re-run Pkg.build(\"PyCall\")")
    end

    global const inspect = pyimport("inspect")
    global const builtin = pyimport(pyversion.major < 3 ? "__builtin__" : "builtins")

    pyexc_initialize() # mappings from Julia Exception types to Python exceptions

    # Python has zillions of types that a function be, in addition
    # to the FunctionType in the C API.  We have to obtain these
    # at runtime and cache them in globals
    types = pyimport("types")
    global const BuiltinFunctionType = types["BuiltinFunctionType"]
    global const TypeType = pybuiltin("type")
    global const MethodType = types["MethodType"]
    global const MethodWrapperType = pytypeof(PyObject(PyObject[])["__add__"])
    global const ufuncType = try
        pyimport("numpy")["ufunc"]
    catch
        PyNULL() # NumPy not available
    end

    # cache Python None -- PyPtr, not PyObject, to prevent it from
    # being finalized prematurely on exit
    global const pynothing = @pyglobalobj(:_Py_NoneStruct)

    # xrange type (or range in Python 3)
    global const pyxrange = @pyglobalobj(:PyRange_Type)

    # cache ctypes.c_void_p type and function if available
    vpt, pvp = try
        (pyimport("ctypes")["c_void_p"],
         p::Ptr -> pycall(c_void_p_Type, PyObject, @compat UInt(p)))
    catch # fallback to CObject
        (@pyglobalobj(:PyCObject_FromVoidPtr),
         p::Ptr -> PyObject(ccall(pycobject_new, PyPtr, (Ptr{Void}, Ptr{Void}), p, C_NULL)))
    end
    global const c_void_p_Type = vpt
    global const py_void_p = pvp

    # traceback.format_tb function, for show(PyError)
    global const format_traceback = pyimport("traceback")["format_tb"]

    # all cfunctions must be compiled at runtime
    global const jl_Function_call_ptr =
        cfunction(jl_Function_call, PyPtr, (PyPtr,PyPtr,PyPtr))
    global const pyio_repr_ptr = cfunction(pyio_repr, PyPtr, (PyPtr,))
    global const pyjlwrap_dealloc_ptr = cfunction(pyjlwrap_dealloc, Void, (PyPtr,))
    global const pyjlwrap_repr_ptr = cfunction(pyjlwrap_repr, PyPtr, (PyPtr,))
    global const pyjlwrap_hash_ptr = cfunction(pyjlwrap_hash, Uint, (PyPtr,))
    global const pyjlwrap_hash32_ptr = cfunction(pyjlwrap_hash32, Uint32, (PyPtr,))

    # similarly, any MethodDef calls involve cfunctions
    global const jl_TextIO_methods = make_io_methods(true)
    global const jl_IO_methods = make_io_methods(false)
    global const jl_IO_getset = PyGetSetDef[
            PyGetSetDef("closed", jl_IO_closed)
            PyGetSetDef("encoding", jl_IO_encoding)
            PyGetSetDef()
    ]

    # PyMemberDef stores explicit pointers, hence must be initialized in __init__
    global const pyjlwrap_members =
        PyMemberDef[ PyMemberDef(pyjlwrap_membername,
                                 T_PYSSIZET, sizeof_PyObject_HEAD, READONLY,
                                 pyjlwrap_doc),
                     PyMemberDef(C_NULL,0,0,0,C_NULL) ]

    init_datetime()
    pyjlwrap_init()

    global const jl_FunctionType = pyjlwrap_type("PyCall.jl_Function",
                                                 t -> t.tp_call =
                                                 jl_Function_call_ptr)

    if !already_inited
        # some modules (e.g. IPython) expect sys.argv to be set
        if pyversion.major < 3
            argv_s = bytestring("")
            argv = unsafe_convert(Ptr{Uint8}, argv_s)
            ccall(@pysym(:PySys_SetArgvEx), Void, (Cint,Ptr{Ptr{Uint8}},Cint), 1, &argv, 0)
        else
            argv_s = Cwchar_t[0]
            argv   = unsafe_convert(Ptr{Cwchar_t}, argv_s)
            ccall(@pysym(:PySys_SetArgvEx), Void, (Cint, Ptr{Ptr{Cwchar_t}}, Cint), 1, &argv, 0)
        end

        # Some Python code checks sys.ps1 to see if it is running
        # interactively, and refuses to be interactive otherwise.
        # (e.g. Matplotlib: see PyPlot#79)
        if isinteractive()
            let sys = pyimport("sys")
                if !haskey(sys, "ps1")
                    sys["ps1"] = ">>> "
                end
            end
        end
    end
end


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# The `main.jl` file is where `openBF` is implemented. Only
# [`openBF`](openBF.html) library is needed to be imported.
# Until the official `openBF` Julia `Pkg` is created, the
# library is loaded locally.
# push!(LOAD_PATH, "src/")
using openBF
# reload("openBF")

# The project name must be specified by the user when launching
# the simulation.
# `openBF` takes care of all the pre-processing and solution steps. A
# simulation
# is started with the following command
#
#     julia main.jl project_name
#
# `project_name` is a string used to initialise all the output files.
project_name = ARGS[1]
verbose = true
no_out = false
no_inputs = true

# In `main.jl` all functions from `openBF` library are called using the dot
# notation: `library.function(parameters)`. Function
# [`projectPreamble`](initialise.html#projectPreamble) checks for all
# the files
# needed by the simulation and warns whether any is missing, it prints
# the initial `openBF` logo, and makes the directory structure to save
# temporary and final results.
openBF.projectPreamble(project_name, no_out, no_inputs)

# `project_constants.jl` file is user defined and must be in the same folder
# where the simulation is started (see [tutorial](../index.html#tutorial)
# page). Here `project_constants.jl` content is loaded in memory for further
# use.
if verbose
    println("Load project $project_name files")
end
include(join([project_name, "_constants.jl"]))

# ### Arterial system
# Arterial model and structure is encoded in a `.csv` file user defined. It
# must be in the same folder in which the simulation is started.
# [`readModelData`](initialise.html#readModelData) reads the
# `.csv` file and fill a 2D matrix with all the informations.
model = openBF.readModelData(join([project_name, ".csv"]))

# Data from `project.csv`, `project_constants.jl`, and `project_inlet.dat` (if
# specified) are used to create instances of [`BTypes`](BTypes.html) data
# structures by [`loadGlobalConstants`](initialise.html#loadGlobalConstants).
heart, blood_prop, total_time = openBF.loadGlobalConstants(project_name,
  inlet_BC_switch, inlet_type, cycles, rho, mu, gamma_profile)

# The arterial tree is represented as a graph by means of `Graphs` library (
# see [grafo.jl](grafo.html) for a detailed example).
# `simple_graph` creates an empty `GenericGraph` with a user defined number of
# nodes.
#
# ----------------------------------------------------------------------------
# Parameters
# ----------- ----------------------------------------------------------------
# `nodes`     `::Int64` number of total nodes in the graph. Here an estimate
#             of this number is taken by counting `model` matrix rows (one
#             for each vessel) and adding one extra node.
# ----------------------------------------------------------------------------
# Returns
# ----------- ----------------------------------------------------------------
# `graph`     `::GenericGraph` empty graph data structure.
# ----------------------------------------------------------------------------
#
#
# grafo = LightGraphs.DiGraph(length(model[:,1])+1)

# Two data structures are used to describe the arterial system. One is a
# collection of [`BTypes`](BTypes.html)`.Vessel` instances, and the second
# is the `grafo` structure. `vessels` collection is filled by using
# [`initialiseVessel`](initialise.html#initialiseVessel). The first
# vessel is inserted by hand.
if verbose
    println("Build arterial network")
end
vessels = [openBF.initialiseVessel(model[1,:], 1, heart, blood_prop,
  initial_pressure, Ccfl)]
edge_list = zeros(Int8, length(model[:,1]), 3)
edge_list[1,1] = vessels[1].ID
edge_list[1,2] = vessels[1].sn
edge_list[1,3] = vessels[1].tn

# The graph is built with `Graphs.add_edge!` function.
#
# ----------------------------------------------------------------------------
# Parameters
# ----------- ----------------------------------------------------------------
# `grafo`     `::GenericGraph` arterial system graph structure.
#
# `sn`        `::Int64` edge source node stored in `Vessel` structure.
#
# `tn`        `::Int64` edge terminal node.
# ----------------------------------------------------------------------------
# LightGraphs.add_edge!(grafo, vessels[1].sn, vessels[1].tn)

# The model matrix is read iteratively starting from the second row, the first
# row contains column headers.
for i in 2:length(model[:,1])

  # Each new `vessel` is instantiated and `push!`ed at the bottom of
  # `vessels` collection.
  push!(vessels, openBF.initialiseVessel(model[i,:], i, heart, blood_prop,
    initial_pressure, Ccfl))
  edge_list[i,1] = vessels[i].ID
  edge_list[i,2] = vessels[i].sn
  edge_list[i,3] = vessels[i].tn
  # LightGraphs.add_edge!(grafo, vessels[end].sn, vessels[end].tn)
end

# # `edge_list` is a list of all the edges in `grafo`
# edge_list = collect(LightGraphs.edges(grafo))
# edge_map = Dict{LightGraphs.SimpleGraphs.SimpleEdge, Int}(e=>i
#   for (i, e) in enumerate(edge_list))
# node_map = Dict{Tuple{Int, Int}, Int}()
# n_i = 1
# for e in edge_list
#   s = LightGraphs.src(e)
#   t = LightGraphs.dst(e)
#   node_map[(s,t)] = n_i
#   n_i += 1
# end

# Before starting the main loop the counter`current_time` is set to zero. It
# will be updated to keep track of time within the simulation.
if verbose
    println("Start simulation \n")
end
current_time = 0

# In order to show the progress bar an initial estimate of the total running
# time is needed. Thus, $\Delta t$ is calculated with system initial
# conditions
# and used to compute the number or total iterations before the end of the
# simulation. See [ProgressMeter](https://github.com/timholy/ProgressMeter.jl)
# documentation for `Progress` options.
dts  = zeros(Float64, length(edge_list[:,1]))
dt = openBF.calculateDeltaT(vessels, dts)

# # ### Venous system
# # Same as the arterial system
# v_model = join([project_name, "_veins.csv"])
# if (isfile(v_model)) == true
#   model_v = openBF.readModelData(join([project_name, "_veins.csv"]))
#
#   grafo_v = Graphs.simple_graph(length(model_v[:,1])+1)
#
#   vessels_v = [openBF.initialiseVessel(model_v[1,:], 1, heart, blood_prop,
#     initial_pressure, Ccfl)]
#
#   Graphs.add_edge!(grafo_v, vessels_v[1].sn, vessels_v[1].tn)
#
#   for i in 2:length(model_v[:,1])
#
#     push!(vessels_v, openBF.initialiseVessel(model_v[i,:], i, heart, blood_prop,
#       initial_pressure, Ccfl))
#
#     Graphs.add_edge!(grafo_v, vessels_v[end].sn, vessels_v[end].tn)
#   end
#
#   edge_list_v = Graphs.edges(grafo_v)
#
#   dts_v  = zeros(Float64, length(edge_list_v))
#   dt_v = openBF.calculateDeltaT(vessels_v, dts_v)
#
#   dt = minimum([dt, dt_v])
#
#   venous_model = true
# else
#   venous_model = false
# end

# The simulation is ran in a `while` loop to be ended whether the simulation
# reached convergence or a user defined finish time.
passed_cycles = 0
if verbose
    @printf("Solving cardiac cycle no: %d", passed_cycles + 1)
    tic()
end
counter = 1
counter_v = 0
jump = 100
timepoints = linspace(0, heart.cardiac_T, jump)
while true
  # At the beginning of each time step the $\Delta t$ is computed with
  # [`calculateDeltaT`](godunov.html#calculateDeltaT). This is because
  # during the simulation the local properties of each vessel will change and,
  # as a consequence, also the $\Delta t$ will change. The `current_time` is
  # then calculated by incrementing it with the new `dt`.
  dt = openBF.calculateDeltaT(vessels, dts)

  # if venous_model
  #   dt_v = openBF.calculateDeltaT(vessels_v, dts_v)
  #   dt = minimum([dt, dt_v])
  # end

  current_time += dt

  # [`solveModel`](godunov.html#solveModel) reads the `grafo` object
  # and runs the solver for each part of it. This function can distinguish
  # between inlet, bifurcation, conjunction, anastomosis, and outlet.
  #Solve arteries
  openBF.solveModel(vessels, heart, edge_list, blood_prop, dt, current_time)

  # [`updateGhostCells`](boundary_conditions.html#updateGhostCells)
  # updates all vessels ghost cells after the solver ends one iteration.
  openBF.updateGhostCells(vessels)

  # All quantities in each vessel are stored by
  # [`saveTempData`](IOutils.html#saveTempData) in `.temp` files until the
  # end of the cardiac cycle.
  if current_time >= timepoints[counter]
    openBF.saveTempData(current_time, vessels)
    counter += 1
  end

  # #Solve veins
  # if venous_model
  #   openBF.solveModel(grafo_v, edge_list_v, vessels_v,
  #                     grafo, edge_list, vessels,
  #                     blood_prop, dt, current_time)
  #   openBF.updateGhostCells(vessels_v)
  #
  #   if counter_v == 100
  #     openBF.saveTempData(current_time, vessels_v)
  #     counter_v = 0
  #   else
  #     counter_v += 1
  #   end
  #
  # end

  # Every time a cardiac cycle has been simulated, this condition returns
  # `true` and data from `.temp` files are transferred to `.out` files (
  # see [IOutils.jl](IOutils.html)).

  if (current_time - heart.cardiac_T*passed_cycles) >= heart.cardiac_T &&
      (current_time - heart.cardiac_T*passed_cycles + dt) > heart.cardiac_T

      openBF.closeTempFiles(vessels)

      err = openBF.checkConvergence(edge_list, vessels, passed_cycles)
    #   println("Iteration: ", passed_cycles, " Error: ", err[1:5],"%")
      if verbose
          @printf(" - Error = %4.2f%%\n", err)
      end

      openBF.transferLastToOut(vessels)
      openBF.openCloseLastFiles(vessels)
      openBF.transferTempToLast(vessels)
      openBF.openTempFiles(vessels)

      if err <= 5.
          break
      end

      # if venous_model
      #   openBF.closeTempFiles(vessels_v)
      #   openBF.transferLastToOut(vessels_v)
      #   openBF.openCloseLastFiles(vessels_v)
      #   openBF.transferTempToLast(vessels_v)
      #   openBF.openTempFiles(vessels_v)
      # end

    passed_cycles += 1
    if verbose
        @printf("Solving cardiac cycle no: %d", passed_cycles + 1)
    end

    timepoints += heart.cardiac_T
    counter = 1
    # # When at least 3 cardiac cycles have been simulated, waveforms are
    # # checked for convergence.
    # if passed_cycles >= 3
    #
    #   # The error is computed for all vessels in the system by
    #   # [`checkAllQuantities`](check_convergence.html#checkAllQuantities)
    #   # function.
    #   #
    #   # err = openBF.checkAllQuantities(vessels, passed_cycles, 1000)
    #   err = openBF.checkConvergence(vessels, 5.)
    #
    #   # The convergence is reached when the difference (the error)
    #   # between two consecutive waveforms is less than 5%. In this case, the
    #   # main loop is exited.
    #   if err < 5.
    #     println("\nConverged in $passed_cycles cycles, end!")
    #     break
    #   end
    # end

    # openBF.openCloseLastFiles(vessels)
    # openBF.transferTempToLast(vessels)
    # openBF.openTempFiles(vessels)
  end

  # When the `current_time` is equal to the `total_time` defined by the user
  # exit from `while` loop. This would also mean that an error smaller than 5%
  # has not been achieved. The main is exited by raising an error containing
  # the error value.
  if current_time >= total_time
    if verbose
        println("Not converged after $passed_cycles cycles, End!")
    end
    break
  end
end
if verbose
@printf "\n"
    toc()
end

# Make sure that data from `.temp` files are transferred.
openBF.closeTempFiles(vessels)
openBF.transferTempToOut(vessels)
# run(`sh cleaner.sh`)
# run(`rm cleaner.sh`)

if no_out == true
    # rm("*.out")
    cleanOuts(vessels)
    cleanTemps(vessels)
    # rm("*.temp")
end

if no_inputs == true
    rm("$project_name.csv")
    rm("$project_name\_inlet.dat")
    rm("$project_name\_constants.jl")
end

# rm("appender.sh")

# if venous_model
#   openBF.closeTempFiles(vessels_v)
#   openBF.transferTempToOut(vessels_v)
# end

cd("..")
# run(`rm main.jl`)


text:
include(joinpath("src", "axelrodnet.jl"))

for sc in 0:3:30
    Axelrodnet.run_experiment(
        experiment_name="grid_sc" * string(sc),
        agentcount=100,
        n_iter=5_000_000,
        nettopology="grid",
        networkprops=Dict("grid_height" => 10),
        stubborncount=sc,
        rndseed=1,
        n_replicates=10,
        export_every_n=1_000_000,
        exportdata=true,
        exportmode="default",
        keep_rawnet=false
    )
end

for sc in 0:3:30, beta in 0.05:0.05:0.2
    Axelrodnet.run_experiment(
        experiment_name="wattsstrogatz_sc" * string(sc) * "_beta" * lpad(trunc(Int, beta * 100), 2, "0"),
        agentcount=100,
        n_iter=5_000_000,
        nettopology="watts_strogatz",
        networkprops=Dict("k" => 10, "beta" => beta),
        stubborncount=sc,
        rndseed=1,
        n_replicates=10,
        export_every_n=1_000_000,
        exportdata=true,
        exportmode="default",
        keep_rawnet=false
    )
end

for sc in 0:3:30
    Axelrodnet.run_experiment(
        experiment_name="barabasialbert_sc" * string(sc),
        agentcount=100,
        n_iter=5_000_000,
        nettopology="barabasi_albert",
        networkprops=Dict("m0" => 10),
        stubborncount=sc,
        rndseed=1,
        n_replicates=10,
        export_every_n=1_000_000,
        exportdata=true,
        exportmode="default",
        keep_rawnet=false
    )
end

Axelrodnet.aggregate_all_experiments("experiments", write_data=true)


text:
module LinearEnrollment
    using CurricularAnalyticsSimulationLibrary: Student

    function enroll!(currentTerm, simulation, max_credits)
        studentProgress = simulation.studentProgress
        terms = simulation.curriculum.terms

        for (termnum, term) in enumerate(terms)
            # Itterate through courses
            for course in term.courses
                # Clear the array of enrolled students for the course
                course.students = Student[]

                for student in simulation.enrolledStudents
                    # Enroll in coreqs first
                    for coreq in course.coreqs
                        if canEnroll(student, coreq, studentProgress, max_credits, currentTerm)
                            # Enroll the student in the course
                            push!(course.students, student)

                            # Increment the course's enrollment counters
                            course.enrolled += 1
                            course.termenrollment[currentTerm] += 1

                            # Increse the student's term credits
                            student.termcredits += course.credits
                        end
                    end

                    # Determine wheter the student can be enrolled in the current course.
                    if canEnroll(student, course, studentProgress, max_credits, currentTerm)

                        # Enroll the student in the course
                        push!(course.students, student)

                        # Increment the course's enrollment counters
                        course.enrolled += 1
                        course.termenrollment[currentTerm] += 1

                        # Increse the student's term credits
                        student.termcredits += course.credits
                    end
                end
            end
        end
    end

    # Function that determines wheter a student can enroll in a course
    function canEnroll(student, course, studentProgress, max_credits, term)
        # Find the prereq ids of the current course
        prereqIds = map(x -> x.id, course.prereqs)

        !in(student, course.students) &&
        (length(course.prereqs) == 0 || sum(studentProgress[student.id, prereqIds]) == length(course.prereqs)) &&       # No Prereqs or the student has completed them
        studentProgress[student.id, course.id] == 0.0 &&                                                                # The student has not already completed the course
        student.termcredits + course.credits <= max_credits &&                                                          # The student will not exceed the maximum number of credit hours
        course.termReq <= term &&                                                                                # The student must wait until the term req has been met
        enrolledInCoreqs(student, course, studentProgress)                                                              # The student is enrolled in or has completed coreqs
    end

    # Determines whether a student is enrolled in or has completed coreqs for a given course
    function enrolledInCoreqs(student, course, studentProgress)
        enrolled = true

        for coreq in course.coreqs
            enrolled =  enrolled && (in(student, coreq.students) || studentProgress[student.id, coreq.id] == 1.0)
        end

        return enrolled
    end
end

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
"""
    refinedldiv(A::AbstractMatrix{T},F,B::AbstractVecOrMat{T}) => X, converged

Compute an accurate solution to a linear system ``A X = B.``

`F` must be a factorization object of `A` (typically produced by `lu`
or `cholesky`). Crucially, `eltype(F)` may have lower precision than
`T`.  Uses iterative refinement with intermediates of `eltype` `T` so
that the backward error is of order `eps(T)`. Note: this function
allocates several working matrices.

# Example
```julia
julia> A = rand(BigFloat,n,n); b = rand(BigFloat,n); F = lu(Float64.(A));
julia> x = refinedldiv(A,F,b)
```
Note that with the high precision of `BigFloat` full convergence is
often impossible, but the result is still quite accurate.

# Optional arguments
- `maxiter::Int`: default 10.
- `tol::Real`: tolerance for relative residual norm, scaled to `eps(T)`
- `verbose::Bool`: whether to report on residuals.
"""
function refinedldiv end

function refinedldiv(A::AbstractMatrix{DT}, F::Factorization{T},
                B::AbstractVecOrMat{DT};
                maxiter=10, tol=10, verbose=false
                ) where {DT <: AbstractFloat, T <: AbstractFloat}
    nrhs = size(B,2)
    n = size(A,1)
    cvtok = true
    anorm = opnorm(A,Inf)
    bnorms = [norm(view(B,:,j),Inf) for j=1:nrhs]
    tolf = tol * eps(DT)
    local X
    try
        X = T.(B)
    catch
        cvtok = false
    end
    cvtok || throw(ArgumentError("unable to convert to "
                                 * "designated narrow type $T"))
    ldiv!(F,X)
    Xd = DT.(X)
    Rd = B - A * Xd
    xnorms = [norm(view(Xd,:,j),1) for j=1:nrhs]
    # ω₂ from Arioli, Duff, Ruiz SIAM J.Matrix Anal.Appl. 13, 138 (1992)
    relresnorms = [norm(view(Rd,:,j),Inf) / (anorm * xnorms[j] + bnorms[j])
                   for j=1:nrhs]
    verbose && println("initial max rel resid norm: ",maximum(relresnorms))
    if all(relresnorms .< tolf)
        return Xd, true
    end
    for iter=1:maxiter
        R = T.(Rd)
        ldiv!(F,R)
        Xd .+= DT.(R)
        Rd = B - A * Xd

        relresnorms = [norm(view(Rd,:,j),Inf) / (anorm * xnorms[j] + bnorms[j])
                   for j=1:nrhs]
        verbose && println("iter $iter max rel resid norm: ",maximum(relresnorms))
        if all(relresnorms .< tolf)
            return Xd, true
        end
    end
    return Xd, false
end

max_stars_repo_name: RalphAS/DoubleBLAS


text:
module StatisticalRethinkingPlots

using Reexport
using Requires

@reexport using StatsPlots, Plots, LaTeXStrings

using Parameters
using KernelDensity
using StatisticalRethinking
using Distributions

using DocStringExtensions: SIGNATURES, FIELDS, TYPEDEF

function __init__()
  @require Turing="fce5fe82-541a-59a6-adf8-730c64b5f9a0" include("require/turing/turing_plots.jl")
  @require StanSample="c1514b29-d3a0-5178-b312-660c88baa699" include("require/stan/stan_plots.jl")
end

include("pairsplot.jl")
include("pk_utilities.jl")
include("plot_bounds.jl")
include("plot_density_interval.jl")
include("plot_lines.jl")
include("plot_models.jl")
include("trankplot.jl")
include("plot_model_coef.jl")
include("plot_logistic_coef.jl")
include("pk_plot.jl")
include("coeftab_plot.jl")

end