text:
# # 1. Basic viscous flow
# In this example, we will simulate various examples of a basic *unbounded* viscous flow---a flow
# without boundaries. Our initial condition will be a distribution of *vorticity*.

#md # ```@setup compileplots
#md # using Plots
#md # ```

using ViscousFlow
#-
using Plots

#=
### The basic steps
To carry out any simulation in `ViscousFlow`, we need to carry out a few basic steps:
* **Specify the problem**: Set the Reynolds number, the free stream, and any other problem parameters
* **Discretize**: Set up a solution domain, choose the grid Reynolds number and the critical time step limits
* **Construct the system structure**: Create the operators that will be used to perform the simulation
* **Initialize**: Set the initial flow field and initialize the integrator
* **Solve**: Solve the flow field
* **Examine**: Examine the results

We will go through all of these here. For the examples we will carry out in this notebook,
the first two steps need only be carried out once.

The package utilizes parameters arranged into a dictionary, or `Dict`.
The dictionary associates names (or *keys*) with values. Most of these
key have standardizes names. We
will describe the entries in this as we go, adding them in one at a time for
expositional purposes. We initializing it first and call it `my_params`
for this example.
=#
my_params = Dict()

#=
### Problem specification
We will set the Reynolds number to be 200 and no free stream. The
Reynolds number key is "Re". If we had a free stream, then we would
set this with two keys: "freestream speed" and "freestream angle".
They default to zero if we don't put anything into the dictionary.
But we always have to put in the Reynolds number.
=#
my_params["Re"] = 200

#=
### Discretize
We will set up a domain from x = -2 to x = 2, and y = -2 to y = 2. The Reynolds number helps us
determine the grid spacing and time step size. To set these, we set a target *grid Reynolds
number*, with the key "grid Re". We will set this to 4 here; if we don't set it, it
defaults to 2. Note that the choice we make here is a compromise:
* smaller grid Reynolds number means smaller grid spacing, but slower simulations
* larger grid Reynolds number means less accurate results
=#
xlim = (-2.0,2.0)
ylim = (-2.0,2.0)
my_params["grid Re"] = 4.0

#=
Then we set up the grid
=#
g = setup_grid(xlim,ylim,my_params)

# ### Construct the system structure
# This part is easy - you supply the parameters you have just set up.
# It returns a structure with all of the necessary mathematical operators:
sys = viscousflow_system(g,phys_params=my_params);

# This is now ready to solve any unbounded viscous flow problems.
# Now, we will solve a few different problems to see how it works.

#=
## A basic example: the Lamb-Oseen vortex
This example starts with a single vortex with a Gaussian distribution of vorticity. To generate this, will use the `SpatialGaussian` function:
The command below creates a Gaussian with radius σ at (0,0) with strength 1.
=#
σ = 0.2
x0 = 0.0
y0 = 0.0
A = 1
gauss = SpatialGaussian(σ,x0,y0,A)

#=
### Initialize
Now, we create an instance of this vorticity distribution on the grid by
calling [`init_sol`](@ref) with this vortex.
=#
u0 = init_sol(gauss,sys)

#=
We use this initial condition to initialize the **integrator**. The integrator is the structure that
holds all of our solution and operator information. With it, we can start the simulation, restart
the simulation later, etc. We specify a range of time over which to advance the solution.

**Note**: This range need only be large enough to contain the whole interval of time we wish to simulate. It does not need to be chosen very precisely.

**Note 2**: There is no need to restart the problem! We can keep applying the `step!` function below as long we need.
=#
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)

# ### Solve
# Now we are ready to solve the problem. Let's advance the solution to $t = 1$:
step!(integrator,1.0)

# We can see now that the solution has been advanced in time:
integrator

#=
### Examine
Let's examine the results. It is important to show a few different ways that we can do this.
The most straightforward way is to just look at the flow fields at the current state of the
integrator. For example, to get the current velocity field, type `velocity(integrator)`. We can
do the same for `vorticity`, `streamfunction` (the streamlines), `scalarpotential`, `convective_derivative`,
and `pressure`.
=#

# We will look at some of these at the current state:
plot(
 plot(vorticity(integrator),sys,title="Vorticity"),
 plot(streamfunction(integrator),sys,title="Streamlines"),
 plot(pressure(integrator),sys,title="Pressure"),
 plot(velocity(integrator),sys))

# For this problem, we can compare with the exact solution. The exact solution is also a Gaussian,
# but with a radius $\sqrt{\sigma^2+4t/Re}$
oseen_exact(t) = SpatialGaussian(sqrt(σ^2+4*t/my_params["Re"]),x0,y0,A)
exactsol(t) = init_sol(oseen_exact(t),sys)

#-
plot(vorticity(integrator)[:,104],label="Numerical")
plot!(vorticity(exactsol(integrator.t),sys,integrator.t)[:,104],label="Exact")
plot!(title=string("Vorticity at t = ",round(integrator.t,digits=2)))

#=
## Second example: co-rotating vortices
The previous example is not very exciting, because the convection of the flow is simply circular.
The next example is more interesting, because we will start with two vortices that influence each
other's motion:
=#
σ = 0.1
x01, y01 = 0.5, 0.0
x02, y02 = -0.5, 0.0
A = 1
twogauss = SpatialGaussian(σ,x01,y01,A) + SpatialGaussian(σ,x02,y02,A)

#=
### Initialize
Now, we create an instance of this vorticity distribution on the grid.
=#
u0 = init_sol(twogauss,sys)

#-
plot(vorticity(u0,sys,0.0),sys)

#-
tspan = (0.0,8.0)
integrator = init(u0,tspan,sys)

# ### Solve
# Now we are ready to solve the problem. Let's advance the solution to $t = 8$:
step!(integrator,8.0)

#=
### Examine
In this case, it is best to view the results as an animation. In the previous example, we just
looked at the final state of the integrator. Here, to animate, we will make use of the solution
history that is held by the integrator, `integrator.sol`. Let's create an alias for this to shorten
our commands:
=#
sol = integrator.sol;

# The vortices orbit each other and then eventually merge together. If we wish to make a nice
# figure, we can arrange snapshots on a grid:
plt = plot(layout = (2,4), size = (800, 400), legend=:false)
tsnap = 0.0:1.0:7.0
for (i, t) in enumerate(tsnap)
    plot!(plt[i],vorticity(sol,sys,t),sys,levels=range(0.1,5,length=31))
end
savefig(plt,"CoRotating.pdf")
plt

#=
If you wish to animate the solution, e.g., plotting the vorticity every 5 steps
you can use

    @gif for t in sol.t
       plot(vorticity(sol,sys,t),sys)
    end every 5
=#

#=
**Try other examples!**
* Make one or both of the vortices into elliptical shapes
* Make one stronger than the other
* Add other vortices into the initial distribution
=#


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module QueryTests

using Test
using Dates

Jan = DatesPlus.DateTime(2013, 1, 1) # Tuesday
Feb = DatesPlus.DateTime(2013, 2, 2) # Saturday
Mar = DatesPlus.DateTime(2013, 3, 3) # Sunday
Apr = DatesPlus.DateTime(2013, 4, 4) # Thursday
May = DatesPlus.DateTime(2013, 5, 5) # Sunday
Jun = DatesPlus.DateTime(2013, 6, 7) # Friday
Jul = DatesPlus.DateTime(2013, 7, 7) # Sunday
Aug = DatesPlus.DateTime(2013, 8, 8) # Thursday
Sep = DatesPlus.DateTime(2013, 9, 9) # Monday
Oct = DatesPlus.DateTime(2013, 10, 10) # Thursday
Nov = DatesPlus.DateTime(2013, 11, 11) # Monday
Dec = DatesPlus.DateTime(2013, 12, 11) # Wednesday
monthnames = ["January", "February", "March", "April",
                "May", "June", "July", "August", "September",
                "October", "November", "December"]
daysofweek = [DatesPlus.Tue, DatesPlus.Sat, DatesPlus.Sun, DatesPlus.Thu, DatesPlus.Sun, DatesPlus.Fri,
              DatesPlus.Sun, DatesPlus.Thu, DatesPlus.Mon, DatesPlus.Thu, DatesPlus.Mon, DatesPlus.Wed]
dows = ["Tuesday", "Saturday", "Sunday", "Thursday", "Sunday", "Friday",
        "Sunday", "Thursday", "Monday", "Thursday", "Monday", "Wednesday"]
daysinmonth = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
@testset "Name functions" begin
    for (i, dt) in enumerate([Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec])
        @test DatesPlus.month(dt) == i
        @test DatesPlus.monthname(dt) == monthnames[i]
        @test DatesPlus.monthname(i) == monthnames[i]
        @test DatesPlus.monthabbr(dt) == monthnames[i][1:3]
        @test DatesPlus.monthabbr(i) == monthnames[i][1:3]
        @test DatesPlus.dayofweek(dt) == daysofweek[i]
        @test DatesPlus.dayname(dt) == dows[i]
        @test DatesPlus.dayname(DatesPlus.dayofweek(dt)) == dows[i]
        @test DatesPlus.dayabbr(dt) == dows[i][1:3]
        @test DatesPlus.dayabbr(DatesPlus.dayofweek(dt)) == dows[i][1:3]
        @test DatesPlus.daysinmonth(dt) == daysinmonth[i]
    end
end
@testset "Customizing locale" begin
    DatesPlus.LOCALES["french"] = DatesPlus.DateLocale(
        ["janvier", "février", "mars", "avril", "mai", "juin",
         "juillet", "août", "septembre", "octobre", "novembre", "décembre"],
        ["janv", "févr", "mars", "avril", "mai", "juin",
         "juil", "août", "sept", "oct", "nov", "déc"],
        ["lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi", "dimanche"],
        [""],
    )
    @test DatesPlus.dayname(Nov; locale="french") == "lundi"
    @test DatesPlus.dayname(Jan; locale="french") == "mardi"
    @test DatesPlus.dayname(Dec; locale="french") == "mercredi"
    @test DatesPlus.dayname(Apr; locale="french") == "jeudi"
    @test DatesPlus.dayname(Jun; locale="french") == "vendredi"
    @test DatesPlus.dayname(Feb; locale="french") == "samedi"
    @test DatesPlus.dayname(May; locale="french") == "dimanche"

    @test DatesPlus.monthname(Jan; locale="french") == "janvier"
    @test DatesPlus.monthname(Feb; locale="french") == "février"
    @test DatesPlus.monthname(Mar; locale="french") == "mars"
    @test DatesPlus.monthname(Apr; locale="french") == "avril"
    @test DatesPlus.monthname(May; locale="french") == "mai"
    @test DatesPlus.monthname(Jun; locale="french") == "juin"
    @test DatesPlus.monthname(Jul; locale="french") == "juillet"
    @test DatesPlus.monthname(Aug; locale="french") == "août"
    @test DatesPlus.monthname(Sep; locale="french") == "septembre"
    @test DatesPlus.monthname(Oct; locale="french") == "octobre"
    @test DatesPlus.monthname(Nov; locale="french") == "novembre"
    @test DatesPlus.monthname(Dec; locale="french") == "décembre"
end
@testset "leap years" begin
    @test DatesPlus.isleapyear(DatesPlus.DateTime(1900)) == false
    @test DatesPlus.isleapyear(DatesPlus.DateTime(2000)) == true
    @test DatesPlus.isleapyear(DatesPlus.DateTime(2004)) == true
    @test DatesPlus.isleapyear(DatesPlus.DateTime(2008)) == true
    @test DatesPlus.isleapyear(DatesPlus.DateTime(0)) == true
    @test DatesPlus.isleapyear(DatesPlus.DateTime(1)) == false
    @test DatesPlus.isleapyear(DatesPlus.DateTime(-1)) == false
    @test DatesPlus.isleapyear(DatesPlus.DateTime(4)) == true
    @test DatesPlus.isleapyear(DatesPlus.DateTime(-4)) == true

    @test DatesPlus.daysinyear(2000) == 366
    @test DatesPlus.daysinyear(2001) == 365
    @test DatesPlus.daysinyear(2000) == 366
    @test DatesPlus.daysinyear(2001) == 365

    @test DatesPlus.daysinyear(DatesPlus.Date(2000)) == 366
    @test DatesPlus.daysinyear(DatesPlus.Date(2001)) == 365
    @test DatesPlus.daysinyear(DatesPlus.DateTime(2000)) == 366
    @test DatesPlus.daysinyear(DatesPlus.DateTime(2001)) == 365
end
@testset "dayofweek" begin
    # Days of week from Monday = 1 to Sunday = 7
    @test DatesPlus.dayofweek(DatesPlus.DateTime(2013, 12, 22)) == 7
    @test DatesPlus.dayofweek(DatesPlus.DateTime(2013, 12, 23)) == 1
    @test DatesPlus.dayofweek(DatesPlus.DateTime(2013, 12, 24)) == 2
    @test DatesPlus.dayofweek(DatesPlus.DateTime(2013, 12, 25)) == 3
    @test DatesPlus.dayofweek(DatesPlus.DateTime(2013, 12, 26)) == 4
    @test DatesPlus.dayofweek(DatesPlus.DateTime(2013, 12, 27)) == 5
    @test DatesPlus.dayofweek(DatesPlus.DateTime(2013, 12, 28)) == 6
    @test DatesPlus.dayofweek(DatesPlus.DateTime(2013, 12, 29)) == 7
end
@testset "daysofweekinmonth" begin
    # There are 5 Sundays in December, 2013
    @test DatesPlus.daysofweekinmonth(DatesPlus.DateTime(2013, 12, 1)) == 5
    # There are 4 Sundays in November, 2013
    @test DatesPlus.daysofweekinmonth(DatesPlus.DateTime(2013, 11, 24)) == 4

    @test DatesPlus.dayofweekofmonth(DatesPlus.DateTime(2013, 12, 1)) == 1
    @test DatesPlus.dayofweekofmonth(DatesPlus.DateTime(2013, 12, 8)) == 2
    @test DatesPlus.dayofweekofmonth(DatesPlus.DateTime(2013, 12, 15)) == 3
    @test DatesPlus.dayofweekofmonth(DatesPlus.DateTime(2013, 12, 22)) == 4
    @test DatesPlus.dayofweekofmonth(DatesPlus.DateTime(2013, 12, 29)) == 5
end
@testset "dayofyear" begin
    @test DatesPlus.dayofyear(2000, 1, 1) == 1
    @test DatesPlus.dayofyear(2004, 1, 1) == 1
    @test DatesPlus.dayofyear(20013, 1, 1) == 1
    # Leap year
    @test DatesPlus.dayofyear(2000, 12, 31) == 366
    # Non-leap year
    @test DatesPlus.dayofyear(2001, 12, 31) == 365

    @test DatesPlus.dayofyear(DatesPlus.DateTime(2000, 1, 1)) == 1
    @test DatesPlus.dayofyear(DatesPlus.DateTime(2004, 1, 1)) == 1
    @test DatesPlus.dayofyear(DatesPlus.DateTime(20013, 1, 1)) == 1
    # Leap year
    @test DatesPlus.dayofyear(DatesPlus.DateTime(2000, 12, 31)) == 366
    # Non-leap year
    @test DatesPlus.dayofyear(DatesPlus.DateTime(2001, 12, 31)) == 365
    # Test every day of a year
    dt = DatesPlus.DateTime(2000, 1, 1)
    for i = 1:366
        @test DatesPlus.dayofyear(dt) == i
        dt += DatesPlus.Day(1)
    end
    dt = DatesPlus.DateTime(2001, 1, 1)
    for i = 1:365
        @test DatesPlus.dayofyear(dt) == i
        dt += DatesPlus.Day(1)
    end
end
@testset "quarterofyear" begin
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 1, 1))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 1, 31))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 2, 1))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 2, 29))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 3, 1))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 3, 31))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 4, 1)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 4, 30)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 5, 1)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 5, 31)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 6, 1)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 6, 30)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 7, 1)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 7, 31)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 8, 1)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 8, 31)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 9, 1)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 9, 30)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 10, 1)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 10, 31)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 11, 1)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 11, 30)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 12, 1)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.Date(2000, 12, 31)) == 4

    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 1, 1))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 1, 31))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 2, 1))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 2, 29))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 3, 1))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 3, 31))  == 1
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 4, 1)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 4, 30)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 5, 1)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 5, 31)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 6, 1)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 6, 30)) == 2
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 7, 1)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 7, 31)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 8, 1)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 8, 31)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 9, 1)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 9, 30)) == 3
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 10, 1)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 10, 31)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 11, 1)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 11, 30)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 12, 1)) == 4
    @test DatesPlus.quarterofyear(DatesPlus.DateTime(2000, 12, 31)) == 4
end
@testset "dayofquarter" begin
    # Non-leap and leap year
    for y in (2014, 2020)
        @test DatesPlus.dayofquarter(DatesPlus.Date(y, 1, 1)) == 1
        @test DatesPlus.dayofquarter(DatesPlus.Date(y, 4, 1)) == 1
        @test DatesPlus.dayofquarter(DatesPlus.Date(y, 7, 1)) == 1
        @test DatesPlus.dayofquarter(DatesPlus.Date(y, 10, 1)) == 1
        @test DatesPlus.dayofquarter(DatesPlus.Date(y, 3, 31)) == 90 + isleapyear(y)
        @test DatesPlus.dayofquarter(DatesPlus.Date(y, 6, 30)) == 91
        @test DatesPlus.dayofquarter(DatesPlus.Date(y, 9, 30)) == 92
        @test DatesPlus.dayofquarter(DatesPlus.Date(y, 12, 31)) == 92
        @test DatesPlus.dayofquarter(DatesPlus.DateTime(y, 1, 1)) == 1
        @test DatesPlus.dayofquarter(DatesPlus.DateTime(y, 4, 1)) == 1
        @test DatesPlus.dayofquarter(DatesPlus.DateTime(y, 7, 1)) == 1
        @test DatesPlus.dayofquarter(DatesPlus.DateTime(y, 10, 1)) == 1
        @test DatesPlus.dayofquarter(DatesPlus.DateTime(y, 3, 31)) == 90 + isleapyear(y)
        @test DatesPlus.dayofquarter(DatesPlus.DateTime(y, 6, 30)) == 91
        @test DatesPlus.dayofquarter(DatesPlus.DateTime(y, 9, 30)) == 92
        @test DatesPlus.dayofquarter(DatesPlus.DateTime(y, 12, 31)) == 92
    end
end

end

ext: jl


text:
"""
`solve(P::Problem, O::Options=DefaultOptions)`

Solve the problem `P` using options `O`.

See `Problem` and `Options` for more info about how to state a problem and change the default options.
"""
function solve(P::Problem, O::Options=DefaultOptions)::Vector

    # Unpack the data structures
    @unpack p, q, f, a, α, A, b, β, B = P
    @unpack n = O

    # Compute the step
    h = (b - a) / n

    # Compute auxiliary variables
    h⁻¹ = 1 / h
    hₕ = h / 2
    ax₀ = h⁻¹ * (hₕ * p(a) + 1)
    axₙ = h⁻¹ * (hₕ * p(b) - 1)

    # Compute diagonals and right-hand column
    dl = OffsetArray([[1 - hₕ * p(a + k * h) for k in 1:n-1]; axₙ], 2:n+1)
    d = [-ax₀ + hₕ * q(a) - α; [h^2 * q(a + k * h) - 2 for k in 1:n-1]; -axₙ - hₕ * q(b) - β]
    du = [ax₀; [1 + hₕ * p(a + k * h) for k in 1:n-1]]
    f = [hₕ * f(a) + A; [h^2 * f(a + k * h) for k in 1:n-1]; B - hₕ * f(b)]

    # Perform the forward sweep
    for i in 2:n+1
        w = dl[i] / d[i-1]
        d[i] = d[i] - w * du[i-1]
        f[i] = f[i] - w * f[i-1]
    end

    # Prepare the result array
    y = Vector{Float64}(undef, n+1)

    # Perform the back substitution
    y[n+1] = f[n+1] / d[n+1]
    for i in n:-1:1
        y[i] = (f[i] - du[i] * y[i+1]) / d[i]
    end

    return y

end


text:
using Images
#Additional required packages to run this example.
#Pkg.add("Images")

using AprilTags
using Dates

# Simple method to show the image with the tags
function showImage(image, tags)
    # Convert image to RGB
    imageCol = RGB.(image)
    #traw color box on tag corners
    foreach(tag->drawTagBox!(imageCol, tag), tags)
    imageCol
end

##
# Create default detector
detector = AprilTagDetector()

#Run against a file
image = load(dirname(Base.source_path()) *"/../data/tagtest.jpg")

@sync begin
    starttime = Dates.value(now())
    #call detector in thread
    @async global t1 = @timed begin
        println("time before detector $(Dates.value(now())-starttime) ms")
        global tags = threadcalldetect(detector, image)
        # ↑ comment --- compare with this --- uncomment ↓
        # global tags = detector(image)
        println("time after detector $(Dates.value(now())-starttime) ms")

        @show length(tags)
    end
    #carry on doing something else
    @async global t2 = @timed begin
        println("time starting other $(Dates.value(now())-starttime) ms")
        imageCol = load(dirname(Base.source_path()) *"/../data/colortag.jpg")
        A = randn(2000,2000)
        randsum = sum(A*A')
        println("time finished other $(Dates.value(now())-starttime) ms")
        @show randsum
    end
end

    println("time for detector: $(t1[2]) seconds")
    println("time for other: $(t2[2]) seconds")
    println()

showImage(image, tags)

## clean up detector once finished
freeDetector!(detector)


text:
using Test

@testset "Efficiency" begin
    sdd = SDDEfficiency(AP33Model())
    sili = SiLiEfficiency(Beryllium(10.0e-4))

    @test isapprox(efficiency(sdd, 300.0), 0.152, atol = 0.001)
    @test isapprox(efficiency(sdd, 3000.0), 0.753, atol = 0.001)
    @test isapprox(efficiency(sdd, 30000.0), 0.091, atol = 0.001)

    @test isapprox(efficiency(sili, 300.0), 0.000, atol = 0.001)
    @test isapprox(efficiency(sili, 3000.0), 0.957, atol = 0.001)
    @test isapprox(efficiency(sili, 30000.0), 0.485, atol = 0.001)
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
importall POMDPs
using Distributions

type D4ICRA <: POMDPs.POMDP{Float64,Int,Int}
	discount_factor::Float64
	step_size::Int64
end

#2second = 50000 or maybe 1 second?

function makeSolver()

	actions(::D4ICRA) = [0,1,2,3]
	n_actions(::D4ICRA) = 4
	pomdp = D4ICRA(0.9,1);
	solver = POMCPSolver(tree_queries=50000, c=10)
	planner = solve(solver, pomdp);

	return [solver,pomdp,planner];

end

#observations, function of state only
function generate_o(p::D4ICRA, s::Array, a::Int, sp::Array, rng::AbstractRNG)

	#ROBOT_VIEW_RADIUS
	if sqrt((s[1]-s[3])*(s[1]-s[3]) + (s[2]-s[4])*(s[2]-s[4])) < 20
		return 1
	else
		return 0
	end

end

#state transitions, function of state and action
function generate_sr(p::D4ICRA, s::Array, a::Int, rng::AbstractRNG)

	delta = 10;
    if a == 0
        s = s + [0,delta,0,0]
    elseif a == 1
        s= s + [0,-delta,0,0]
    elseif a == 2
    	s= s + [delta,0,0,0]
    elseif a == 3
    	s= s + [-delta,0,0,0]
    else
    	s= s
    end

    s[1] = max(0,s[1])
    s[1] = min(437,s[1]);
    s[2] = max(0,s[2])
    s[2] = min(754,s[2]);

    r = 0;

    if sqrt((s[1]-s[3])*(s[1]-s[3]) + (s[2]-s[4])*(s[2]-s[4])) < 25
		r = 10
	else
		r = 0
	end

    return s,r

end

function initial_state_distribution(pomdp::D4ICRA)
    sig_0 = [10,10,sqrt(0.00001),sqrt(0.00001)]
    return MvNormal(sig_0)
end;

using BasicPOMCP
using POMDPToolbox
importall Base
using Dates

function getPOMCPAction(planner,b)
	return action(planner,b);
end

actions(::D4ICRA) = [0,1,2,3]
n_actions(::D4ICRA) = 4
#solver,pomdp,planner = makeSolver();
#b = initial_state_distribution(pomdp);
#act = getPOMCPAction(planner,b);
#println(act);


text:
module test_pkgs_dataflow

using Base.Test
using DataFlow # IVertex DataFlow.il
               # DVertex DataFlow.dl
using DataFlow: prewalk, walk

# https://github.com/MikeInnes/DataFlow.jl/blob/master/docs/vertices.md

# input-linked
iv = IVertex(2)
@test Vector{IVertex{Int}}() == iv.inputs
@test 2 == iv.value

@test prewalk(identity, iv) == iv
@test walk(iv, identity, identity) == iv

@test contains(iv, 2)

@test DataFlow.il(iv) == iv

# doubly-linked
dv = DVertex(2)
@test Vector{DVertex{Int}}() == dv.inputs
@test 2 == dv.value
@test DataFlow.dl(iv) == dv
@test DataFlow.il(dv) == iv

# Syntax → Graph
import DataFlow: syntax, graphm

ex = syntax(dv)
@test ex isa Expr
@test ex.args[1].args[1] == 2
@test graphm(Dict(), ex) == dv

end # module test_pkgs_dataflow


text:
struct CameraFps
    position::Vector{Float64}
    up::Vector{Float64}
    along::Vector{Float64}
    forward::Vector{Float64}
end

function CameraFps()
    return CameraFps([0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, -1.0])
end

function view_matrix(camera)
    x = dot(camera.along, camera.position)
    y = dot(camera.up, camera.position)
    z = dot(camera.forward, camera.position)

    return [
        camera.along  camera.up -camera.forward zeros(3);
       -x            -y          z              1.0
    ]'
end

function walk!(camera, delta)
    camera.position .-= delta .* camera.forward
    return camera
end

function strafe!(camera, delta)
    camera.position .+= delta .* camera.along
    return camera
end

function yaw!(camera, angle)
    @. camera.forward = cos(angle) * camera.forward + sin(angle) * camera.along

    normalize!(camera.forward)

    camera.along .= cross(camera.forward, camera.up)

    return camera
end


text:
# paste code in the REPL:
using WebSockets, JSON

url = "wss://www.bitmex.com/realtime"

payload = Dict(                         # 1
               :op => "subscribe",
               :args => "trade:XBT"
           )

function open_websocket()
   WebSockets.open(url) do ws
     if isopen(ws)
       write(ws, JSON.json(payload))    # 2
     end

     while isopen(ws)
       data, success = readguarded(ws)
       !success && break
       data = JSON.parse(String(data))
       print(data, "\n")
     end

     if !isopen(ws)
       @async open_websocket()
     end

   end
 end

@async open_websocket()


text:
import J4Org: extract_links

@testset "extract_links" begin
    s="# some text without link [c]  [[a][b] "
    a=extract_links(s)
    @test isempty(a)

    s="# some text [[some_target][]]"
    a=extract_links(s)
    @test length(a)==1
    @test a[1]=="some_target"

    s="# some text [[some_target][]] another one [[link_target][link_name]] and a last one [[a4][]]"
    a=extract_links(s)
    @test length(a)==2
    @test a[1]=="some_target"
    @test a[2]=="a4"

    s="# some text [[some_target][]] another one \"[[link_target][]]\" and a last one [[a4][]]"
    a=extract_links(s)
    @test length(a)==2
    @test a[1]=="some_target"
    @test a[2]=="a4"
end;


text:
@testset "Sensors" begin

    fname = joinpath(dirname(@__FILE__), "..", "..", "data", "test.sfp")
    s = read_sfp(fname)

    @testset "Read file" begin

        s = read_sfp(fname)
        @test length(s) == 3
    end

    @testset "Show" begin

        show(s)
        show(s[1])
    end

    @testset "Info" begin

        @test label(s) == ["Fp1", "Fpz", "Fp2"]
        @test label(s[1]) == "Fp1"
        @test labels(s) == ["Fp1", "Fpz", "Fp2"]
        @test labels(s[1]) == "Fp1"

        @test EEG.x(s) == [-27.747648, -0.085967, 27.676888]
        @test EEG.x(s[1]) == -27.747648

        @test EEG.y(s) == [98.803864, 103.555275, 99.133354]
        @test EEG.y(s[1]) == 98.803864

        @test EEG.z(s) == [34.338360, 34.357265, 34.457005]
        @test EEG.z(s[1]) == 34.338360
    end

    @testset "Matching" begin

        a, b = match_sensors(s, ["Fp1", "Fp2"])

    end

    @testset "Leadfield matching" begin

        L = ones(4, 3, 5)
        L[:, :, 2] = L[:, :, 2] * 2
        L[:, :, 3] = L[:, :, 3] * 3
        L[:, :, 4] = L[:, :, 4] * 4
        L[:, :, 5] = L[:, :, 5] * 5

        L, idx = match_sensors(L, ["meh", "Cz", "Fp1", "Fp2", "eh"], ["Fp2", "Cz"])

        @test size(L) == (4, 3, 2)
        @test L[:, :, 1] == ones(4, 3) * 4
        @test L[:, :, 2] == ones(4, 3) * 2
        @test idx == [4, 2]

    end

    @testset "Standard sets" begin

        @test length(EEG_64_10_20) == 64
        @test length(EEG_Vanvooren_2014) == 18
        @test length(EEG_Vanvooren_2014_Left) == 9
        @test length(EEG_Vanvooren_2014_Right) == 9

    end

end


text:
module FluidIntonation

import MIDI
import PyCall
import Match

import FunctionalCollections as FC
import Rocket as R

include("stream.jl")
include("loss.jl")
include("midi_io.jl")

TICK = 1

NOTE_ON_SIGNAL= 144
NOTE_OFF_SIGNAL= 128
SUSTAIN_ON_SIGNAL = 176

channel(pitch) = pitch % 12

centtoratio(cent) = 10^((log10(2) / 1200) * cent)
pitchbend(bend_ratio) = round.(8192 .+ 4096 * 12 * log2.(bend_ratio))

function possible_adjustments(event, state)
    new_state = Match.@match event begin
        MIDI.NoteOnEvent(_, _, note, _)  => haskey(state, note) ? state : FC.assoc(state, note, 0)
        MIDI.NoteOffEvent(_, _, note, _) => haskey(state, note) ? FC.dissoc(state, note) : state
        _                                => state
    end

    if (length(new_state) > 1)
        note_adjustment_pairs = [(a, b) for (a, b) in new_state]
        new_adjustments = adjust(note_adjustment_pairs)
        new_adjustment_pairs = collect(zip(map(first, note_adjustment_pairs), new_adjustments))
        bends = map(new_adjustment_pairs) do (note, adj)
            MIDI.PitchBendEvent(TICK, Int(round(pitchbend(centtoratio(adj)))), channel=channel(note))
        end
        return bends, FC.PersistentHashMap(new_adjustment_pairs)
    else
        return MIDI.MIDIEvent[], new_state
    end
end

function make_input_event(msg)
    (signal, pitch, velocity), _ = msg
    if signal == NOTE_ON_SIGNAL
        MIDI.NoteOnEvent(TICK, pitch, velocity, channel=channel(pitch))
    elseif signal == NOTE_OFF_SIGNAL
        MIDI.NoteOffEvent(TICK, pitch, velocity, channel=channel(pitch))
    else
        error("Could not recognize signal")
    end
end

function make_output_event(evt)
    out = Match.@match evt begin
        MIDI.NoteOnEvent(_, status, note, velocity)  => [status, note, velocity]
        MIDI.NoteOffEvent(_, status, note, velocity) => [status, note, velocity]
        MIDI.PitchBendEvent(_, status, pitch)        => [status, midibits(UInt16(pitch))...]
        other                                        => error(other)
    end
    PyCall.PyVector(out)
end

function midibits(value::UInt16)
    # split into two 7 bit segments
    mask = UInt16((1 << 7) - 1)
    lsb = value & mask
    msb = (value & (mask << 7)) >> 7
    (lsb, msb)
end

source = R.Subject(Tuple{Vector{Int64}, Float64}, scheduler = R.AsyncScheduler())

input_stream = source |>
    R.map(MIDI.MIDIEvent, make_input_event) |>
    R.take_until(source |> filter(x -> typeof(x) == MIDI.ControlChangeEvent))

intervals = R.interval(100)

adjustment_stream = R.merged((input_stream, intervals)) |>
    stateful_flatmap(MIDI.MIDIEvent, possible_adjustments, FC.PersistentHashMap{Int, Float64}())

output_stream = R.merged((input_stream, adjustment_stream)) |>
    R.map(PyCall.PyVector, make_output_event)

function test()
    println("Starting test")
    R.subscribe!(output_stream, (msg) -> @show msg)
    R.next!(source, ([NOTE_ON_SIGNAL, 100, 100], 1.0))
    R.next!(source, ([NOTE_ON_SIGNAL, 120, 100], 1.0))
    R.next!(source, ([NOTE_ON_SIGNAL, 123, 80], 1.0))
    R.next!(source, ([NOTE_OFF_SIGNAL, 100, 100], 1.0))
    println("Done")
end

function start()
    rtmidi = PyCall.pyimport("rtmidi")
    with_midi(port=8, driver_class=rtmidi.MidiOut) do midi_out
        R.subscribe!(output_stream, (msg) -> midi_out.send_message(msg))

        with_midi(transmit_midi_input_to(source), port=0, driver_class=rtmidi.MidiIn)
    end
end

end # module

max_stars_repo_name: miksu/FluidIntonation

max_issues_repo_name: miksu/FluidIntonation

max_forks_repo_name: miksu/FluidIntonation


text:
"""Definition and functions working with the Lattice type. Lattice is a lightweight
    wrapper that stores lattice geometry in absolute coordinates for fast site indexing,
    along with a box size for easy collection of all sites.
   Lattice is arbitrary-dimensional, but contains no symmetry information.
"""

"""
    Lattice

A type holding geometry information about a 3D lattice in a simulation box.
"""
struct Lattice <: AbstractArray{SVector{3, Float64}, 4}
    lat_vecs      :: SMatrix{3, 3, Float64, 9}       # Columns are lattice vectors
    basis_vecs    :: Vector{SVector{3, Float64}}     # Each SVector gives a basis vector
    types         :: Vector{String}                  # Indices labeling atom types
    size          :: SVector{3, Int}                 # Number of cells along each dimension

    @doc """
        Lattice(lat_vecs, basis_vecs, types, latsize)

    Construct a 3-dimensional `Lattice`.
    # Arguments
    - `lat_vecs`: A matrix where the lattice vectors form the columns.
                  Must be `convert`-able into `SMatrix{3, 3, Float64, 9}`.
    - `basis_vecs`: A `Vector` of basis positions of sites within the unit cell,
                     given in fractional coordinates.
                    Each element must be `convert`-able into `SVector{3, Float64}`.
    - `types::Vector{String}`: A list of atomic types identifiers for each site.
                                 Equivalent sites should have the same identifier.
    - `latsize`: Specifies the number of unit cells extending along each lattice vector.
                 Must be `convert`-able into `SVector{3, Int}`.
    """
    function Lattice(lat_vecs, basis_vecs, types, latsize)
        @assert all(isequal(3), size(lat_vecs))          "All `lat_vecs` dims must be length 3"
        @assert all(isequal(3), map(length, basis_vecs)) "All basis_vecs should be length 3 to match lat_vecs"
        @assert length(basis_vecs) > 0                   "At least one basis atom required"
        @assert length(basis_vecs) == length(types)      "Length of basis_vecs and types should match"
        @assert all(v->all(0 .<= v .< 1), basis_vecs)    "All basis_vecs should be given in fractional coordinates [0, 1)"
        @assert length(latsize) == 3                     "latsize should be length 3"
        lat_vecs = convert(SMatrix{3, 3, Float64, 9}, lat_vecs)
        basis_vecs = [lat_vecs * convert(SVector{3, Float64}, b) for b in basis_vecs]
        latsize = SVector{3, Int}(latsize)
        new(lat_vecs, basis_vecs, types, latsize)
    end
end

"""
    Lattice(lat_vecs, basis_vecs, latsize)

Construct a `Lattice`, with the dimension inferred by the shape of `lat_vecs`,
and all sites assumed to be the same type (labeled `"A"`).
"""
Lattice(lvecs, bvecs, latsize) = Lattice(lvecs, bvecs, fill("A", length(bvecs)), latsize)

"""
    lattice_params(lattice::Lattice)
"""

function Base.show(io::IO, ::MIME"text/plain", lattice::Lattice)
    println(io, join(size(lattice), "x"), " Lattice")
    # Print out lattice vectors, types, basis?
end

# Constructors converting Lattice -> Crystal, and Crystal -> Lattice

"""
    Crystal(lattice::Lattice)

Construct a `Crystal` using geometry information in `lattice`, inferring symmetry information.
"""
function Crystal(lattice::Lattice)
    L = lattice.lat_vecs
    # Convert absolute basis positions to fractional coordinates
    basis_coords = map(b -> inv(L) * b, lattice.basis_vecs)
    Crystal(L, basis_coords, types=lattice.types)
end

function Lattice(cryst::Crystal, latsize) :: Lattice
    Lattice(cryst.lat_vecs, cryst.positions, cryst.types, latsize)
end

# Number of sublattices.
nbasis(lat::Lattice) = length(lat.basis_vecs)
# Volume of a unit cell.
cell_volume(lat::Lattice) = abs(det(lat.lat_vecs))
# Volume of the full simulation box.
volume(lat::Lattice) = cell_volume(lat) * prod(lat.size)
lattice_vectors(lat::Lattice) = lat.lat_vecs
lattice_params(lat::Lattice) = lattice_params(lat.lat_vecs)

"Produce an iterator over all unit cell indices."
@inline function eachcellindex(lat::Lattice) :: CartesianIndices{3}
    return CartesianIndices(Tuple(lat.size))
end

@inline function Base.size(lat::Lattice)
    return (nbasis(lat), lat.size...)
end

#=== Indexing returns absolute coordinates of lattice points ===#

@inline function Base.getindex(lat::Lattice, b::Int64, brav::CartesianIndex{3})
    return lat.lat_vecs * convert(SVector{3, Float64}, brav) + lat.basis_vecs[b]
end

@inline function Base.getindex(lat::Lattice, b::Int64, brav::NTuple{3, Int64})
    return lat.lat_vecs * convert(SVector{3, Float64}, brav) + lat.basis_vecs[b]
end

@inline function Base.getindex(lat::Lattice, b::Int64, brav::Vararg{Int64, 3})
    return lat.lat_vecs * convert(SVector{3, Float64}, brav) + lat.basis_vecs[b]
end

# TODO: Should just be another Lattice
# TODO: Potentially confusing. `gen_reciprocal` computes the reciprocal lattice
#        vectors of the unit cell, not the entire simulation box. Most use cases
#        want the former, but Ewald summations want the latter.
#       These differ just by a scale factor, which is currently explicitly handled
#        in Ewald summation calculations.
"Defines a reciprocal lattice structure"
struct ReciprocalLattice
    lat_vecs  :: SMatrix{3, 3, Float64, 9}     # Columns of this are the reciprocal lattice vectors
    size      :: SVector{3, Int}
end

function ReciprocalLattice(lat::Lattice)
    recip_vecs = 2π * transpose(inv(lat.lat_vecs))
end

function Base.eachindex(lat::ReciprocalLattice) :: CartesianIndices
    return CartesianIndices(Tuple(lat.size))
end

"Generate a `ReciprocalLattice` for a given `Lattice`"
function gen_reciprocal(lat::Lattice) :: ReciprocalLattice
    recip_vecs = 2π * transpose(inv(lat.lat_vecs))
    return ReciprocalLattice(recip_vecs, lat.size)
end

"Returns just the underlying Bravais lattice"
function brav_lattice(lat::Lattice) :: Lattice
    return Lattice(
        lat.lat_vecs,
        [@SVector zeros(3)],
        ["A"],
        lat.size
    )
end

cell_type(lattice::Lattice) = cell_type(lattice.lat_vecs)

function distance(lat::Lattice, b::Bond)
    return norm(lat.lat_vecs * b.n + (lat.basis_vecs[b.j] - lat.basis_vecs[b.i]))
end


text:
using HTTP
using JSON

function anki_add_tag(card_id, tag)
	if length(card_id) == 1
		card_id = "[$card_id]"
	end

	response = HTTP.request(:POST, "http://localhost:8765", [], """{
				\"action\": \"addTags\",
				\"version\": 6,
				\"params\": {
					\"notes\": $card_id,
					\"tags\": \"$tag\"
					}
				}""")
	JSON.parse(String(response.body))
end

function anki_remove_tag(card_id, tag)
	if length(card_id) == 1
		card_id = "[$card_id]"
	end

	response = HTTP.request(:POST, "http://localhost:8765", [], """{
				\"action\": \"removeTags\",
				\"version\": 6,
				\"params\": {
					\"notes\": $card_id,
					\"tags\": \"$tag\"
					}
				}""")
	JSON.parse(String(response.body))
end


text:
import JSON
import Nemo: gen, terms, monomials

include("myeval.jl")

#------------------------------------------------------------------------------

# returns a model parsed from the given `filename`
# (with matrices in dense format)
# each returned matrix is represented with a dense array
function load_dense_set(filename::String)
    ans = []

    filepath = "$(normpath(joinpath(@__FILE__, "..", "..")))/src/data/matrices/$(filename)"

    open(filepath, "r") do input_stream
        for arr in JSON.Parser.parse(input_stream)
            push!(ans, arr)
        end
    end

    return ans
end

#------------------------------------------------------------------------------

# returns a model parsed from the given `filename`
# (with matrices in coordinate format)
# each returned matrix is represented with an array of
# items (i, j, x) where typeof(x) is Rational
function load_COO_set(filename::String)
    model = []

    filepath = "$(normpath(joinpath(@__FILE__, "..", "..")))/src/data/matrices/$(filename)"

    str_to_nom = (x -> parse(BigInt, first(split(x, '/'))))
    str_to_den = (x -> (if occursin('/', x) parse(BigInt, last(split(x, '/'))) else 1 end))

    open(filepath, "r") do input_stream
        for arr in JSON.Parser.parse(input_stream)
            normal_repr = [
                arr[1], arr[2],
                map(x -> (x[1], x[2], str_to_nom(x[3]) // str_to_den(x[3])), arr[3])
            ]

            push!(model, normal_repr)
        end
    end

    return model
end

#------------------------------------------------------------------------------

# returns an array of all models from data/matrices that satisfy
# the given criterias:
#   `from_dim`  - lower bound for the model's matrices one-side dimension
#   `to_dim`    - upper bound -//-
#   `from_size` - lower bound for the model's total number of matrices
#   `to_size`   - upper bound -//-
#
# each returned model is represented with an array:
#   [filename, dim, size, matrices-data]
# each returned matrix is represented with an array of
# items (i, j, x) of nonzeroes where typeof(x) is Rational
function load_COO_if(;from_dim, to_dim, from_size=0, to_size=Inf)
    models = []

    for filename in readdir("src/data/matrices/")
        model_data = load_COO_set(filename)
        model_dim = first(first(model_data))
        model_size = length(model_data)

        if from_dim <= model_dim <= to_dim &&
                from_size <= model_size <= to_size
            push!(models, [filename, model_dim, model_size, model_data])
        end
    end

    @info "loaded models: $(length(models))"

    return models
end

#------------------------------------------------------------------------------

# returns a model parsed from the given `group\name` file
# (with matrices in MTX format)
# each returned matrix is represented with an array of
#   [filename, dim, size, data]
# where data is stored as items (i, j, x) where typeof(x) is Rational
#
# In case of any parse error occurred, e.g inexact rationals encountered,
# returs an empty array []
function load_MTX(group, name)
    model = []

    # seems to be a default separator
    sep = r"%-+\r?\n"

    filepath = replace(
        "$(normpath(joinpath(@__FILE__, "..", "..")))src/data/sscollection/$group/$name",
        "\\" => "/"
    )

    open(filepath, "r") do input_stream
        # some models have >3 sections, taking the last
        data = last(split(read(input_stream, String), sep))

        # should be ok
        data = split(data, r"\r?\n")
        meta = map(x -> parse(Int, x), split(first(data), ' '))
        # strange matrices, discard them
        # NEED TO BE HANDLED IN ANOTHER WAY
        if length(meta) != 3
            return []
        end

        m, n, nnz = meta
        # we want only square ones
        if m != n
            return []
        end

        try
            model = [
                group * '/' * name, # filename
                m,                  # dim
                1,                  # size
                collect(            # data
                    map(
                        # `Any` disallows auto-promoting Int to Rational
                        X -> Any[
                            parse(Int, X[1]),
                            parse(Int, X[2]),
                            length(X) > 2 ? Rational(parse(Float64, X[3])) : Rational(1)
                        ] ,
                        filter(!isempty, map(split, view(data, 1:length(data))))
                    )
                )
            ]
        catch e
            if isa(e, InexactError)
                # occurrs if floats are impossibly small/big,
                # discard the model then
                return []
            else
                rethrow(e)
            end
        end
    end

    return model
end

#------------------------------------------------------------------------------

# returns an array of all models from data/sscollection that satisfy
# the given criterias:
#   `from_dim`  - lower bound for the model's matrices one-side dimension
#   `to_dim`    - upper bound -//-
#
# each returned model is represented with an array:
#   [filename, dim, size, matrices-data]
# each returned matrix is represented with an array of
# items (i, j, x) of nonzeroes where typeof(x) is Rational
function load_MTX_if(;from_dim=0, to_dim=Inf)
    models = []

    extension = ".mtx"
    dir = "src/data/sscollection/"

    for group in readdir(dir)
        for name in readdir("$dir/$group/")
            if endswith(name, extension)
                model = load_MTX(group, name)
                if ! isempty(model)
                    fn, dim, sz, data = model
                    if from_dim <= model[2] <= to_dim
                        push!(models, [ fn, dim, sz, [ dim, dim, data ] ] )
                    end
                end
            end
        end
    end

    @info "loaded models: $(length(models))"

    return models
end

#------------------------------------------------------------------------------

function load_ODEs(filename; pathabs=false)
    # *.ode --> polynomial form of ode system

    filepath = filename
    if !pathabs
        filepath = replace(
            "$(normpath(joinpath(@__FILE__, "..", "..")))src/data/ODEs/$filename",
            "\\" => "/"
        )
    end

    lines = []
    open(filepath, "r") do inputs
        lines = map(strip, readlines(inputs))
    end

    # drop init out
    reactions = filter(!isempty, map(strip, lines[
        (findfirst(startswith("begin reactions"), lines) + 1
        :
        findlast(startswith("end reactions"), lines) - 1)
    ]))

    # scan reactions to discover all variables
    sep1 = r"(->)|(,)|([+-/ \t\*])"
    # assuming there is at least one utf letter in any variable
    strings = map(String, unique(filter!(
        x -> !isempty(findall(isletter, x)),
        split(join(reactions, ' '), sep1)
    )))

    S, xs = QQ[strings...]

    # symbol :x to x from QQ[x]
    mapping = Dict{Symbol, fmpq_mpoly}(
        Meta.parse(x) => gen(S, i)
        for (i, x) in enumerate(strings)
    )

    ODEs = Dict(x => S(0) for x in xs)

    sep2 = r"(->)|(,)"
    for reaction in reactions
        # divide each line into parts,
        # expected to be three of them
        lhs, rhs, speed = map(
            Meta.parse ∘ String ∘ strip,
            split(reaction, sep2)
        )

        reagents, products, speed = map(
            x -> myeval(x, mapping),
            (lhs, rhs, speed)
        )

        concentration = prod(terms(reagents))

        for reagent in monomials(reagents)
            ODEs[reagent] -= concentration * speed
        end
        for product in monomials(products)
            ODEs[product] += concentration * speed
        end

    end

    for key in collect(keys(ODEs))
        if iszero(ODEs[key])
            delete!(ODEs, key)
        end
    end

    @info "loaded a system of $(length(ODEs)) ODEs from $filename"

    return ODEs
end


text:
#=
持橋，大羽 著　「ガウス過程と機械学習」（機械学習プロフェッショナルシリーズ）
多変量ガウス分布のサンプリングのJuliaによる実装
=#

using LinearAlgebra
using Plots
using LaTeXStrings

μ1 = μ2 = 0.0

𝛍 = [
    μ1
    μ2
]

σ11 = σ22 = 1.0
σ12 = σ21 = -0.7

Σ = [
    σ11 σ12
    σ21 σ22
]

function 𝒩(𝐱...; 𝛍 = 𝛍, Σ = Σ)
    Λ = inv(Σ)
    D = length(𝛍)
    return √(det(Λ)) * exp(-0.5 * (𝐱 .- 𝛍)' * Λ * (𝐱 .- 𝛍)) / √(2π)^D
end

x1min, x1max = -4.0, 4.0
x2min, x2max = -4.0, 4.0
x1 = range(x1min, x1max, length = 50)
x2 = range(x2min, x2max, length = 50)
z = 𝒩.(x1', x2)
contour(x1, x2, z)
plot!(
    xlabel = L"x_1",
    ylabel = L"x_2",
    aspect_ratio = 1,
    xlims = [-4.0, 4.0],
    ylims = [-4.0, 4.0],
)
strΣ="\\sigma_{11}=$σ11, \\sigma_{12}=$σ12, \\sigma_{21}=$σ21, \\sigma_{22}=$σ22"
str𝛍="\\mu = [$μ1, $μ2]"
title!(latexstring("\\mathcal{N}(x \\mid \\mu, \\Sigma),\\ $str𝛍 , \\atop $strΣ"))

# Box–Muller transform を用いた方法
L = cholesky(Σ).L
cL = size(L, 2)

boxmuller(r1, r2) = √(-2 * log(r1)) * sin(2π * r2)
univariate(n...) = boxmuller.(rand(n...), rand(n...))

function multivariate(n)
    x = univariate(cL, n) # cL 次元ベクトルを n 個生成する
    return L * x # cL 次元ベクトル を変換して 𝒩(0, Σ) に従うサンプルを生成
end

xy = multivariate(100)
scatter!(xy[1, :], xy[2, :],color=:blue, label="boxmuller")

# Distributions.jl でサンプルを作成
d = MvNormal(𝛍, Σ)
xy = rand(d,100)
scatter!(xy[1,:], xy[2,:],color=:red, label="mvnormal")

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function cc_mod(x::Neuro1DAnalyticSignal)::Neuro1DRealSignal
    Neuro1DRealSignal(
        signal = abs.(x.signal),
        sample_rate = x.sample_rate
    )
end

function cc_arg(x::Neuro1DAnalyticSignal)::Neuro1DRealSignal
    Neuro1DRealSignal(
        signal = angle.(x.signal),
        sample_rate = x.sample_rate
    )
end

function cc_rel(x::Neuro1DAnalyticSignal)::Neuro1DRealSignal
    Neuro1DRealSignal(
        signal = real.(x.signal),
        sample_rate = x.sample_rate
    )
end


text:
include("../../src/Transform/dct.jl")

@testset "Transform/dct" begin
	a=rand(10, 20, 30)
	b=copy(a)
	dct!(b)
	idct!(b)
	@test isapprox(a, b)
	a=rand(12, 22, 32)
	b=copy(a)
	dct!(b)
	idct!(b)
	@test isapprox(a, b)
end


text:
@testset "pointer" begin

# constructors
voidptr_a = CuPtr{Cvoid}(Int(0xDEADBEEF))
@test reinterpret(Ptr{Cvoid}, voidptr_a) == Ptr{Cvoid}(Int(0xDEADBEEF))

# getters
@test eltype(voidptr_a) == Cvoid

# comparisons
voidptr_b = CuPtr{Cvoid}(Int(0xCAFEBABE))
@test voidptr_a != voidptr_b

@testset "conversions" begin

# between regular and CUDA pointers
@test_throws ArgumentError convert(Ptr{Cvoid}, voidptr_a)

# between CUDA pointers
intptr_a = CuPtr{Int}(Int(0xDEADBEEF))
@test convert(typeof(intptr_a), voidptr_a) == intptr_a

# convert back and forth from UInt
intptr_b = CuPtr{Int}(Int(0xDEADBEEF))
@test convert(UInt, intptr_b) == 0xDEADBEEF
@test convert(CuPtr{Int}, Int(0xDEADBEEF)) == intptr_b
@test Int(intptr_b) == Int(0xDEADBEEF)

# pointer arithmetic
intptr_c = CuPtr{Int}(Int(0xDEADBEEF))
intptr_d = 2 + intptr_c
@test isless(intptr_c, intptr_d)
@test intptr_d - intptr_c == 2
@test intptr_d - 2 == intptr_c
end

@testset "GPU or CPU integration" begin

a = [1]
ccall(:clock, Nothing, (Ptr{Int},), a)
@test_throws Exception ccall(:clock, Nothing, (CuPtr{Int},), a)
ccall(:clock, Nothing, (CUDAdrv.PtrOrCuPtr{Int},), a)

ptr = convert(CuPtr{eltype(a)}, CU_NULL)
b = CuTestArray{eltype(a), ndims(a)}(ptr, size(a))
ccall(:clock, Nothing, (CuPtr{Int},), b)
@test_throws Exception ccall(:clock, Nothing, (Ptr{Int},), b)
ccall(:clock, Nothing, (CUDAdrv.PtrOrCuPtr{Int},), b)

end

end


text:
module PreProcess
using Distributions

function isBadRow(x::Vector{T}; thresh::Float64=-6.0) where {T <: Number}
  return any(x .< thresh)
end

function removeBadRows!(y::Vector{Matrix{T}};
                        thresh::Float64=-6.0) where {T <: Number}
  I = length(y)
  N = size.(y, 1)

  goodRows = [[!isBadRow(y[i][n, :], thresh=thresh)
               for n in 1:N[i]] for i in 1:I]

  for i in 1:I
    y[i] = y[i][goodRows[i], :]
  end
end

#= Idea:
  - if across all samples, a marker has > 90% missing or negative, then remove
  - if across all samples, a marker has > 90% positive, then remove
=#

function isBadColumn(x::Vector{T};
                     maxNanOrNegProp::Float64=.9,
                     maxPosProp::Float64=.9)::Bool where {T <: Number}
  n = length(x)

  missOrNegProp = sum((isnan.(x)) .| (x .< 0)) / n
  posProp = sum(x .> 0) / n

  return missOrNegProp > maxNanOrNegProp || posProp > maxPosProp
end

function isBadColumnForAllSamples(j::Int, y::Vector{Matrix{T}};
                                  maxNanOrNegProp::Float64=.9,
                                  maxPosProp::Float64=.9) where {T}
  I = length(y)
  results = [isBadColumn(y[i][:, j],
                         maxNanOrNegProp=maxNanOrNegProp,
                         maxPosProp=maxPosProp) for i in 1:I]
  result = all(results)
  return result
end

function markerIsNearZero(j::Int, yi::Matrix{T};
                          thresh=0.5, prop=.25) where {T}
  out = (isnan.(yi[:, j]) .== false) .& (abs.(yi[:, j]) .< thresh)
  return mean(out) > prop
end

function markerIsNearZeroInAnySample(j::Int, y::Vector{Matrix{T}};
                                     thresh=0.5, prop=.25) where T
  I = length(y)

  badMarker = [markerIsNearZero(j, y[i], thresh=thresh, prop=prop)
               for i in 1:I]
  return any(badMarker)
end

"""
subsample data if 0 < `subsample` < 1
returns goodColumns, new_J
"""
function preprocess!(y::Vector{Matrix{T}}; maxNanOrNegProp::Float64=.9,
                     maxPosProp::Float64=.9, subsample::Float64=0.0,
                     rowThresh::Float64=-Inf,
                     noisyThresh=0.5, noisyProp=0.25,
                     rmMarkersNearZeroInAnySample::Bool=true) where {T}
  Js = size.(y, 2)
  J = Js[1]
  I = length(y)

  @assert all(J .== Js)
  @assert 0 <= subsample <= 1

  # Remove bad rows
  if !isinf(rowThresh)
    removeBadRows!(y, thresh=rowThresh)
  end

  goodColumns = [!isBadColumnForAllSamples(j, y,
                                           maxNanOrNegProp=maxNanOrNegProp,
                                           maxPosProp=maxPosProp) for j in 1:J]

  if rmMarkersNearZeroInAnySample
    goodColumns2 = [!markerIsNearZeroInAnySample(j, y,
                                                 thresh=noisyThresh,
                                                 prop=noisyProp)
                    for j in 1:J]
    goodColumns = goodColumns .& goodColumns2
  end

  if 0 < subsample < 1
    for i in 1:I
      Ni = size(y[i], 1)
      new_Ni = Int(round(Ni * subsample))
      # random_indices = rand(1:Ni, new_Ni)
      random_indices = Distributions.sample(1:Ni, new_Ni, replace=false)
      y[i] = y[i][random_indices, :]
    end
  end

  for i in 1:I
    y[i] = y[i][:, goodColumns]
  end

  return goodColumns, J
end

end # PreProcess


text:
using pfnet

case = ARGS[1]

parser = Parser(case)

net = parse_case(parser, case)

show_components(net)


text:
module Types

using CommonMark
export BasicSlide, SectionSlide, SlideView, Presentation, AbstractSlide
export meta, segment

abstract type AbstractSlide end

struct BasicSlide <: AbstractSlide
    meta::Union{Dict, Nothing}
    path::String
    root::CommonMark.Node
end

struct SectionSlide <: AbstractSlide
    meta::Union{Dict, Nothing}
    path::String
    slides::Vector{BasicSlide}
end

struct Presentation
    meta::Union{Dict, Nothing}
    path::String
    slides::Vector{AbstractSlide}
    standalone::Bool
end

meta(x::AbstractSlide) = x.meta
meta(x::Presentation) = x.meta

function Base.show(io::IO, m::MIME"text/plain", x::BasicSlide, env=Dict{String,Any}())
    # printstyled(io, "Slide:\n\n"; bold=true)
    show(io, m, x.root, env)
    return
end

function Base.show(io::IO, m::MIME"text/plain", x::SectionSlide, env=Dict{String,Any}())
    printstyled(io, "Section:\n\n"; bold=true)
    print_slides(io, m, x.slides, env, "---")
    return
end

function Base.show(io::IO, m::MIME"text/plain", x::Presentation, env=Dict{String,Any}())
    printstyled(io, "Presentation:\n\n"; bold=true)

    if haskey(x.meta, "authors")
        printstyled(io, "authors: ", join(x.meta["authors"], ", "))
        println(io, "\n")
    end

    printstyled(io, "="^(displaysize()[2]÷2); color=:light_black)
    println(io, "\n")
    print_slides(io, m, x.slides, env, "="^(displaysize()[2]÷2))
    return
end

function print_slides(io::IO, m::MIME"text/plain", slides, env, div="---")
    for (k, slide) in enumerate(slides)
        show(io, m, slide, env)

        if k != length(slides)
            println(io)
            printstyled(io, div * "\n\n"; color=:light_black)
        end
    end
end

end


text:
using Base.Test
using SQSChannels

@testset "test sqs channel" begin
    queueName = "testSQSChannel"
    c = SQSChannel( queueName )

    println("test purge queue ...")
    empty!(c)

    testmsg = "test message"

    println("test batch mode of put! ...")
    msgs = repeat([testmsg], outer = 24)
    put!(c, msgs)

    println("test put! a set of messages ...")
    msgCollection = Set("$i" for i in 1:5)
    put!(c, msgCollection)
    for i in 1:5
        msg = take!(c)
    end

    println("test put! a single message ...")
    put!(c, testmsg)
    sleep(10)

    println("test fetch message ...")
    handle, message = fetch(c)
    println("fetched message: $(message)")
    @test message == testmsg
end


text:
using RDatasets, Econometrics # аналог import numpy...
using Gadfly
using CSV
using GLM
using LinearAlgebra

psid = dataset("Ecdat", "PSID") # Panel Study of Income Dynamics
psid

psid.Age
psid."Age"
psid[!, "Age"]
psid[2:5, ["Age", "Educatn"]]

model_0 = fit(EconometricModel,
    @formula(log(1 + Earnings) ~ Age + Educatn), psid)

model_1 = fit(EconometricModel,
    @formula(log(1 + Earnings) ~ Age), psid)

# ftest(model_0)

vcov(model_0)
vcov(model_0, HC0)
vcov(model_0, HC2)

confint(model_0, level=0.9, se=sqrt.(diag(vcov(model_0, HC2))))

diag(vcov(model_0, HC2))

residuals(model_0)
fitted(model_0)
coef(model_0)

?vcov

diamonds = dataset("ggplot2", "diamonds") # diamonds

demo = DataFrame(x=[4, 5, 6], y=[6, missing, 9], z=5:7, w=9:2:13)

demo[2:3, ["x", "y"]]

missing

x = [4, 5, missing]

describe(demo)

demo[!, "x"] *= 10;

demo[!, "x"]

cos(sin(0))
0 |> sin |> cos

iris = dataset("datasets", "iris")

plot(iris, x=:SepalLength, y=:SepalWidth)
plot(iris, x=:SepalLength, Geom.bar)
plot(iris, x=:SepalLength, Geom.bar, Guide.xlabel("Длина чашелистика в мм"))

# import Pkg; Pkg.add("Gadfly")

CSV.write("demo.csv", demo)

aaa = CSV.read("demo.csv")

pwd()
cd("/home/boris/Documents/metrics_hse_2020_21/")

# using DataFramesMeta

# import Pkg; Pkg.add("DataFramesMeta")


text:
switch(x::T) where {T <: Edge} = T(x.dst, x.src)
switch(x::T) where {T <: ModifiedEdge} = T(switch(x.edge), x.modifiers)
undirect(x::DirectedEdge) = UndirectedEdge(x.src, x.dst)
undirect(x::ModifiedEdge) = ModifiedEdge(undirect(x.edge), x.modifiers)

@testset "Single Arrow" begin
    result = Edge(:a, :b)
    @test unarrow(:a → :b)[1] == result
    @test unarrow(:a ⇒ :b)[1] == result
    @test unarrow(:a ← :b)[1] == switch(result)
    @test unarrow(:a ⇐ :b)[1] == switch(result)
    @test unarrow(:a ↔ :b)[1] == undirect(result)
    @test unarrow(:a ⇔ :b)[1] == undirect(result)
    @test unarrow(:a ⇔ :b)[1] == undirect(switch(result))
end

@testset "Single Chain" begin
    result = [Edge(:a, :b), Edge(:b, :c)]
    @test unarrow(:a → :b → :c) == result
    @test unarrow(:a ⇒ :b ⇒ :c) == result
    @test unarrow(:a → :b ⇒ :c) == result

    @test unarrow(:a ← :b ← :c) == switch.(result)
    @test unarrow(:a ⇐ :b ← :c) == switch.(result)
    @test unarrow(:a ⇐ :b ⇐ :c) == switch.(result)

    @test unarrow(:a ↔ :b ↔ :c) == undirect.(switch.(result))
    @test unarrow(:a ⇔ :b ↔ :c) == undirect.(switch.(result))
    @test unarrow(:a ⇔ :b ⇔ :c) == undirect.(switch.(result))
end

@testset "Multiple Nodes Arrow" begin
    bresult = [Edge(:a, :c); Edge(:b, :d)]
    cresult = [Edge(:a, :c); Edge(:a, :d); Edge(:b, :c); Edge(:b, :d)]
    @test unarrow([:a :b] → [:c :d]) ⊆ bresult
    @test unarrow([:a :b] ← [:c :d]) ⊆ switch.(bresult)

    @test unarrow([:a :b] ⇒ [:c :d]) ⊆ cresult
    @test unarrow([:a :b] ⇐ [:c :d]) ⊆ switch.(cresult)

    @test unarrow([:a :b] ↔ [:c :d]) ⊆ undirect.(bresult)
    @test unarrow([:a :b] ⇔ [:c :d]) ⊆ undirect.(cresult)
end

@testset "Multiple Nodes Chain" begin
    bresult = [Edge(:a, :c); Edge(:b, :d); Edge(:c, :e); Edge(:d, :f); Edge(:e, :g); Edge(:f, :h)]
    cresult = [Edge(:a, :c); Edge(:b, :d); Edge(:c, :e); Edge(:d, :f); Edge(:e, :g); Edge(:f, :h);
    Edge(:b, :c); Edge(:a, :d); Edge(:d, :e); Edge(:c, :f); Edge(:f, :g); Edge(:e, :h)]

    @test unarrow([:a :b] → [:c :d] → [:e :f] → [:g :h]) ⊆ bresult
    @test unarrow([:a :b] ← [:c :d] ← [:e :f] ← [:g :h]) ⊆ switch.(bresult)

    @test unarrow([:a :b] ⇒ [:c :d] ⇒ [:e :f] ⇒ [:g :h]) ⊆ cresult
    @test unarrow([:a :b] ⇐ [:c :d] ⇐ [:e :f] ⇐ [:g :h]) ⊆ switch.(cresult)

    @test unarrow([:a :b] ↔ [:c :d] ↔ [:e :f] ↔ [:g :h]) ⊆ undirect.(bresult)
    @test unarrow([:a :b] ⇔ [:c :d] ⇔ [:e :f] ⇔ [:g :h]) ⊆ undirect.(cresult)
end

struct Start <: EdgeModifier
    s
end

struct Weight <: EdgeModifier
    w
end

@testset "Modified  Chain" begin
    # correct placement of modifiers
    result = ModifiedEdge(Edge(:a, :b), [Start(4)])
    @test unarrow(:a → Start(4)*:b)[1] == result
    @test unarrow(Start(4) * :a ← :b)[1] == switch(result)
    @test unarrow(:a ↔ :b * Start(4))[1] == undirect(result)

    #incorrect
    @test unarrow(Start(4) * :a → :b)[1] == result.edge
    @test unarrow(:a ← :b * Start(4))[1] == switch(result).edge
    @test unarrow(Start(4) * :a ↔ :b)[1] == undirect(result).edge

end


text:
using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

# Create a
function pendulum(u, p, t)
    x = u[2]
    y = -9.81sin(u[1]) - 0.1u[2]^3 -0.2*cos(u[1])
    return [x;y]
end

u0 = [0.99π; -1.0]
tspan = (0.0, 20.0)
prob = ODEProblem(pendulum, u0, tspan)
sol = solve(prob, Tsit5(), saveat = 0.3)

plot(sol)

# Create the differential data
DX = similar(sol[:,:])
for (i, xi) in enumerate(eachcol(sol[:,:]))
    DX[:,i] = pendulum(xi, [], 0.0)
end

# Create a basis
@variables u[1:2]
@parameters w[1:2]
# Lots of polynomials
polys = Operation[1]
for i ∈ 1:5
    push!(polys, u.^i...)
    for j ∈ 1:i-1
        push!(polys, u[1]^i*u[2]^j)
    end
end

# And some other stuff
h = [w[1]*cos(u[1]); w[2]*sin(u[1]); u[1]*u[2]; u[1]*sin(u[2]); u[2]*cos(u[2]); polys...]

basis = Basis(h, u, parameters = w)
println(basis)

# Get the reduced basis via the sparse regression
# Thresholded Sequential Least Squares, works fine for more data
# than assumptions, converges fast but fails sometimes with too much noise
opt = STRRidge(1e-2)
# Enforce all 100 iterations
Ψ = SINDy(sol[:,1:25], DX[:, 1:25], basis, opt, p = [1.0; 1.0], maxiter = 100, convergence_error = 1e-5)
println(Ψ)
print_equations(Ψ)

# Lasso as ADMM, typically needs more information, more tuning
opt = ADMM(1e-2, 1.0)
Ψ = SINDy(sol[:,1:50], DX[:, 1:50], basis, opt, p = [1.0; 1.0], maxiter = 5000, convergence_error = 1e-3)
println(Ψ)
print_equations(Ψ)

# Get the associated parameters out of the result
parameters(Ψ)

# SR3, works good with lesser data and tuning
opt = SR3(1e-2, 1.0)
Ψ = SINDy(sol[:,1:end], DX[:, 1:end], basis, opt, p = [0.5; 0.5], maxiter = 5000, convergence_error = 1e-5)
println(Ψ)
print_equations(Ψ, show_parameter = true)

# Vary the sparsity threshold -> gives better results
λs = exp10.(-5:0.1:-1)
# Use SR3 with high relaxation (allows the solution to diverge from LTSQ) and high iterations
opt = SR3(1e-2, 5.0)
Ψ = SINDy(sol[:,1:10], DX[:, 1:10], basis, λs, opt, p = [1.0; 1.0], maxiter = 15000)
println(Ψ)
print_equations(Ψ)

# Transform into ODE System
sys = ODESystem(Ψ)
dudt = ODEFunction(sys)
ps = parameters(Ψ)

# Simulate
estimator = ODEProblem(dudt, u0, tspan, ps)
sol_ = solve(estimator, Tsit5(), saveat = sol.t)

# Yeah! We got it right
scatter(sol.t[1:10], sol[:,1:10]', color = :red, label = nothing)
scatter!(sol.t[11:end], sol[:,11:end]', color = :blue, label = nothing)
plot!(sol_.t, sol_[:, :]', color = :green, label = "Estimation")

plot(sol.t, abs.(sol-sol_)')
norm(sol[:,:]-sol_[:,:], 2)


text:
# Run tests on functions in src/by_reference_doc/search.jl

@testset verbose = true "GET-request endpoints for search" begin

    @testset "For query = $(query)" for query in ["Coldplay", "Greenday", "Adele"]

        @test ~isempty(Spotify.search_get(;q = "$query", item_type = "album")[1])

    end

    @testset "For item type = $(item_type)" for item_type in ["album", "playlist", "track"]

        @test ~isempty(Spotify.search_get(;q = "Hans Zimmer", item_type = "$item_type")[1])

    end

end

# Test if 400 exception handling works by using an empty search query

@testset "Test 400 Not Found exception" begin

    @test_logs (:info, "Check if the search query is present") match_mode=:any Spotify.search_get(;q = "", item_type = "album")[1]

end


text:
using BinaryBuilder, Pkg

name = "MUMPS"
version = v"5.4.2" # <-- This is a lie, we're bumping to 5.4.2 to create a Julia v1.6+ release with experimental platforms

sources = [
  ArchiveSource("http://mumps.enseeiht.fr/MUMPS_5.4.1.tar.gz",
                "93034a1a9fe0876307136dcde7e98e9086e199de76f1c47da822e7d4de987fa8"),
  DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
mkdir -p ${libdir}
cd $WORKSPACE/srcdir/MUMPS*
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mumps_int32.patch

makefile="Makefile.G95.PAR"
cp Make.inc/${makefile} Makefile.inc

if [[ "${target}" == aarch64-apple-darwin* ]]; then
    # Fix the error:
    #     Type mismatch in argument ‘s’ at (1); passed INTEGER(4) to LOGICAL(4)
    FFLAGS=("-fallow-argument-mismatch")
fi

if [[ "${target}" == *-apple* ]]; then
  CFLAGS=("-fno-stack-check")
fi

make_args+=(OPTF=-O \
            CDEFS=-DAdd_ \
            LMETISDIR=${prefix} \
            IMETIS="-I${includedir}" \
            LMETIS="-L${libdir} -lparmetis -lmetis" \
            ORDERINGSF="-Dpord -Dparmetis" \
            CC="mpicc -fPIC ${CFLAGS[@]}" \
            FC="mpif90 -fPIC ${FFLAGS[@]}" \
            FL="mpif90 -fPIC" \
            SCALAP="${libdir}/scalapack32.${dlext} ${libdir}/libopenblas.${dlext}" \
            INCPAR= \
            LIBPAR=-lscalapack32 \
            LIBBLAS=-lopenblas)

# Options for SCOTCH
# LSCOTCHDIR=${prefix}
# ISCOTCH="-I${includedir}"
# LSCOTCH="-L${libdir} -lesmumps -lscotch -lscotcherr"
# ORDERINGSF="-Dpord -Dparmetis -Dscotch"

if [[ "${target}" == *-apple* ]]; then
  make_args+=(RANLIB=echo)
fi

# NB: parallel build fails
make all "${make_args[@]}"

# build shared libs
all_load="--whole-archive"
noall_load="--no-whole-archive"
extra=""
if [[ "${target}" == *-apple-* ]]; then
    all_load="-all_load"
    noall_load="-noall_load"
    extra="-Wl,-undefined -Wl,dynamic_lookup -headerpad_max_install_names"
fi

cd lib
libs=(-lparmetis -lmetis -lscalapack32 -lopenblas)
mpif90 -fPIC -shared -Wl,${all_load} libpord.a ${libs[@]} -Wl,${noall_load} ${extra[@]} -o libpord.${dlext}
cp libpord.${dlext} ${libdir}

libs+=(-L${libdir} -lpord)
mpif90 -fPIC -shared -Wl,${all_load} libmumps_common.a ${libs[@]} -Wl,${noall_load} ${extra[@]} -o libmumps_common.${dlext}
cp libmumps_common.${dlext} ${libdir}

libs+=(-lmumps_common)
for libname in cmumps dmumps smumps zmumps
do
  mpif90 -fPIC -shared -Wl,${all_load} lib${libname}.a ${libs[@]} -Wl,${noall_load} ${extra[@]} -o lib${libname}.${dlext}
done
cp *.${dlext} ${libdir}
cd ..

cp include/* ${prefix}/include
"""

# OpenMPI and MPICH are not precompiled for Windows
# MUMPS doesn't build on PowerPC
platforms = expand_gfortran_versions(filter!(p -> !Sys.iswindows(p) && arch(p) != "powerpc64le", supported_platforms()))

# The products that we will ensure are always built
products = [
    LibraryProduct("libsmumps", :libsmumps),
    LibraryProduct("libdmumps", :libdmumps),
    LibraryProduct("libcmumps", :libcmumps),
    LibraryProduct("libzmumps", :libzmumps)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="MPICH_jll", uuid="7cb0a576-ebde-5e09-9194-50597f1243b4")),
    Dependency(PackageSpec(name="METIS_jll", uuid="d00139f3-1899-568f-a2f0-47f597d42d70")),
    Dependency(PackageSpec(name="SCOTCH_jll", uuid="a8d0f55d-b80e-548d-aff6-1a04c175f0f9")),
    Dependency(PackageSpec(name="PARMETIS_jll", uuid="b247a4be-ddc1-5759-8008-7e02fe3dbdaa")),
    Dependency(PackageSpec(name="SCALAPACK32_jll", uuid="aabda75e-bfe4-5a37-92e3-ffe54af3c273")),
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2")),
    # Dependency(PackageSpec(name="SCOTCH_jll", uuid="a8d0f55d-b80e-548d-aff6-1a04c175f0f9"))
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
#Minimum cover EMS method
using JuMP,Gurobi

#define the model
m =Model(Gurobi.Optimizer)

#Let n represents the number of decision variables
n=10

#declaring variables
@variable(m, x[1:n], Bin) # binary constraint

#define the objective function
@objective(m , Min, sum(x[i] for i =1:n)) # objective is to minimize the coverage

#Since decision variables x2,x3,x8 and x10>=1, I declared them in a for-loop together
for i in [2,3,8,10]
	@constraint(m,x[i]>=1)    #(x[i] for i in [2,3,8,10]) >=1)
end

#additional constraint
@constraint(m, x[1]+x[2]>=1)    #district 2
@constraint(m, x[1]+x[3]>=1)    #district 3
@constraint(m, x[2]+x[4]>=1)	#district 7
@constraint(m, x[3]+x[4]>=1)	#district 8
@constraint(m, x[4]+x[6]>=1)	#district 10
@constraint(m, x[4]+x[5]>=1)	#district 11
@constraint(m, x[4]+x[5]+x[7]>=1)	#distict 13
@constraint(m, x[4]+x[5]+x[6]>=1)	#distict 12
@constraint(m, x[8]+x[9]>=1)		#distict 14, district 18
@constraint(m, x[6]+x[9]>=1)		#distict 15
@constraint(m, x[5]+x[6]>=1)		#distict 16
@constraint(m, x[5]+x[7]+x[10]>=1)	#distict 17
@constraint(m, x[8]+x[10]>=1)		#distict 19

print(m)

#Printing the optimal solutions obtained
println("Optimal Solutions:")
status = JuMP.optimize!(m)
println(" Objective value: ",JuMP.objective_value(m))

println("Decision variable values:")
for i=1:10
  println("x[$i] = ",   JuMP.value(x[i]))
end

max_stars_repo_name: Oyekanmi/Optimization-Julia

max_issues_repo_name: Oyekanmi/Optimization-Julia

max_forks_repo_name: Oyekanmi/Optimization-Julia


text:
@testset "harmonic transformer model" begin
    @testset "YN - yn - 0" begin
        path  = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_xfmr_YNyn0.m")
        data  = _PMs.parse_file(path)
        hdata = _HPM.replicate(data)

        run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
    end

    @testset "D - yn - 11" begin
        path  = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_xfmr_Dyn11.m")
        data  = _PMs.parse_file(path)
        hdata = _HPM.replicate(data)

        run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
    end

    @testset "Y - d - 11" begin
        path  = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_xfmr_Yd11.m")
        data  = _PMs.parse_file(path)
        hdata = _HPM.replicate(data)

        run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
    end

    @testset "YN - d - 11" begin
        path  = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_xfmr_YNd11.m")
        data  = _PMs.parse_file(path)
        hdata = _HPM.replicate(data)

        run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
    end
end


text:
function assemble_cell_body_force!(
    systemrhs,
    rhsfunc,
    basis,
    cellquads,
    mesh,
    cellsign,
    cellid,
)

    dim = dimension(mesh)
    detjac = determinant_jacobian(mesh)
    cellmap = cell_map(mesh, cellsign, cellid)
    quad = cellquads[cellsign, cellid]
    rhs = linear_form(rhsfunc, basis, quad, cellmap, dim, detjac)
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_rhs!(systemrhs, nodeids, dim, rhs)
end

function assemble_body_force!(systemrhs, rhsfunc, basis, cellquads, mesh)

    ncells = number_of_cells(mesh)
    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        check_cellsign(cellsign)

        if cellsign == +1 || cellsign == 0
            assemble_cell_body_force!(
                systemrhs,
                rhsfunc,
                basis,
                cellquads,
                mesh,
                +1,
                cellid,
            )
        end
        if cellsign == -1 || cellsign == 0
            assemble_cell_body_force!(
                systemrhs,
                rhsfunc,
                basis,
                cellquads,
                mesh,
                -1,
                cellid,
            )
        end
    end
end


text:
using Test
using EllZeroTrendFiltering

for problem_fcn in ["discontinuous1",
                    "discontinuous2",
                    "OEIS1",
                    "super_exponential1",
                    "super_exponential2",
                    "super_exponential3",
                    "super_exponential4",
                    "super_exponential5",
                    "square_wave",
                    "white_noise",
                    "exponential",
                    "straight_line"]

    include(joinpath(dirname(@__FILE__),"problems", problem_fcn * ".jl"))
    local g, I_sols  # TODO, can remove in julia 1.0?
    g, ζ_vec, I_sols, f_sols = @eval $(Symbol(problem_fcn))()

    M = length(I_sols)
    I_vec, _, f_vec = fit_pwl_constrained(g, M)

    @testset "Data set: $problem_fcn, constrained fit with m=$m" for m in 1:length(f_sols)
        @test f_vec[m] ≈ f_sols[m]  atol=1e-10
        if !isempty(I_sols[m])
            @test I_vec[m] == I_sols[m]
        end
    end

    # Iterate only over problems where solution exists
    @testset "Data set: $problem_fcn, subsampled, constrained fit with m=$m" for m in 1:length(f_sols)
        # Generate max(length(g)/5, m) random gridpoints in range 1:length(g)
        randI = shuffle(1:length(g))[1:max(floor(Int,length(g)/5), m+1)]
        # Merge with correct solution, and [1,length(g)] if solution doesn't exist
        t_grid = sort(union(randI, I_sols[m], [1, length(g)]))
        # Solve on subsampled grid,
        I_vec, _, f_vec = fit_pwl_constrained(g, m, t = t_grid)
        # Should be same if we had a solution
        if !isempty(I_sols[m])
            @test f_vec[m] ≈ f_sols[m]  atol=1e-10
            @test t_grid[I_vec[m]] == I_sols[m]
        else # Or higher cost if not
            @test f_vec[m] >= f_sols[m] || isapprox(f_vec[m], f_sols[m], atol=1e-10)
        end
    end
end


text:
##  Continuous wavelet analysis in julia
##  Author Bernd Blasius (bernd.blasius@gmail.com)

# for a more involved example, including plotting, wavelet power spectrum
# and significance levels, see test/cwt_test_elnino.jl

using Statistics
include("wavelet.jl")

# continuoous wavelet transform (cwt)
function cwt(x,dt;
             pad=1,        # pad the time series with zeros
             nvoices = 32, # total number of voices
             noctave = 4,  # total number of octaves
             s0 = 2*dt,    # starting scale
             mother = "MORLET",
             param = 6)
    # x : vector of time series
    # dt: sampling time interval
    # this function assumes that the time series is sampled equally
    # no checks are made to ensure this!!

    dj = 1 ./ nvoices
    j1 = noctave * nvoices

    # normalize time series, but store the variance
    variance = var(x)
    x = x .- mean(x)
    x = x ./ std(x)

    wave, period, scale, coi = wavelet(x,dt,pad,dj,s0,j1,mother,param)

    # wavelet power spetrum
    power = variance * abs2.(wave)

    # Global wavelet spectrum
    global_ws = variance * mean(power,dims=2)   # time-average over all times

    return wave, period, power, global_ws, variance, scale, coi
end

max_stars_repo_name: berndblasius/WaveletAnalysis

max_issues_repo_name: berndblasius/WaveletAnalysis

max_forks_repo_name: berndblasius/WaveletAnalysis


text:
import Interpolations:
    BoundaryCondition,
    Flat

import ..UncertainDatasets:
    AbstractUncertainValueDataset,
    UncertainIndexValueDataset

import ..SamplingConstraints:
    TruncateQuantiles,
    constrain

const EXTRAP_BC = Union{BoundaryCondition, Float64}

##############################################################################
# Rectangular grid composite types
##############################################################################
const BCC = Union{BoundaryCondition, Float64}
"""
    RegularGrid

## Fields
- **min**: The minimum value of the grid.
- **max**: The maximum value of the grid.
- **step**: The interval size.
- **extrapolation_bc**: The extrapolation condition. Can also be NaN.
"""
struct RegularGrid <: InterpolationGrid
    min::Number
    max::Number
    step::Number
    extrapolation_bc::BCC
end

function get_edgepoints(x::RegularGrid)
    x.min:x.step:x.max
end

Base.length(x::RegularGrid) = length(x.min:x.step:x.max)

##############################################################################
# Rectangular grid constructors for different types of uncertain datasets
##############################################################################

function RegularGrid(x::AbstractUncertainValueDataset, n_steps::Int;
    extrapolation_bc::EXTRAP_BC = Flat(), round_digits::Int = 7,
    trunc_finite::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

    # Constrain to some large quantile to avoid taking the minima and
    # maximuma of distributions with infinite supports
    constrained_indices = constrain(x, TruncateQuantiles(0.0001, 0.9999))

    min_x = floor(minimum(constrained_indices), digits = round_digits)
    max_x = ceil(maximum(constrained_indices), digits = round_digits)

    RegularGrid(min_x, max_x, (max_x - min_x)/n_steps, extrapolation_bc)
end

function RegularGrid(x::UncertainIndexValueDataset, n_steps::Int;
    extrapolation_bc::EXTRAP_BC = Flat(), round_digits::Int = 7,
    trunc_finite::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

    # Constrain to some large quantile to avoid taking the minima and
    # maximuma of distributions with infinite supports
    constrained_indices = constrain(x.indices, TruncateQuantiles(0.0001, 0.9999))

    min_x = floor(minimum(constrained_indices), digits = round_digits)
    max_x = ceil(maximum(constrained_indices), digits = round_digits)

    RegularGrid(min_x, max_x, (max_x - min_x)/n_steps, extrapolation_bc)
end

function RegularGrid(x::AbstractUncertainValueDataset, step::Float64;
    extrapolation_bc::EXTRAP_BC = Flat(),
    round_digits = 7,
    trunc_finite::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

    # Constrain to some large quantile to avoid taking the minima and
    # maximuma of distributions with infinite supports
    constrained_data = constrain(x, trunc_finite)

    min_x = floor(minimum(constrained_data), digits = round_digits)
    max_x = ceil(maximum(constrained_data), digits = round_digits)

    RegularGrid(min_x, max_x, step, extrapolation_bc)
end

function RegularGrid(x::UncertainIndexValueDataset, step::Float64;
    extrapolation_bc::EXTRAP_BC = Flat(), round_digits::Int = 7,
    trunc_finite::TruncateQuantiles = TruncateQuantiles(0.001, 0.999))

    # Constrain to some large quantile to avoid taking the minima and
    # maximuma of distributions with infinite supports
    constrained_indices = constrain(x.indices, TruncateQuantiles(0.0001, 0.9999))

    min_x = floor(minimum(constrained_indices), digits = round_digits)
    max_x = ceil(maximum(constrained_indices), digits = round_digits)

    RegularGrid(min_x, max_x, step, extrapolation_bc)
end

RegularGrid(xmin, xmax, step; extrapolation_bc::EXTRAP_BC = Flat()) =
    RegularGrid(xmin, xmax, step, extrapolation_bc)

function RegularGrid(xvals::AbstractRange; extrapolation_bc::EXTRAP_BC = Flat())
    RegularGrid(minimum(xvals), maximum(xvals), step(xvals), extrapolation_bc)
end

function RegularGrid(xvals::Vector{T}; extrapolation_bc::EXTRAP_BC = Flat()) where {T}
    RegularGrid(minimum(xvals), maximum(xvals), xvals[2] - xvals[1], extrapolation_bc)
end

export RegularGrid


text:
module MKL

const debug = true

if debug println("MKL : Intel MKL wrapper") end

Libdl.dlopen(:libmkl_rt) # necessary to preload - don't know why
const libmkl   = :libmkl_rt
const libvml   = :libmkl_vml_avx

# VML
include("vml_service.jl")
include("vml_arithmetic.jl")
include("vml_power_root.jl")

# BLAS
include("blas_level_1.jl")
include("blas_level_2.jl")

# LAPACK

end


text:
function calc_Etot( pw::PWGrid, Vpot, psi::Array{ComplexF64,2} )

    Ω = pw.Ω
    Npoints = pw.Npoints
    Nstates = size(psi)[2]

    Kpsi = op_K( pw, psi )
    Ekin = 0.0
    for is = 1:Nstates
        Ekin = Ekin + real( dot( psi[:,is], Kpsi[:,is] ) )
    end

    # Calculate in real space
    rho = calc_rho( pw, psi )
    Epot = dot( rho, Vpot ) * Ω/Npoints

    Etot = Ekin + Epot
    return Etot
end