text:
function verify(imagelist)
  dims = ""
  depth = ""
  mode = ""
  ptr = 1
  while ptr <= length(imagelist)
    wand = MagickWand()
    ImageMagick.pingimage(wand, imagelist[ptr])
    dmode = getimagecolorspace(wand)
    ddepth = getimagedepth(wand)
    ddims = size(wand)
    if mode == ""
      mode = dmode
    elseif mode != dmode
      println("Mode of ", imagelist[ptr], " is ", dmode, ", which doesn't match mode ", mode, " of the first image.")
      println("Exiting")
      exit(-1)
    end

    if depth == ""
      depth = dmode
    elseif depth != dmode
      println("Depth of ", imagelist[ptr], " is ", dmode, ", which doesn't match depth ", depth, " of the first image.")
      println("Exiting")
      exit(-1)
    end

    if dims == ""
      dims = ddims
    elseif dims != ddims
      println("Shape of ", imagelist[ptr], " is ", ddims, ", which doesn't match shape ", shape, " of the first image.")
      println("Exiting")
      exit(-1)
    end
    if verbose
      println("verified ", imagelist[ptr])
    end
    ptr += 1
  end
end

function getdepth(str)
  str = string(str)
  depth = match(r"(UInt\d+)", str)
  if depth == nothing
    println("Error: failed to determine depth from ", str)
    exit(-1)
  end
  if depth[1] == "UInt8"
    rslt = UInt8
  elseif depth[1] == "UInt16"
    rslt = UInt16
  else
    println("Can't handle color depth ", depth[1])
    exit(-1)
  end
  return rslt
end

arglist = ARGS

approach = 1 # mean by default
nomen = "mean"
verbose = false
oname = ""
filelist = String[]
while length(arglist) > 0
  opt = shift!(arglist)
  if opt == "-mean"
    approach = 1
    nomen = "mean"
  elseif opt == "-max"
    approach = 2
    nomen = "max"
  elseif opt == "-min" # why?
    approach = 3
    nomen = "min"
  elseif opt == "-sigma" # sigma
    approach = 4
    nomen = "sigma clipped"
    sigma = 1 # will overwrite in a minute if specified in the command line
    if length(arglist) > 0 # see if they gave an argument
      text = shift!(arglist)
      x = tryparse(Float16, text)
      if isnull(x)
        unshift!(arglist, text)
      else
        sigma = x
        approach = 4
      end
    end
  elseif opt == "-o"
    if length(arglist) > 0
      oname = shift!(arglist)
    else
      println("-o specified but no filename")
      exit(-1)
    end
  elseif opt == "-v"
    verbose = true
  elseif opt == "-"
    println("Unknown option $opt")
    exit(-1)
  else
    if length(arglist) == 0 && oname == "" # last name is output file if not otherwise specified
      oname = opt
    else
      filelist = push!(filelist, opt)
    end
  end
end

if verbose
  println("Verifying that input files exist... ")
end

for fname in filelist
  if !isfile(fname)
    println("File does not exist: $fname")
    exit(-1)
  end
end

if verbose
  println("OK!")
end

using ImageMagick

if verbose
  print("\nVerifying that input file formats match... ")
end

verify(filelist)

if verbose
  println("OK!")
end

println("\nApplying ", nomen, " to ", length(filelist), " input files and saving to ", oname)

using Images
names = filelist
first = true
what = nothing
sums = nothing
count_ = length(filelist)
while length(names) > 0
  fname = shift!(names)
  if verbose
    print("Loading ", fname , "... ")
  end
  img =  load(fname)
  if isnull(img)
    println("$fname failed to load")
    exit(-1)
  end
  sz = size(img)
  if first
    szx = sz[1]
    szy = sz[2]
    what = getdepth(typeof(img))
    if verbose
      println("First image is $szx x $szy")
    end
  elseif sz[1] != szx || sz[2] != szy
    println("$name has size ", sz[1], "x", sz[2], ", which doesn't match initial size $szx x $szy")
    exit(-1)
  end

  img = reinterpret(what, img)
  img = img[1:3, :, :]
  println(size(img), " ", typeof(img))
  if verbose
    print("converting to floats... ")
  end
  img = map(Float32, img)
  if verbose
    print("merging... ")
  end
  if approach == 1 || approach == 4
    if first
      sums = zeros(Float32, size(img))
    end
    sums .+= img
  elseif approach == 2
    if first
      maxes = img
    else
      mat = img .> maxes
      maxes[mat] = img[mat]
    end
  elseif approach == 3
    if first
      mins = img
    else
      mat = img .< mins
      mins[mat] = img[mat]
    end
  end
  if approach == 4
    if first
      sumsq = zeros(Float32, size[img])
    end
    sumsq .+= img .^ 2
  end
  first = false
  if verbose
    println("Done!")
  end
end

if approach == 1
  println(size(sums), " ", typeof(sums))
  sums ./= count_
  if what == UInt16
    sums = map(x -> floor(UInt16, x), sums)
    img = colorview(RGB{N0f16}, sums)
  elseif what == UInt8
    img = map(x -> floor(UInt8, x), sums)
    img = colorview(RGB{N0f8}, img)
  end
  save(oname, img)
elseif approach == 2
  save(oname, maxes)
elseif approach == 3 # why?
  save(oname, mins)
elseif approach == 4 #
  println("Um ... this doesn't work yet")
else
  println("How did approach get set to $approach?")
  exit(-1)
end

max_stars_repo_name: FatherGeoffHorton/AstroTools

max_issues_repo_name: FatherGeoffHorton/AstroTools

max_forks_repo_name: FatherGeoffHorton/AstroTools


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Foundations of machine learning")

## Run this cell to load the graphics packages
using Plots; gr()

# ## Adding a function parameter

#-

# In the last notebook, we saw an example of adding **parameters** to functions. These are values that control the behavior of a function. Let's look at that in some more detail.

#-

# Let's go back to our original version of the σ function:

σ(x) = 1 / (1 + exp(-x))

# Instead of working with a single function, we can work with a whole class (set) of functions that look similar but differ in the value of a **parameter**. Let's make a new function that uses the previous $\sigma$ function, but also has a parameter, $w$. Mathematically, we could write
#
# $$f_w(x) = f(x; w) = \sigma(w \, x).$$
#
# (Here, $w$ and $x$ are multiplied in the argument of $\sigma$; we could write $w \times x$ or $w * x$, or even $w \cdot x$, but usually the symbols are not written.)
#
# Mathematically speaking, we can think of $f_w$ as a different function for each different value of the parameter $w$.

#-

# In Julia, we write this as follows:

f(x, w) = σ(w * x)

# Note that Julia just treats parameters as additional *arguments* to the function; the function `f` has two arguments, the value of `x` and the value of `w` that we want to use.

#-

# We can now investigate the effect of $w$ interactively. To do so, we need a way of writing in Julia "the function of one variable $x$ that we obtain when we fix the value of $w$". We write this as an "anonymous function", as we saw in the notebook on functions:
#
#     x -> f(x, w)
#
# We can read this as "the function that maps $x$ to the value of $f(x, w)$, for a value of $w$ that was previously given".

#-

# Now we are ready to draw the function. For each plot, we *fix* a value of the parameter $w$ and draw the resulting function as a function of $x$. We can interactively modify the value of $w$, and plot the new function that comes out:

## Try manipulating `w` between -2 and 2
w = 0.0
plot(x->f(x, w), -5, 5, ylims=(0,1), label="sigmoid")
plot!(x->(x>0), -5,5, label="Square Wave")

# #### Exercise 1
#
# Try writing your own function that takes a parameter. Start by copying and executing
#
# ```julia
# square(x) = x^2
# ```
#
# Then use `square` to declare a new function `square_and_scale` that takes two inputs, `a` and `x` such that
#
# $$\mathrm{square\_and\_scale}(x; a) := a \cdot x^2$$

#-

# #### Solution
#
# You should have declared the following function:
#
# ```julia
# square_and_scale(x, a) = a * x^2
# ```
#
# We can write some tests as follows, using the approximate equality operator, `≈`, written as `\approx<TAB>`:
#
# ```julia
# square_and_scale(1, 1) ≈ 1
# square_and_scale(1, 2) ≈ 2
# square_and_scale(2, 1) ≈ 4
# square_and_scale(2, 2) ≈ 8
# ```

#-

# #### Exercise 2

#-

# Once you have declared `square_and_scale`, uncomment the code below and see how the parameter `a` scales the function `square` :

## x = -10:10
## a = 5.0 # Try manipulating a here
## plot(x, square.(x), label="x^2")
## plot!(x, square_and_scale.(x, a), ls=:dash, label="ax^2")

# ## Fitting a function to data

#-

# As we saw in the previous notebook, what we would like to do is use the fact that we now have a parameter in our function in order to do something useful! Namely, we want to model data with it.

#-

# So suppose that we are given a single data point $(x_0, y_0) = (2, 0.8)$. We can try to "fit" the function $f_w$ by adjusting the parameter $w$ until the function passes through the data.
#
# **Game**: Change `w` until the graph of the function hits the data point. Which value of $w$ does that correspond to?

x0, y0 = 2, 0.8

w = 0.0 # try mainipulatoing `w` between -2 and 2
plot(x->f(x, w), -5, 5, ylims=(0, 1), label="f")
scatter!([x0], [y0], label="data")

# ## Quantifying how far we are from the goal: the *loss function*

#-

# We can see visually when the graph of the function passes through the data point. But the goal is to be able to automate this process so that the computer can do it unaided.
#
# So we will need a more precise way of deciding and quantifying (i.e. measuring with a number) *how far away we are from the goal*; here, the goal means hitting the data point with the function.

#-

# #### Exercise 3
#
# Can you think of a way of measuring how far away our function is from the data point?

#-

# ### Defining the loss function

#-

# We need to measure how far away the curve is from the data point when we choose a particular value of $w$.
# One way to do this is by finding the vertical distance $d$ from the curve to the data point.
#
# Instead of just taking the distance, it is common to take the *square* of the distance, $d^2$.
#
# Since we are taking the vertical distance, we need the distance at the given value of $x_0$ where the data point lies. For a given value of the parameter $w$, the height of the point on the curve with that value of $x_0$ is $f(x_0, w)$.
#
# So we take
# $$d := y_0 - f(x_0, w)$$
#
# and
# $$d^2 = [y_0 - f(x_0, w)]^2.$$
#
# This is our measure of distance. It will change when $w$ changes -- in other words, it is itself a *function of $w$*; we will denote this function by $L(w)$, and call it the **loss function**:
#
# $$L(w) := [y_0 - f(x_0, w)]^2.$$
#
# So the goal is to find the value $w^*$ of $w$ where the loss function is *least*; in other words, we need to *minimize* the loss function!
#
# (Another name for a loss function is a *cost function*.)

#-

# #### Exercise 4
#
# (a) Define the loss function `L(w)` in Julia.
#
# (b) Draw the data point and the function `x -> f(x, w)`. Also draw a vertical line from the data point to the function `x -> f(x, w)`.
#
# (c) Make the plot interactive.
#
# (d) Add as the plot title the value of the loss function for the current value of $w$.
#
# (e) Use the slider to find the value $w^*$ of $w$ for which the loss function reaches its minimum value. What is $w^*$? What is the value of the loss function there, $L(w^*)$?

#-

# #### Solution

#-

# (a) We have data

x0, y0 = 2, 0.8

# We define the loss function as follows:

L(w) = (y0 - f(x0, w))^2

# (b-d): The visualization can be written as follows

w = 0.0 # try manipulating w here between -2 and 2

plot(x->f(x, w), -5, 5, ylims=(0, 1), label="f_w", lw=3)  # function

scatter!([x0], [y0], label="data")
plot!([x0, x0], [y0, f(x0, w)], label="loss")  # vertical line

title!("L(w) = $(round(L(w); sigdigits = 5))")

# ## What does the loss function look like?

#-

# The loss function $L(w)$ tells us how far away the function $f_w$ is from the data when the parameter value is $w$, represented visually as the vertical line in the previous plot. When the data are fixed, this is a function only of the parameter $w$. What does this function look like as a function of $w$? Let's draw it!

#-

# #### Exercise 5
#
# Draw the function $L(w)$ as a function of $w$.

#-

# #### Solution

plot(L, -3, 3, xlabel="w", ylabel="L(w)", ylims=(0, 0.7), label="L")

# ### Features of the loss function

#-

# This graph quantifies how far we are from the data point for a given value of $w$.
# What features can we see from the graph?
#
# Firstly, we see that $L(w)$ is always bigger than $0$, for any value of $w$. This is because we want $L$ to be some kind of measure of *distance*, and distances cannot be negative.
#
# Secondly, we see that there is a special value $w^*$ of $w$ where the function $L$ reaches its minimum value. In this particular case, it actually reaches all the way down to $0$!
# This means that the original function $f$ (the one we manipulated above) passes exactly through the data point $(x_0, y_0)$.

#-

# #### Exercise 6
#
# Draw a zoomed-in version of the graph to find the place $w^*$ where the function hits $0$ more precisely.

#-

# ### A different way of defining the loss function

#-

# **Why did we use such a complicated function $L$ with those squares inside?** We could instead just have used the absolute distance, instead of the distance squared, using the *absolute value* function, written mathematically as $| \cdot |$, and in Julia as `abs`.

#-

# #### Exercise 7
#
# Define a new loss function, `L_abs`, using the absolute value, and see what it looks like.

#-

# #### Solution

L_abs(w) = abs(y0 - f(x0, w))

plot(L_abs, -3, 3, xlabel="w", ylabel="L_abs(w)", ylims=(0, 0.7), label="L_abs")

# Now we see why it was previously generally preferred to use squares: using the absolute value gives a cost function that is *not smooth*. This makes it difficult to use methods from calculus to find the minimum. Nonetheless, using non-smooth functions is very common in machine learning nowadays.

ext: jl

max_stars_repo_path: Courses/Foundations of machine learning/0500.Building-models-quiz.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# Differential calculus mimetic operators

"""
    curl!(q::Edges{Primal},w::Nodes{Dual})

Evaluate the discrete curl of `w` and return it as `q`.

# Example

```jldoctest
julia> w = Nodes(Dual,(8,6));

julia> w[3,4] = 1.0;

julia> q = Edges(Primal,w);

julia> curl!(q,w)
Edges{Primal,8,6,Float64} data
u (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  -1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
v (in grid orientation)
6×7 Array{Float64,2}:
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  -1.0  1.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function curl!(edges::Edges{Primal, NX, NY},
               s::Nodes{Dual,NX, NY}) where {NX, NY}

    @inbounds for y in 1:NY-1, x in 1:NX
        edges.u[x,y] = s[x,y+1] - s[x,y]
    end

    @inbounds for y in 1:NY, x in 1:NX-1
        edges.v[x,y] = s[x,y] - s[x+1,y]
    end
    edges
end

"""
    curl(w::Nodes{Dual}) --> Edges{Primal}

Evaluate the discrete curl of `w`. Another way to perform this operation is
to construct a `Curl` object and apply it with `*`.

# Example

```jldoctest
julia> C = Curl();

julia> w = Nodes(Dual,(8,6));

julia> w[3,4] = 1.0;

julia> C*w
Edges{Primal,8,6,Float64} data
u (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  -1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0  0.0
v (in grid orientation)
6×7 Array{Float64,2}:
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0  -1.0  1.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0
```
"""
curl(nodes::Nodes{Dual,NX,NY,T}) where {NX,NY,T} = curl!(Edges(Primal, nodes), nodes)

function curl!(nodes::Nodes{Dual,NX, NY},
               edges::Edges{Primal, NX, NY}) where {NX, NY}

    u, v = edges.u, edges.v
    @inbounds for y in 2:NY-1, x in 2:NX-1
        nodes[x,y] = u[x,y-1] - u[x,y] - v[x-1,y] + v[x,y]
    end
    nodes
end

function curl(edges::Edges{Primal, NX, NY}) where {NX, NY}
    curl!(Nodes(Dual,edges), edges)
end

struct Curl end

(*)(::Curl,w::Union{Nodes{Dual,NX,NY},Edges{Primal,NX,NY}}) where {NX,NY} = curl(w)

"""
    divergence!(w::Nodes,q::Edges)

Evaluate the discrete divergence of edge data `q` and return it as nodal data `w`.
Note that `q` can be either primal or dual edge data, but `w` must be of the same
cell type.

# Example

```jldoctest
julia> q = Edges(Primal,(8,6));

julia> q.u[3,2] = 1.0;

julia> w = Nodes(Primal,(8,6));

julia> divergence!(w,q)
Nodes{Primal,8,6,Float64} data
Printing in grid orientation (lower left is (1,1))
5×7 Array{Float64,2}:
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  1.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
```
"""
function divergence!(nodes::Nodes{Primal, NX, NY},
                     edges::Edges{Primal, NX, NY}) where {NX, NY}

    u, v = edges.u, edges.v

    @inbounds for y in 1:NY-1, x in 1:NX-1
        nodes[x,y] = - u[x,y] + u[x+1,y] - v[x,y] + v[x,y+1]
    end
    nodes
end

function divergence!(nodes::Nodes{Dual, NX, NY},
                     edges::Edges{Dual, NX, NY}) where {NX, NY}

    u, v = edges.u, edges.v

    @inbounds for y in 2:NY-1, x in 2:NX-1
        nodes[x,y] = - u[x-1,y] + u[x,y] - v[x,y-1] + v[x,y]
    end
    nodes
end

"""
    divergence(q::Edges) --> Nodes

Evaluate the discrete divergence of edge data `q`. Can also perform this operation
by creating an object of Divergence type and applying it with `*`.

# Example

```jldoctest
julia> D = Divergence();

julia> q = Edges(Primal,(8,6));

julia> q.u[3,2] = 1.0;

julia> D*q
Nodes{Primal,8,6,Float64} data
Printing in grid orientation (lower left is (1,1))
5×7 Array{Float64,2}:
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  1.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
```
"""
function divergence(edges::Edges{T, NX, NY}) where {T <: CellType, NX, NY}
    divergence!(Nodes(T, edges), edges)
end

"""
    divergence!(w::XEdges/YEdges,q::NodePair)

Evaluate the discrete divergence of node pair data `q` and return it as data `w`.
Note that `q` can be either primal/dual or dual/primal edge data, and `w`
must be, respectively, primal x-edges/dual y-edges or primal y-edges/dual x-edges type.
"""
function divergence!(nodes::Union{XEdges{Primal, NX, NY},YEdges{Dual, NX, NY}},
                     edges::NodePair{Primal, Dual,NX, NY}) where {NX, NY}

    u, v = edges.u, edges.v

    @inbounds for y in 1:NY-1, x in 2:NX-1
        nodes[x,y] = - u[x-1,y] + u[x,y] - v[x,y] + v[x,y+1]
    end
    nodes
end

function divergence!(nodes::Union{YEdges{Primal, NX, NY},XEdges{Dual, NX, NY}},
                     edges::NodePair{Dual, Primal,NX, NY}) where {NX, NY}

    u, v = edges.u, edges.v

    @inbounds for y in 2:NY-1, x in 1:NX-1
        nodes[x,y] = - u[x,y] + u[x+1,y] - v[x,y-1] + v[x,y]
    end
    nodes
end

struct Divergence end

(*)(::Divergence,w::Edges{T,NX,NY}) where {T<:CellType,NX,NY} = divergence(w)

"""
    grad!(q::Edges{Primal},w::Nodes{Primal})

Evaluate the discrete gradient of primal nodal data `w` and return it as primal
edge data `q`.

# Example

```jldoctest
julia> w = Nodes(Primal,(8,6));

julia> w[3,4] = 1.0;

julia> q = Edges(Primal,(8,6));

julia> grad!(q,w)
Edges{Primal,8,6,Float64} data
u (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
v (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
```
"""
function grad!(edges::Edges{Primal, NX, NY},
                   p::Nodes{Primal, NX, NY}) where {NX, NY}

    @inbounds for y in 1:NY-1, x in 2:NX-1
        edges.u[x,y] = - p[x-1,y] + p[x,y]
    end
    @inbounds for y in 2:NY-1, x in 1:NX-1
        edges.v[x,y] = - p[x,y-1] + p[x,y]
    end
    edges
end

"""
    grad(w::Nodes{Primal}) --> Edges{Primal}

Evaluate the discrete gradient of primal nodal data `w`. Can also perform this
operation by creating an object of Grad type and applying it with `*`.

# Example

```jldoctest
julia> w = Nodes(Primal,(8,6));

julia> w[3,4] = 1.0;

julia> G = Grad();

julia> G*w
Edges{Primal,8,6,Float64} data
u (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
v (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  -1.0  0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0   0.0  0.0  0.0  0.0  0.0
```
"""
function grad(p::Nodes{Primal, NX, NY,T}) where {NX, NY,T}
  grad!(Edges(Primal,p),p)
end

"""
    grad!(q::Edges{Dual},w::Nodes{Dual})

Evaluate the discrete gradient of dual nodal data `w` and return it as dual
edge data `q`.
"""
function grad!(edges::Edges{Dual, NX, NY},
                   p::Nodes{Dual, NX, NY}) where {NX, NY}

    @inbounds for y in 1:NY, x in 1:NX-1
        edges.u[x,y] = - p[x,y] + p[x+1,y]
    end
    @inbounds for y in 1:NY-1, x in 1:NX
        edges.v[x,y] = - p[x,y] + p[x,y+1]
    end
    edges
end

"""
    grad(w::Nodes{Dual}) --> Edges{Dual}

Evaluate the discrete gradient of dual nodal data `w`. Can also perform this
operation by creating an object of Grad type and applying it with `*`.
"""
function grad(p::Nodes{Dual, NX, NY}) where {NX, NY}
  grad!(Edges(Dual,p),p)
end

"""
    grad!(d::EdgeGradient{Primal,Dual},q::Edges{Primal})

Evaluate the discrete gradient of primal edge data `q` and return it as edge
gradient data `d`, where the diagonal entries of the gradient lie on primal
nodes and the off-diagonal entries lie at dual nodes.
"""
function grad!(d::EdgeGradient{Primal, Dual, NX, NY},
                     edges::Edges{Primal, NX, NY}) where {NX, NY}

    @inbounds for y in 1:NY-1, x in 1:NX-1
        d.dudx[x,y] = - edges.u[x,y] + edges.u[x+1,y]
        d.dvdy[x,y] = - edges.v[x,y] + edges.v[x,y+1]
    end
    @inbounds for y in 2:NY-1, x in 2:NX-1
        d.dudy[x,y] = - edges.u[x,y-1] + edges.u[x,y]
        d.dvdx[x,y] = - edges.v[x-1,y] + edges.v[x,y]
    end
    d
end

"""
    grad!(d::EdgeGradient{Dual,Primal},q::Edges{Dual})

Evaluate the discrete gradient of dual edge data `q` and return it as edge
gradient data `d`, where the diagonal entries of the gradient lie on dual
nodes and the off-diagonal entries lie at primal nodes.
"""
function grad!(d::EdgeGradient{Dual, Primal, NX, NY},
                     edges::Edges{Dual, NX, NY}) where {NX, NY}

    @inbounds for y in 2:NY-1, x in 2:NX-1
        d.dudx[x,y] = - edges.u[x-1,y] + edges.u[x,y]
        d.dvdy[x,y] = - edges.v[x,y-1] + edges.v[x,y]
    end
    @inbounds for y in 1:NY-1, x in 1:NX-1
        d.dudy[x,y] = - edges.u[x,y] + edges.u[x,y+1]
        d.dvdx[x,y] = - edges.v[x,y] + edges.v[x+1,y]
    end
    d
end

"""
    grad(q::Edges{Primal/Dual}) --> EdgeGradient{Dual/Primal,Primal/Dual}

Evaluate the discrete gradient of primal or dual edge data `q`. Can also perform this
operation by creating an object of Grad type and applying it with `*`.
"""
function grad(edges::Edges{C, NX, NY}) where {C<:CellType,NX,NY}
  grad!(EdgeGradient(C,edges),edges)
end

struct Grad end

(*)(::Grad,w::Union{Nodes{Primal, NX, NY},Nodes{Dual, NX, NY},Edges{Dual,NX,NY},Edges{Primal,NX,NY}}) where {NX,NY} = grad(w)

include("differencing1d.jl")


text:
module FearlessPoker
using FearlessSQL, StatsBase, Memoize

# types to export
export Hand

export RANKS, SUITS, DECK, STARTING_HANDS

# functions to export
export evaluate, equity

# package code goes here
# source all the separate files
srcdir = dirname(@__FILE__)
srcfiles = setdiff(readdir(srcdir), [basename(@__FILE__)])
map(include, srcfiles)

end # module


text:
"""
Matplotlib's style settings [1]. Note that directly editing `mpl.rcParams` has no effect
[2]. Editing this object however does work.

- [1] https://matplotlib.org/stable/tutorials/introductory/customizing.html#the-default-matplotlibrc-file
- [2] https://github.com/JuliaPy/PyPlot.jl#modifying-matplotlibrcparams
"""
rcParams = nothing

"""
A copy of the initial `mpl.rcParams`. Note that we do not use `mpl.rcParamsDefault` or
`mpl.rcParamsOrig`, as these are different to what's actually used by default (e.g. in a
Jupyter notebook).
"""
rcParams_original = nothing

function __init__()
    global rcParams_original = copy(mpl.rcParams)
    global rcParams = PyPlot.PyDict(mpl."rcParams")  # String quotes prevent conversion from
                                                     # Python to Julia dict.
    set_mpl_style!(sciplotlib_style)
end


text:
# This does MCMC, either using raw statistic, or using NN transform,
# depending on the argument usenn
using Flux, Econometrics, LinearAlgebra, Statistics, DelimitedFiles
using BSON:@load
include("lib/transform.jl")
include("lib/lnL.jl")
function MCMC(m, usenn, info)
    # get the trained net
    @load "best.bson" model
    S = nSimulationDraws # number of simulations
    if usenn
        m = transform(m', info)
        m = Float64.(model(m'))
        θinit = (ub+lb)./2.0 # prior mean as initial θ
        lnL = θ -> LL(θ, m, S, model, info)
    else
        θinit = (ub+lb)./2.0 # prior mean as initial θ
        lnL = θ -> LL(θ, m, S)
    end
    # use a rapid SAMIN to get good initialization values for chain
    obj = θ -> -1.0*lnL(θ)
    θinit, junk, junk, junk = samin(obj, θinit, lb, ub; coverage_ok=0, maxevals=100000, verbosity = 0, rt = 0.5)
    Prior = θ -> prior(θ, lb, ub) # uniform, doesn't matter
    # define things for MCMC
    verbosity = false
    ChainLength = 200
    tuning = [0.1, 0.1, 0.1] # fix this somehow
    Proposal = θ -> proposal1(θ, tuning, lb, ub)
    chain = mcmc(θinit, ChainLength, burnin, Prior, lnL, Proposal, verbosity)
    # now use a MVN random walk proposal with updates of covariance and longer chain
    # on final loop
    Σ = NeweyWest(chain[:,1:3])
    tuning = 1.0
    MC_loops = 4 # one round to adjust covariance, then the final
    for j = 1:MC_loops
        P = try
            P = (cholesky(Σ)).U
        catch
            P = diagm(diag(Σ))
        end
        Proposal = θ -> proposal2(θ,tuning*P, lb, ub)
        if j == MC_loops
            ChainLength = 1000
        end
        θinit = mean(chain[:,1:3],dims=1)[:]
        chain = mcmc(θinit, ChainLength, 0, Prior, lnL, Proposal, verbosity)
        if j < MC_loops
            accept = mean(chain[:,end])
            if accept > 0.35
                tuning *= 2.0
            elseif accept < 0.25
                tuning *= 0.25
            end
            Σ = 0.5*Σ + 0.5*NeweyWest(chain[:,1:3])
        end
    end
    # plain MCMC fit
    chain = chain[:,1:3]
    return chain
end


text:
using Documenter
using IntervalArithmetic, Turmeric, NumberIntervals

makedocs(
    modules = [Turmeric],
    doctest = true,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Turmeric.jl",
    authors = "David P. Sanders and Luis Benet, and Josua Grawitter",
    pages = Any[
        "Home" => "index.md",
        "`roots` interface" => "roots.md",
        "Internals" => "internals.md",
        "Bibliography" => "biblio.md",
        "API" => "api.md"
        ]
    )

deploydocs(
    repo = "github.com/gwater/Turmeric.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
export convert_seis, convert_seis!

@doc """
    convert_seis!(S[, chans=CC, units_out=UU, v=V])
    convert_seis(S, chans=CC, units_out=UU, v=V])
    convert_seis!(C[, units_out=UU, v=V])
    convert_seis(CC, units_out=UU, v=V)

Convert all seismic data channels in `S` to velocity seismograms, differentiating or integrating as needed.

### Keywords
* `units_out=UU` specifies output units.
  + Default: "m/s".
  + Allowed: "m", "m/s", or "m/s2". (SeisIO uses Unicode (UTF-8) UCUM units.)
* `v=V` sets verbosity.
* `chans=CC` restricts seismogram conversion to seismic data channels within `CC`.
  + `chans` can be an Integer, UnitRange, or Array{Int64,1}.
  + By default, all seismic data channels in `S` are converted (if needed).
  + This does not allow `convert_seis!` to work on non-seismic data.

!!! warning

    `convert_seis!` becomes less reversible as seismograms lengthen, particularly at Float32 precision.

### References
[^1] Neumaier, A. (1974). "Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen" [Rounding Error Analysis of Some Methods for Summing Finite Sums]. Zeitschrift für Angewandte Mathematik und Mechanik (in German). 54 (1): 39–51. doi:10.1002/zamm.19740540106.

""" convert_seis!
function convert_seis!(S::GphysData;
  chans::ChanSpec=Int64[],
  units_out::String="m/s",
  v::Integer=KW.v)

  # sanity check
  if (units_out in ("m", "m/s", "m/s2")) == false
    error("units_out must be in (\"m\", \"m/s\", \"m/s2\")!")
  end

  # get seismic data channels
  chans = mkchans(chans, S, keepirr=false)
  filt_seis_chans!(chans, S)

  # now loop over all seismic channels
  for i in chans
    if v > 2
      println(stdout, "Begin channel ", i)
    end

    # get units
    units_in = lowercase(S.units[i])
    units_in == units_out && continue

    # differentiate or integrate
    if v > 0
      println(stdout, "Converting channel ", i, ": ", units_in, " => ", units_out)
    end

    T = eltype(S.x[i])
    fs = T(S.fs[i])
    δ = one(T)/fs
    dt = round(Int64, 1.0/(S.fs[i]*μs))
    d2 = div(dt, 2)
    t = S.t[i]
    Nt = size(t, 1) - 1

    if units_out == "m/s"

      # differentiate from m to m/s
      if units_in == "m"
        diff_x!(S.x[i], t[:,1], fs)
        for j = 1:Nt
          t[j,2] -= d2
        end

      # integrate from m/s2 to m/s
      else
        int_x!(S.x[i], t[:,1], δ)
        for j = 1:Nt
          t[j,2] -= d2
        end
      end

    elseif units_out == "m"
      int_x!(S.x[i], t[:,1], δ)
      for j = 1:Nt
        t[j,2] -= d2
      end
      if units_in == "m/s2"
        int_x!(S.x[i], t[:,1], δ)
        for j = 1:Nt
          t[j,2] -= d2
        end
      end
    else # units == "m/s2"
      diff_x!(S.x[i], t[:,1], fs)
      for j = 1:Nt
        t[j,2] += d2
      end
      if units_in == "m"
        diff_x!(S.x[i], t[:,1], fs)
        for j = 1:Nt
          t[j,2] += d2
        end
      end
    end

    # change units
    S.units[i] = units_out

    # log processing
    proc_note!(S, i, string("convert_seis!(S, chans=", i,
                            ", units_out = ", units_out, ")"),
                     string("converted units from ", units_in, " to ", units_out))
    if v > 2
      println(stdout, "Done channel ", i)
    end
  end

  return nothing
end

function convert_seis!(C::GphysChannel;
  units_out::String="m/s",
  v::Integer=KW.v)

  # sanity check
  if (units_out in ("m", "m/s", "m/s2")) == false || (C.units in ("m", "m/s", "m/s2")) == false
    error("units must be in (\"m\", \"m/s\", \"m/s2\")!")
  end

  units_in = lowercase(C.units)
  units_in == units_out && return nothing

  # differentiate or integrate
  if v > 0
    println(stdout, "Converting ", units_in, " => ", units_out)
  end

  T = eltype(C.x)
  fs = T(C.fs)
  δ = one(T)/fs
  dt = round(Int64, 1.0/(C.fs*μs))
  d2 = div(dt, 2)
  t = C.t
  Nt = size(t, 1) - 1

  if units_out == "m/s"

    # differentiate from m to m/s
    if units_in == "m"
      diff_x!(C.x, t[:,1], fs)
      for j = 1:Nt
        t[j,2] -= d2
      end

    # integrate from m/s2 to m/s
    else
      int_x!(C.x, t[:,1], δ)
      for j = 1:Nt
        t[j,2] -= d2
      end
    end

  elseif units_out == "m"
    int_x!(C.x, t[:,1], δ)
    for j = 1:Nt
      t[j,2] -= d2
    end
    if units_in == "m/s2"
      int_x!(C.x, t[:,1], δ)
      for j = 1:Nt
        t[j,2] -= d2
      end
    end
  else # units == "m/s2"
    diff_x!(C.x, t[:,1], fs)
    for j = 1:Nt
      t[j,2] += d2
    end
    if units_in == "m"
      diff_x!(C.x, t[:,1], fs)
      for j = 1:Nt
        t[j,2] += d2
      end
    end
  end

  # change units
  C.units = units_out

  # log processing
  proc_note!(C, string("convert_seis!(C, units_out = ", units_out, ")"),
                string("converted units from ", units_in, " to ", units_out))
  return nothing
end

@doc (@doc convert_seis!)
function convert_seis(S::GphysData;
  chans::ChanSpec=Int64[],
  units_out::String="m/s",
  v::Integer=KW.v)

  U = deepcopy(S)
  convert_seis!(U, chans=chans, units_out=units_out, v=v)
  return U
end

function convert_seis(C::GphysChannel;
  units_out::String="m/s",
  v::Integer=KW.v)

  U = deepcopy(C)
  convert_seis!(U, units_out=units_out, v=v)
  return U
end


text:
module FastFloat16s

    export FastFloat16

    include("fast_float16.jl")

end


text:
## Plot results from email experiment
using MAT
mat = matread("../data/emailEUcore.mat")
label = vec(round.(Int64,mat["labels"]))

mat3  = matread("../data/EmailTriCore.mat")
coreinds = mat3["core_inds"]
label = label[coreinds]

mat = matread("Output/Louvain_Email_clusterings.mat")
Lams = mat["Lams"]

## Get runtimes and clusterings
# Get clusters from running triangle motif clustering...
triC = mat["taC"]
tt = mat["totaltimes"]
TriRun = tt[:,2]

# And standard Louvain clustering with resolution parameter lambda
SimC = mat["SimC"]
SimRun = tt[:,1]

# Results from running other algorithms
grac = matread("Output/Graclus_clusterings_2to340.mat")
GracC = grac["GracCtri"]
GracRun = grac["times"][:,2]

met = matread("Output/Metis_2to340.mat")
MetC = met["MetCtri"]
MetRun = met["times"][:,2]

hmet = matread("Output/hmetis_2to340.mat")
hMetC = hmet["hmetC"]
hMetRun = hmet["times"][:,1]

spec = matread("Output/RecursiveSpectral_Email.mat")
specC = spec["SpecTri"]
SpecRun = spec["times"][:,2]

# Check ARI scores
numpoints = length(Lams)
SimARI = zeros(numpoints)
TriARI = zeros(numpoints)
SimSize = zeros(numpoints)
TriSize = zeros(numpoints)
GracARI= zeros(numpoints)
GracSize = zeros(numpoints)
MetARI= zeros(numpoints)
MetSize = zeros(numpoints)
hMetARI= zeros(numpoints)
hMetSize = zeros(numpoints)
SpecARI= zeros(numpoints)
SpecSize = zeros(numpoints)

include("../src/Graph_Louvain.jl")
for i = 1:length(Lams)
    lam = Lams[i]
    cSim = round.(Int64,SimC[:,i])
    SimARI[i] = ari(cSim,label)
    SimSize[i] = maximum(cSim)

    cTri = round.(Int64,triC[:,i])
    TriARI[i] = ari(cTri,label)
    TriSize[i] = maximum(cTri)

    cGrac = round.(Int64,GracC[:,i])
    GracARI[i] = ari(cGrac,label)
    GracSize[i] = maximum(cGrac)

    cMet = round.(Int64,MetC[:,i])
    MetARI[i] = ari(cMet,label)
    MetSize[i] = maximum(cMet)

    chMet = round.(Int64,hMetC[:,i])
    hMetARI[i] = ari(chMet,label)
    hMetSize[i] = maximum(chMet)

    cSpec = round.(Int64,specC[:,i])
    SpecARI[i] = ari(cSpec,label)
    SpecSize[i] = maximum(cSpec)

end

# After checking, rec-spec is the only alg that
# produces cluster sizes in different orders.
# We adjust for clearest comparison in plots
p = sortperm(SpecSize)
SpecSize = SpecSize[p]
SpecARI = SpecARI[p]
SpecRun = SpecRun[p]

## Make Plots
using Plots

lw = 3
msw = 0
s1 = 300
s2 = 250
ms = 5
lw = 2

## Runtime vs. ARI
p=scatter(TriRun, TriARI,color=:green,linewidth= lw,label = "HypLam",markerstrokecolor = :green,
markerstrokewidth = msw, markersize = ms, markershape = :circle, legend = false)

scatter!(SimRun, SimARI,color=:black,linewidth= lw,grid = false, xscale = :log10,
xaxis = "Runtime", yaxis= "ARI Scores", label = "LamLouv",markerstrokecolor = :black,
markerstrokewidth = msw, markersize = ms, markershape = :circle)

scatter!(GracRun, GracARI,color=:blue,linewidth= lw, label = "Graclus",markerstrokecolor = :blue,
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

scatter!(MetRun, MetARI,color=:orange,linewidth= lw, label = "Metis",markerstrokecolor = :orange,
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

scatter!(hMetRun, hMetARI,color=:purple,linewidth= lw, label = "hMetis",markerstrokecolor = :purple,
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

scatter!(SpecRun, SpecARI,color=:brown,linewidth= lw, label = "Spectral",markerstrokecolor = :brown,
markerstrokewidth = msw, markersize = ms,markershape = :circle, size = (s1,s2),
foreground_color_legend = nothing,background_color_legend=nothing)

##
ms = 10
scatter!([TriRun[findmax(TriARI)[2]]], [TriARI[findmax(TriARI)[2]]],color=:green,linewidth= lw,label = "HypLam",
markerstrokewidth = msw, markersize = ms, markershape = :circle,markerstrokecolor = :green)

scatter!([SimRun[findmax(SimARI)[2]]], [SimARI[findmax(SimARI)[2]]],color=:black,linewidth= lw,grid = false, xscale = :log10,
xaxis = "Runtime", yaxis= "ARI Scores", label = "LamLouv",markerstrokecolor = :black,
markerstrokewidth = msw, markersize = ms, markershape = :circle)

scatter!([GracRun[findmax(GracARI)[2]]], [GracARI[findmax(GracARI)[2]]],color=:blue,linewidth= lw, label = "Graclus",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2),markerstrokecolor = :blue)

scatter!([MetRun[findmax(MetARI)[2]]], [MetARI[findmax(MetARI)[2]]],color=:orange,linewidth= lw, label = "Metis",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2),markerstrokecolor = :orange)

scatter!([hMetRun[findmax(hMetARI)[2]]], [hMetARI[findmax(hMetARI)[2]]],color=:purple,linewidth= lw, label = "hMetis",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2),markerstrokecolor=:purple)

scatter!([SpecRun[findmax(SpecARI)[2]]], [SpecARI[findmax(SpecARI)[2]]],color=:brown,linewidth= lw, label = "Spectral",
markerstrokewidth = msw, markersize = ms,markershape = :circle, size = (s1,s2),markerstrokecolor=:brown)

savefig(p,"Plots/Email-tradeoff.pdf")

## ARI vs. ClusterNumbers
ms = 5
rightmax = minimum( [maximum(SimSize), maximum(TriSize)])
p=plot(TriSize, TriARI,color=:green,markerstrokecolor = :green,linewidth= lw,label = "HypLam",
markerstrokewidth = msw, markersize = ms, markershape = :circle, legend = :topright)

plot!(SimSize, SimARI,color=:black,markerstrokecolor = :black,linewidth= lw,
xlim = [0, rightmax],grid = false,
xaxis = "Number of Clusters", yaxis= "ARI Scores", label = "LamLouv",
markerstrokewidth = msw, markersize = ms, markershape = :circle)

plot!(GracSize, GracARI,color=:blue,markerstrokecolor = :blue,linewidth= lw, label = "Graclus",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

plot!(MetSize, MetARI,color=:orange,markerstrokecolor = :orange,
linewidth= lw, label = "Metis",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

plot!(hMetSize, hMetARI,color=:purple,markerstrokecolor = :purple,
linewidth= lw, label = "hMetis",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

plot!(SpecSize, SpecARI,color=:brown,markerstrokecolor = :brown,
linewidth= lw, label = "Spectral",
markerstrokewidth = msw, markersize = ms,markershape = :circle, size = (s1,s2),
foreground_color_legend = nothing,background_color_legend=nothing)

savefig(p,"Plots/Email-ari.pdf")

## Runtime vs cluster-number

p=plot(TriSize, TriRun,color=:green,markerstrokecolor = :green,
linewidth= lw,label = "HyperLam",yscale = :log10,
markerstrokewidth = msw, markersize = ms, markershape = :circle, legend = :topright)

plot!(SimSize, SimRun,color=:black,linewidth= lw,
xlim = [0, rightmax],grid = false,markerstrokecolor = :black,
xaxis = "Number of Clusters", yaxis= "Runtime (seconds)", label = "LamLouv",
markerstrokewidth = msw, markersize = ms, markershape = :circle)

plot!(GracSize, GracRun,color=:blue,markerstrokecolor = :blue,
linewidth= lw, label = "Graclus",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

plot!(MetSize, MetRun,color=:orange,markerstrokecolor = :orange,
linewidth= lw, label = "Metis",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

plot!(hMetSize, hMetRun,color=:purple,markerstrokecolor = :purple,
linewidth= lw, label = "hMetis", legend = false,
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

plot!(SpecSize, SpecRun,color=:brown,markerstrokecolor = :brown,
linewidth= lw, label = "Spectral",
markerstrokewidth = msw, markersize = ms, markershape = :circle, size = (s1,s2))

savefig(p,"Plots/Email-runtime.pdf")

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
library( Matrix )

str_path <- "C:\\CBIG\\ZD\\"

myd <- read.csv( paste0( str_path, "logistic_regression_ex2.csv" ) )

x <- as.matrix( myd[, -c( 1 )] )
y <- myd[, c(1)]

## make sure its rank = # columns
ncol( x )
rankMatrix( x )

myglm <- glm( y ~ -1 + x, family = binomial( link = logit ) )

myglm_coeff <- summary( myglm )$coefficients

rownames( myglm_coeff ) <- c( paste0( "beta", seq( 1, 3 ) ), paste0( "creative_break_", seq( 1, 4 ) ) )

print( myglm_coeff )

## using optim to reproduce the estimates

Loglikelihood <- function( beta )
{
	xbeta <- x %*% beta
	Li <- -1 * y * xbeta + log( 1 + exp( xbeta ))
	sum( Li )
}

beta_start_value <- as.matrix( rep( 0.0, ncol( x ) ), ncol(x), 1 )

myoptim <- optim( beta_start_value, fn = Loglikelihood, gr = NULL, method = "BFGS", hessian = TRUE )

var_cov_m <- solve( myoptim$hessian )

beta_se <- sqrt( diag( var_cov_m ) )

zvalue <- myoptim$par / beta_se

pvalue <- 2 * ( 1 - pnorm(abs( zvalue ) ) )

coeff_outputs <- cbind( myoptim$par, beta_se, zvalue, round( pvalue, 6 ) )

colnames( coeff_outputs ) <- c( "Estimates", "Std. Error", "z value", "Pr(>|z|)" )
rownames( coeff_outputs ) <- c( paste0( "beta", seq( 1, 3 ) ), paste0( "creative_break_", seq( 1, 4 ) ) )
print( coeff_outputs )

print( summary( myglm )$coefficients )

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

## this is the residual deviance
2 * myoptim$value

########## since colmn 4 to column 7 are indicators, assuming they are the compaign breaks :)
##### so we want to add penalty function to make the estimates >= 0, using the penalty that Zhiyuan recommended

## assuming the penalty function is :

a <- c( 1000, 9, 60 )

Penalized_Loglikelihood <- function( beta )
{
	xbeta <- x %*% beta
	Li <- -1 * y * xbeta + log( 1 + exp( xbeta ))
	## constraining beta[3] to beta[7]
	myPenValue <- 0
	for ( i in 4 : 7 )
	{
		myPenValue <- myPenValue + a[1] * exp( -1.0 * a[2] * exp( a[3] * beta[ i ] ) )
	}
	sum( Li ) + myPenValue
}

beta_start_value <- as.matrix( rep( 0.0, ncol( x ) ), ncol(x), 1 )

myoptim_Penalized <- optim( beta_start_value, fn = Penalized_Loglikelihood, gr = NULL, method = "BFGS", hessian = TRUE)

myoptim_Penalized

var_cov_m <- solve( myoptim_Penalized$hessian )

beta_se <- sqrt( diag( var_cov_m ) )

zvalue <- myoptim_Penalized$par / beta_se

pvalue <- 2 * ( 1 - pnorm(abs( zvalue ) ) )

coeff_outputs <- cbind( myoptim_Penalized$par, beta_se, zvalue, round( pvalue, 6 ) )

colnames( coeff_outputs ) <- c( "Estimates", "Std. Error", "z value", "Pr(>|z|)" )
rownames( coeff_outputs ) <- c( paste0( "beta", seq( 1, 3 ) ), paste0( "creative_break_", seq( 1, 4 ) ) )

print( coeff_outputs )

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
# Reactors
#
# On every MessageCreate event, we can apply a reaction to the message.
# The code below can be easily extends by defining a subtype of
# `AbstractReactor` and implement the respective interface methods.

# ----------------------------------------------------------------------
# Interface
# ----------------------------------------------------------------------

abstract type AbstractReactor end

function reactions(r::AbstractReactor, m::Message)::Vector{Char} end

const NO_REACTION = Char[]

# ----------------------------------------------------------------------
# Reactors
# ----------------------------------------------------------------------

const REACT_WORDS = Dict(
    :happy => [
        "happy",
        "nice",
        "great",
        "awesome",
        "cheers",
        "yay",
        "congrat",
        "congrats",
        "congratulations",
        "it helped",
        "it helps",
        "appreciate",
        "appreciated",
        "noice",
        "thank",
        "thanks",
    ],
    :disappointed => ["disappointed", "unhappy", "sad", "aw shucks", "yeow"],
    :excited => [
        "excited",
        "fantastic",
        "fabulous",
        "wonderful",
        "looking forward to",
        "love",
        "learned",
        "saved me",
        "beautiful",
    ],
    :goodbye => ["cya", "bye", "goodbye", "ciao", "adios", "brb"],
    :dog => ["dog", "dogs", "doggie", "shiba", "corgi", "chihuahua", "retriever"],
    :cat => ["cat", "cats", "feline", "kitten", "kittens"],
    :snake => ["snake", "snakes", "rattle", "python", "pythons"],
    :crab => ["crab", "crabs", "rust"],
)

function contains_any(s::AbstractString, words::AbstractVector{String})
    is_thing(x) = x !== nothing
    regexes = [Regex("\\b" * w * "\\b") for w in words]
    matches = match.(regexes, lowercase(s))
    return any(is_thing(x) for x in matches)
end

struct HappyReactor <: AbstractReactor end

function reactions(::HappyReactor, m::Message)
    if contains_any(m.content, REACT_WORDS[:happy]) &&
       !contains_any(m.content, REACT_WORDS[:disappointed])
        return ['😄']
    end
    return NO_REACTION
end

struct DisappointedReactor <: AbstractReactor end

function reactions(::DisappointedReactor, m::Message)
    if contains_any(m.content, REACT_WORDS[:disappointed]) &&
       !contains_any(m.content, REACT_WORDS[:happy]) &&
       !contains_any(m.content, REACT_WORDS[:excited])
        return ['😞']
    end
    return NO_REACTION
end

struct ExcitedReactor <: AbstractReactor end

function reactions(::ExcitedReactor, m::Message)
    if contains_any(m.content, REACT_WORDS[:excited]) &&
       !contains_any(m.content, REACT_WORDS[:disappointed])
        return ['🤩']
    end
    return NO_REACTION
end

struct GoodbyeReactor <: AbstractReactor end

function reactions(::GoodbyeReactor, m::Message)
    if contains_any(m.content, REACT_WORDS[:goodbye])
        return ['👋']
    end
    return NO_REACTION
end

struct AnimalReactor <: AbstractReactor end

function reactions(::AnimalReactor, m::Message)
    result = Char[]
    if contains_any(m.content, REACT_WORDS[:dog])
        push!(result, '🐕')
    end
    if contains_any(m.content, REACT_WORDS[:cat])
        push!(result, '🐈')
    end
    if contains_any(m.content, REACT_WORDS[:snake])
        push!(result, '🐍')
    end
    if contains_any(m.content, REACT_WORDS[:crab])
        push!(result, '🦀')
    end
    return result
end

# ----------------------------------------------------------------------
# Main logic
# ----------------------------------------------------------------------

const REACTORS = AbstractReactor[
    HappyReactor(),
    DisappointedReactor(),
    ExcitedReactor(),
    GoodbyeReactor(),
    AnimalReactor(),
]

function handler(c::Client, e::MessageCreate, ::Val{:reaction})
    # @info "react_handler called"
    # @info "Received message" e.message.channel_id e.message.id e.message.content
    username = e.message.author.username
    discriminator = e.message.author.discriminator
    !get_opt!(username, discriminator, :reaction) && return nothing
    for reactor in REACTORS
        rs = reactions(reactor, e.message)
        foreach(rs) do emoji
            create(c, Reaction, e.message, emoji)
        end
    end
    return nothing
end


text:
@testset "quadrature" begin

end


text:
## TYPES AND CONSTRUCTOR

const ERR_SPECIFY_MODEL = ArgumentError(
"You need to specify `model=...`, unless `tuning=Explicit()`. ")
const ERR_SPECIFY_RANGE = ArgumentError(
    "You need to specify `range=...`, unless `tuning=Explicit` and "*
    "and `models=...` is specified instead. ")
const ERR_SPECIFY_RANGE_OR_MODELS = ArgumentError(
    "No `model` specified. Either specify an explicit iterator "*
    "of MLJ models over which to optimize (using `models=...`) or "*
    "specify a `model`. ")
const ERR_NEED_EXPLICIT = ArgumentError(
    "You have specified an explicit "*
    "iterator `models` of MLJModels and so cannot "*
    "specify any `tuning` strategy except `Explicit`. Either omit the "*
    "`tuning=...` specification, or specify a *single* model using "*
    "`model=...` instead. ")
const ERR_BOTH_DISALLOWED = ArgumentError(
    "You cannot specify both `model` and `models`. ")
const ERR_MODEL_TYPE = ArgumentError(
    "Only `Deterministic` and `Probabilistic` model types supported.")
const ERR_UNINSTANTIATED_MODEL = AssertionError(
    "Type encountered where model instance expected. (Tuning evaluates "*
    "models by mutating clones of the provided instance, as specified "*
    "by `range`.) ")
const INFO_MODEL_IGNORED =
    "`model` being ignored. Using `model=first(range)`. "
const ERR_TOO_MANY_ARGUMENTS =
    ArgumentError("At most one non-keyword argument allowed. ")
warn_double_spec(arg, model) =
    "Using `model=$arg`. Ignoring keyword specification `model=$model`. "

const ProbabilisticTypes = Union{Probabilistic, MLJBase.MLJModelInterface.ProbabilisticDetector}
const DeterministicTypes = Union{Deterministic, MLJBase.MLJModelInterface.DeterministicDetector}

mutable struct DeterministicTunedModel{T,M<:DeterministicTypes} <: MLJBase.Deterministic
    model::M
    tuning::T  # tuning strategy
    resampling # resampling strategy
    measure
    weights::Union{Nothing,AbstractVector{<:Real}}
    class_weights::Union{Nothing,AbstractDict}
    operation
    range
    selection_heuristic
    train_best::Bool
    repeats::Int
    n::Union{Int,Nothing}
    acceleration::AbstractResource
    acceleration_resampling::AbstractResource
    check_measure::Bool
    cache::Bool
end

mutable struct ProbabilisticTunedModel{T,M<:ProbabilisticTypes} <: MLJBase.Probabilistic
    model::M
    tuning::T  # tuning strategy
    resampling # resampling strategy
    measure
    weights::Union{Nothing,AbstractVector{<:Real}}
    class_weights::Union{Nothing,AbstractDict}
    operation
    range
    selection_heuristic
    train_best::Bool
    repeats::Int
    n::Union{Int,Nothing}
    acceleration::AbstractResource
    acceleration_resampling::AbstractResource
    check_measure::Bool
    cache::Bool
end

const EitherTunedModel{T,M} =
    Union{DeterministicTunedModel{T,M},ProbabilisticTunedModel{T,M}}

MLJBase.caches_data_by_default(::Type{<:EitherTunedModel}) = false

"""
    tuned_model = TunedModel(; model=<model to be mutated>,
                             tuning=RandomSearch(),
                             resampling=Holdout(),
                             range=nothing,
                             measure=nothing,
                             n=default_n(tuning, range),
                             operation=nothing,
                             other_options...)

Construct a model wrapper for hyper-parameter optimization of a
supervised learner, specifying the `tuning` strategy and `model` whose
hyper-parameters are to be mutated.

    tuned_model = TunedModel(; models=<models to be compared>,
                             resampling=Holdout(),
                             measure=nothing,
                             n=length(models),
                             operation=nothing,
                             other_options...)

Construct a wrapper for multiple `models`, for selection of an optimal
one (equivalent to specifying `tuning=Explicit()` and `range=models`
above). Elements of the iterator `models` need not have a common type,
but they must all be `Deterministic` or all be `Probabilistic` *and
this is not checked* but inferred from the first element generated.

See below for a complete list of options.

### Training

Calling `fit!(mach)` on a machine `mach=machine(tuned_model, X, y)` or
`mach=machine(tuned_model, X, y, w)` will:

- Instigate a search, over clones of `model`, with the hyperparameter
  mutations specified by `range`, for a model optimizing the specified
  `measure`, using performance evaluations carried out using the
  specified `tuning` strategy and `resampling` strategy. In the case
  `models` is explictly listed, the search is instead over the models
  generated by the iterator `models`.

- Fit an internal machine, based on the optimal model
  `fitted_params(mach).best_model`, wrapping the optimal `model`
  object in *all* the provided data `X`, `y`(, `w`). Calling
  `predict(mach, Xnew)` then returns predictions on `Xnew` of this
  internal machine. The final train can be supressed by setting
  `train_best=false`.

### Search space

The `range` objects supported depend on the `tuning` strategy
specified. Query the `strategy` docstring for details. To optimize
over an explicit list `v` of models of the same type, use
`strategy=Explicit()` and specify `model=v[1]` and `range=v`.

The number of models searched is specified by `n`. If unspecified,
then `MLJTuning.default_n(tuning, range)` is used. When `n` is
increased and `fit!(mach)` called again, the old search history is
re-instated and the search continues where it left off.

### Measures (metrics)

If more than one `measure` is specified, then only the first is
optimized (unless `strategy` is multi-objective) but the performance
against every measure specified will be computed and reported in
`report(mach).best_performance` and other relevant attributes of the
generated report. Options exist to pass per-observation weights or
class weights to measures; see below.

*Important.* If a custom measure, `my_measure` is used, and the
measure is a score, rather than a loss, be sure to check that
`MLJ.orientation(my_measure) == :score` to ensure maximization of the
measure, rather than minimization. Override an incorrect value with
`MLJ.orientation(::typeof(my_measure)) = :score`.

### Accessing the fitted parameters and other training (tuning) outcomes

A Plots.jl plot of performance estimates is returned by `plot(mach)`
or `heatmap(mach)`.

Once a tuning machine `mach` has bee trained as above, then
`fitted_params(mach)` has these keys/values:

key                 | value
--------------------|--------------------------------------------------
`best_model`        | optimal model instance
`best_fitted_params`| learned parameters of the optimal model

The named tuple `report(mach)` includes these keys/values:

key                 | value
--------------------|--------------------------------------------------
`best_model`        | optimal model instance
`best_history_entry`| corresponding entry in the history, including performance estimate
`best_report`       | report generated by fitting the optimal model to all data
`history`           | tuning strategy-specific history of all evaluations

plus other key/value pairs specific to the `tuning` strategy.

### Complete list of key-word options

- `model`: `Supervised` model prototype that is cloned and mutated to
  generate models for evaluation

- `models`: Alternatively, an iterator of MLJ models to be explicitly
  evaluated. These may have varying types.

- `tuning=RandomSearch()`: tuning strategy to be applied (eg, `Grid()`). See
  the [Tuning
  Models](https://alan-turing-institute.github.io/MLJ.jl/dev/tuning_models/#Tuning-Models)
  section of the MLJ manual for a complete list of options.

- `resampling=Holdout()`: resampling strategy (eg, `Holdout()`, `CV()`),
  `StratifiedCV()`) to be applied in performance evaluations

- `measure`: measure or measures to be applied in performance
  evaluations; only the first used in optimization (unless the
  strategy is multi-objective) but all reported to the history

- `weights`: per-observation weights to be passed the measure(s) in performance
  evaluations, where supported. Check support with `supports_weights(measure)`.

- `class_weights`: class weights to be passed the measure(s) in
  performance evaluations, where supported. Check support with
  `supports_class_weights(measure)`.

- `repeats=1`: for generating train/test sets multiple times in
  resampling ("Monte Carlo" resampling); see [`evaluate!`](@ref) for details

- `operation`/`operations` - One of
  $(MLJBase.PREDICT_OPERATIONS_STRING), or a vector of these of the
  same length as `measure`/`measures`. Automatically inferred if left
  unspecified.

- `range`: range object; tuning strategy documentation describes
  supported types

- `selection_heuristic`: the rule determining how the best model is
  decided. According to the default heuristic,
  `NaiveSelection()`, `measure` (or the first
  element of `measure`) is evaluated for each resample and these
  per-fold measurements are aggregrated. The model with the lowest
  (resp. highest) aggregate is chosen if the measure is a `:loss`
  (resp. a `:score`).

- `n`: number of iterations (ie, models to be evaluated); set by
  tuning strategy if left unspecified

- `train_best=true`: whether to train the optimal model

- `acceleration=default_resource()`: mode of parallelization for
  tuning strategies that support this

- `acceleration_resampling=CPU1()`: mode of parallelization for
  resampling

- `check_measure=true`: whether to check `measure` is compatible with the
  specified `model` and `operation`)

- `cache=true`: whether to cache model-specific representations of
  user-suplied data; set to `false` to conserve memory. Speed gains
  likely limited to the case `resampling isa Holdout`.

"""
function TunedModel(args...; model=nothing,
                    models=nothing,
                    tuning=nothing,
                    resampling=MLJBase.Holdout(),
                    measures=nothing,
                    measure=measures,
                    weights=nothing,
                    class_weights=nothing,
                    operations=nothing,
                    operation=operations,
                    ranges=nothing,
                    range=ranges,
                    selection_heuristic=NaiveSelection(),
                    train_best=true,
                    repeats=1,
                    n=nothing,
                    acceleration=default_resource(),
                    acceleration_resampling=CPU1(),
                    check_measure=true,
                    cache=true)

    # user can specify model as argument instead of kwarg:
    length(args) < 2 || throw(ERR_TOO_MANY_ARGUMENTS)
    if length(args) == 1
        arg = first(args)
        model === nothing ||
            @warn warn_double_spec(arg, model)
        model = arg
    end

    # either `models` is specified and `tuning` is set to `Explicit`,
    # or `models` is unspecified and tuning will fallback to `RandomSearch()`
    # unless it is itself specified:
    if models !== nothing
        if tuning === nothing
            tuning = Explicit()
        elseif !(tuning isa Explicit)
            throw(ERR_NEED_EXPLICIT)
        end
    else
        tuning === nothing && (tuning = RandomSearch())
    end

    # either a `model` is specified or we are in the case
    # `tuning=Explicit()`. In the latter case, either `models` or
    # `range` stores the iterator, but ultimately `range` must do so:
    if model === nothing
        if tuning isa Explicit
            if models === nothing
                range === nothing && throw(ERR_SPECIFY_RANGE_OR_MODELS)
            else
                range = models
            end
        else
            throw(ERR_SPECIFY_MODEL)
        end
    else
        models !== nothing && throw(ERR_BOTH_DISALLOWED)
        range === nothing && throw(ERR_SPECIFY_RANGE)
    end

    # get the model type parameter:
    if tuning isa Explicit
        model !== nothing && @info INFO_MODEL_IGNORED
        model = first(range)
        if model isa Deterministic
            M = Deterministic
        elseif model isa Probabilistic
            M = Probabilistic
        else
            throw(ERR_MODEL_TYPE)
        end
    elseif model isa Type
        throw(ERR_UNINSTANTIATED_MODEL)
    else
        # Model is probably an instantiated model.
        M = typeof(model)
    end

    # get the tuning type parameter:
    T = typeof(tuning)

    args = (
        model,
        tuning,
        resampling,
        measure,
        weights,
        class_weights,
        operation,
        range,
        selection_heuristic,
        train_best,
        repeats,
        n,
        acceleration,
        acceleration_resampling,
        check_measure,
        cache
    )

    if M <: DeterministicTypes
        tuned_model = DeterministicTunedModel{T,M}(args...)
    elseif M <: ProbabilisticTypes
        tuned_model = ProbabilisticTunedModel{T,M}(args...)
    else
        throw(ERR_MODEL_TYPE)
    end

    message = clean!(tuned_model)
    !isempty(message) && @info message

    return tuned_model

end

function MLJBase.clean!(tuned_model::EitherTunedModel)
    message = ""
    if tuned_model.measure === nothing
        tuned_model.measure = default_measure(tuned_model.model)
        if tuned_model.measure === nothing
            error("Unable to deduce a default measure for specified model. "*
                  "You must specify `measure=...`. ")
        else
            message *= "No measure specified. "*
            "Setting measure=$(tuned_model.measure). "
        end
    end

    message *= MLJBase.clean!(tuned_model.tuning)

    if !supports_heuristic(tuned_model.tuning, tuned_model.selection_heuristic)
        message *= "`selection_heuristic=$(tuned_model.selection_heuristic)` "*
        "is not supported by $(tuned_model.tuning). Resetting to "*
        "`NaiveSelectionment()`."
        tuned_model.selection_heuristic = NaiveSelection()
    end

    if (tuned_model.acceleration isa CPUProcesses &&
        tuned_model.acceleration_resampling isa CPUProcesses)
        message *=
        "The combination acceleration=$(tuned_model.acceleration) and"*
        " acceleration_resampling=$(tuned_model.acceleration_resampling) is"*
        "  not generally optimal. You may want to consider setting"*
        " `acceleration = CPUProcesses()` and"*
        " `acceleration_resampling = CPUThreads()`."
    end

    if (tuned_model.acceleration isa CPUThreads &&
        tuned_model.acceleration_resampling isa CPUProcesses)
        message *=
        "The combination acceleration=$(tuned_model.acceleration) and"*
        " acceleration_resampling=$(tuned_model.acceleration_resampling) isn't"*
        " supported. \n Resetting to"*
        " `acceleration = CPUProcesses()` and"*
        " `acceleration_resampling = CPUThreads()`."

        tuned_model.acceleration = CPUProcesses()
        tuned_model.acceleration_resampling = CPUThreads()
    end

    tuned_model.acceleration =
        _process_accel_settings(tuned_model.acceleration)

    return message
end

## FIT AND UPDATE METHODS

# A *metamodel* is either a `Model` instance, `model`, or a tuple
# `(model, s)`, where `s` is extra data associated with `model` that
# the tuning strategy implementation wants available to the `result`
# method for recording in the history.

_first(m::MLJBase.Model) = m
_last(m::MLJBase.Model) = nothing
_first(m::Tuple{Model,Any}) = first(m)
_last(m::Tuple{Model,Any}) = last(m)

# returns a (model, result) pair for the history (called by one of the
# `assemble_events!` methods):
function event!(metamodel,
               resampling_machine,
               verbosity,
               tuning,
               history,
               state)
    model = _first(metamodel)
    metadata = _last(metamodel)
    resampling_machine.model.model = model
    verb = (verbosity >= 2 ? verbosity - 3 : verbosity - 1)
    fit!(resampling_machine, verbosity=verb)
    E = evaluate(resampling_machine)
    entry0 = (model       = model,
              measure     = E.measure,
              measurement = E.measurement,
              per_fold    = E.per_fold,
              metadata    = metadata)
    entry = merge(entry0, extras(tuning, history, state, E))
    if verbosity > 2
        println("hyperparameters: $(params(model))")
    end

    if verbosity > 1
        println("measurement: $(E.measurement[1])")
    end
    return entry
end

function assemble_events!(metamodels,
                         resampling_machine,
                         verbosity,
                         tuning,
                         history,
                         state,
                         acceleration::CPU1)

     n_metamodels = length(metamodels)

     p = Progress(n_metamodels,
         dt = 0,
         desc = "Evaluating over $(n_metamodels) metamodels: ",
         barglyphs = BarGlyphs("[=> ]"),
         barlen = 25,
         color = :yellow)

    verbosity !=1 || update!(p,0)

    entries = map(metamodels) do m
        r = event!(m, resampling_machine, verbosity, tuning, history, state)
        verbosity < 1 || begin
                  p.counter += 1
                  ProgressMeter.updateProgress!(p)
                end
        r
      end

    return entries
end

function assemble_events!(metamodels,
                         resampling_machine,
                         verbosity,
                         tuning,
                         history,
                         state,
                         acceleration::CPUProcesses)

    n_metamodels = length(metamodels)

    entries = @sync begin
        channel = RemoteChannel(()->Channel{Bool}(min(1000, n_metamodels)), 1)
        p = Progress(n_metamodels,
                     dt = 0,
                     desc = "Evaluating over $n_metamodels metamodels: ",
                     barglyphs = BarGlyphs("[=> ]"),
                     barlen = 25,
                     color = :yellow)

        # printing the progress bar
        verbosity < 1 || begin
            update!(p,0)
            @async while take!(channel)
                p.counter +=1
                ProgressMeter.updateProgress!(p)
            end
        end

        ret = @distributed vcat for m in metamodels
            r = event!(m, resampling_machine, verbosity, tuning, history, state)
            verbosity < 1 || begin
                put!(channel, true)
            end
            r
        end
        verbosity < 1 || put!(channel, false)
        ret
    end

    return entries
end

@static if VERSION >= v"1.3.0-DEV.573"
# one machine for each thread; cycle through available threads:
function assemble_events!(metamodels,
                         resampling_machine,
                         verbosity,
                         tuning,
                         history,
                         state,
                         acceleration::CPUThreads)

    if Threads.nthreads() == 1
        return assemble_events!(metamodels,
                         resampling_machine,
                         verbosity,
                         tuning,
                         history,
                         state,
                         CPU1())
   end

    n_metamodels = length(metamodels)
    ntasks = acceleration.settings
    partitions = chunks(1:n_metamodels, ntasks)
    #tasks = Vector{Task}(undef, length(partitions))
    entries = Vector(undef, length(partitions))
    p = Progress(n_metamodels,
         dt = 0,
         desc = "Evaluating over $(n_metamodels) metamodels: ",
         barglyphs = BarGlyphs("[=> ]"),
         barlen = 25,
         color = :yellow)
    ch = Channel{Bool}(min(1000, length(partitions)) )

    @sync begin
        # printing the progress bar
        verbosity < 1 || begin
            update!(p,0)
            @async while take!(ch)
                p.counter +=1
                ProgressMeter.updateProgress!(p)
            end
        end
        # One resampling_machine per task
         machs = [resampling_machine,
                 [machine(Resampler(
                     model= resampling_machine.model.model,
                     resampling    = resampling_machine.model.resampling,
                     measure       = resampling_machine.model.measure,
                     weights       = resampling_machine.model.weights,
                     class_weights = resampling_machine.model.class_weights,
                     operation     = resampling_machine.model.operation,
                     check_measure = resampling_machine.model.check_measure,
                     repeats       = resampling_machine.model.repeats,
                     acceleration  = resampling_machine.model.acceleration,
                     cache         = resampling_machine.model.cache),
                          resampling_machine.args...; cache=false) for
                  _ in 2:length(partitions)]...]

        @sync for (i, parts) in enumerate(partitions)
            Threads.@spawn begin
                entries[i] =  map(metamodels[parts]) do m
                    r = event!(m, machs[i],
                              verbosity, tuning, history, state)
                    verbosity < 1 || put!(ch, true)
                    r
                end
            end
        end
        verbosity < 1 || put!(ch, false)
    end
    reduce(vcat, entries)
end

end # of if VERSION ...

# history is intialized to `nothing` because it's type is not known.
_vcat(history, Δhistory) = vcat(history, Δhistory)
_vcat(history::Nothing, Δhistory) = Δhistory
_length(history) = length(history)
_length(::Nothing) = 0

# builds on an existing `history` until the length is `n` or the model
# supply is exhausted (method shared by `fit` and `update`). Returns
# the bigger history. Called by `fit` and `update`.
function build!(history,
               n,
               tuning,
               model,
               model_buffer,
               state,
               verbosity,
               acceleration,
               resampling_machine)
    j = _length(history)
    models_exhausted = false

    # before generating new models be sure to exhaust the model
    # buffer:
    if isready(model_buffer)
        metamodels = [take!(model_buffer),]
        j += 1
        while isready(model_buffer) && j < n
            push!(metamodels, take!(model_buffer))
            j += 1
        end
        Δhistory = assemble_events!(metamodels,
                                    resampling_machine,
                                    verbosity,
                                    tuning,
                                    history,
                                    state,
                                    acceleration)
        history = _vcat(history, Δhistory)
    end

    while j < n && !models_exhausted
        metamodels, state  = MLJTuning.models(tuning,
                                              model,
                                              history,
                                              state,
                                              n - j,
                                              verbosity)
        Δj = _length(metamodels)
        Δj == 0 && (models_exhausted = true)
        shortfall = n - j - Δj
        if models_exhausted && shortfall > 0 && verbosity > -1
            @info "Only $j (of $n) models evaluated.\n"*
            "Model supply exhausted. "
        end
        Δj == 0 && break
        if shortfall < 0 # ie, we have a surplus of models
            # add surplus to buffer:
            for i in (n - j + 1):length(metamodels)
                put!(model_buffer, metamodels[i])
            end
            # and truncate:
            metamodels = metamodels[1:n - j]
        end

        Δhistory = assemble_events!(metamodels,
                                   resampling_machine,
                                   verbosity,
                                   tuning,
                                   history,
                                   state,
                                   acceleration)
        history = _vcat(history, Δhistory)
        j += Δj

    end
    return history, state
end

# given complete history, pick out best model, fit it on all data and
# generate report and cache (meta_state):
function finalize(tuned_model,
                  model_buffer,
                  history,
                  state,
                  verbosity,
                  rm,
                  data...)
    model = tuned_model.model
    tuning = tuned_model.tuning

    user_history = map(history) do entry
        delete(entry, :metadata)
    end

    entry =  best(tuned_model.selection_heuristic, history)
    best_model = entry.model
    best_history_entry = delete(entry, :metadata)
    fitresult = machine(best_model, data...)

    report0 = (best_model         = best_model,
               best_history_entry = best_history_entry,
               history            = user_history)

    if tuned_model.train_best
        fit!(fitresult, verbosity=verbosity - 1)
        report1 = merge(report0, (best_report=MLJBase.report(fitresult),))
    else
        report1 = merge(report0, (best_report=missing,))
    end

    report = merge(report1, tuning_report(tuning, history, state))
    meta_state = (history, deepcopy(tuned_model), model_buffer, state, rm)

    return fitresult, meta_state, report
end

function MLJBase.fit(tuned_model::EitherTunedModel{T,M},
                     verbosity::Integer, data...) where {T,M}
    tuning = tuned_model.tuning
    model = tuned_model.model
    _range = tuned_model.range
    n = tuned_model.n === nothing ?
        default_n(tuning, _range) : tuned_model.n

    verbosity < 1 || @info "Attempting to evaluate $n models."

    acceleration = tuned_model.acceleration

    state = setup(tuning, model, _range, tuned_model.n, verbosity)
    model_buffer = Channel(Inf)

    # instantiate resampler (`model` to be replaced with mutated
    # clones during iteration below):
    resampler = Resampler(model=model,
                          resampling    = deepcopy(tuned_model.resampling),
                          measure       = tuned_model.measure,
                          weights       = tuned_model.weights,
                          class_weights  = tuned_model.class_weights,
                          operation     = tuned_model.operation,
                          check_measure = tuned_model.check_measure,
                          repeats       = tuned_model.repeats,
                          acceleration  = tuned_model.acceleration_resampling,
                          cache         = tuned_model.cache)
    resampling_machine = machine(resampler, data...; cache=false)
    history, state = build!(nothing, n, tuning, model, model_buffer, state,
                           verbosity, acceleration, resampling_machine)

    rm = resampling_machine
    return finalize(tuned_model, model_buffer,
                    history, state, verbosity, rm, data...)

end

function MLJBase.update(tuned_model::EitherTunedModel,
                        verbosity::Integer,
                        old_fitresult, old_meta_state, data...)

    history, old_tuned_model, model_buffer, state, resampling_machine =
        old_meta_state
    acceleration = tuned_model.acceleration

    tuning = tuned_model.tuning
    range = tuned_model.range
    model = tuned_model.model

    # exclamation points are for values actually used rather than
    # stored:
    n! = tuned_model.n === nothing ?
        default_n(tuning, range) : tuned_model.n

    old_n! = old_tuned_model.n === nothing ?
        default_n(tuning, range) : old_tuned_model.n

    if MLJBase.is_same_except(tuned_model, old_tuned_model, :n) &&
        n! >= old_n!

        verbosity < 1 || @info "Attempting to add $(n! - old_n!) models "*
        "to search, bringing total to $n!. "

        history, state = build!(history, n!, tuning, model, model_buffer, state,
                               verbosity, acceleration, resampling_machine)

        rm = resampling_machine
        return finalize(tuned_model, model_buffer,
                        history, state, verbosity, rm, data...)
    else
        return  fit(tuned_model, verbosity, data...)
    end
end

MLJBase.predict(tuned_model::EitherTunedModel, fitresult, Xnew) =
    predict(fitresult, Xnew)

function MLJBase.fitted_params(tuned_model::EitherTunedModel, fitresult)
    if tuned_model.train_best
        return (best_model=fitresult.model,
                best_fitted_params=fitted_params(fitresult))
    else
        return (best_model=fitresult.model,
                best_fitted_params=missing)
    end
end

## SUPPORT FOR MLJ ITERATION API

MLJBase.iteration_parameter(::Type{<:EitherTunedModel}) = :n
MLJBase.supports_training_losses(::Type{<:EitherTunedModel}) = true

function MLJBase.training_losses(tuned_model::EitherTunedModel, _report)
    _losses = MLJTuning.losses(tuned_model.selection_heuristic, _report.history)
    MLJTuning._length(_losses) == 0 && return nothing

    ret = similar(_losses)
    lowest = first(_losses)
    for i in eachindex(_losses)
        current = _losses[i]
        lowest = min(current, lowest)
        ret[i] = lowest
    end
    return ret
end

## METADATA

MLJBase.is_wrapper(::Type{<:EitherTunedModel}) = true
MLJBase.supports_weights(::Type{<:EitherTunedModel{<:Any,M}}) where M =
    MLJBase.supports_weights(M)
MLJBase.supports_class_weights(::Type{<:EitherTunedModel{<:Any,M}}) where M =
    MLJBase.supports_class_weights(M)
MLJBase.load_path(::Type{<:ProbabilisticTunedModel}) =
    "MLJTuning.ProbabilisticTunedModel"
MLJBase.load_path(::Type{<:DeterministicTunedModel}) =
    "MLJTuning.DeterministicTunedModel"
MLJBase.package_name(::Type{<:EitherTunedModel}) = "MLJTuning"
MLJBase.package_uuid(::Type{<:EitherTunedModel}) =
    "03970b2e-30c4-11ea-3135-d1576263f10f"
MLJBase.package_url(::Type{<:EitherTunedModel}) =
    "https://github.com/alan-turing-institute/MLJTuning.jl"
MLJBase.package_license(::Type{<:EitherTunedModel}) = "MIT"
MLJBase.is_pure_julia(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.is_pure_julia(M)
MLJBase.input_scitype(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.input_scitype(M)
MLJBase.target_scitype(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.target_scitype(M)

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
struct EPOCHSimulation{P,B}
    dir::String
    files::Vector{SDFFile{P,B}}
    param::P
end

function read_simulation(dir)
    file_list = joinpath(dir, "normal.visit")
    if isfile(file_list)
        paths = readlines(file_list)
    else
        @debug "No normal.visit found in $dir."
        paths = filter(f->endswith(f, ".sdf"), readdir(dir))
    end

    input_deck = joinpath(dir, "input.deck")
    if isfile(input_deck)
        p = parse_input(input_deck)
    else
        @warn "No input.deck found in $dir."
        p = missing
    end

    files = read_file.(joinpath.((dir,), paths), (Ref(p),))

    EPOCHSimulation(dir, files, p)
end

get_parameter(sim::EPOCHSimulation, p::Symbol) = getindex(sim.param, p)
get_parameter(sim::EPOCHSimulation, p::Symbol, c::Symbol) = getindex(get_parameter(sim, p), c)

# Indexing
Base.getindex(sim::EPOCHSimulation, i::Int) = sim.files[i]
Base.firstindex(sim::EPOCHSimulation) = firstindex(sim.files)
Base.lastindex(sim::EPOCHSimulation) = lastindex(sim.files)

# Iteration
Base.iterate(sim::EPOCHSimulation, state...) = iterate(sim.files, state...)
Base.eltype(::Type{EPOCHSimulation}) = SDFFile
Base.length(sim::EPOCHSimulation) = length(sim.files)

# Statistics
function Statistics.mean(f::Function, sim::EPOCHSimulation; cond=x->true)
    ThreadsX.map(sim.files) do file
        println("Loading $(file.name)")
        qunatity = f(file)
        z = zero(eltype(qunatity))
        (qunatity ./ length(qunatity)) |>
            Filter(cond) |>
            foldxt(+, simd=true, init=z)
    end
end

# FileTrees
# FileTrees._maketree(node::SDFFile) = File(nothing, basename(node.name), node)


text:
@testset "Threading" begin
    @testset "Check number of threads" begin
        @test ccall(Arblib.@libflint(flint_get_num_threads), Cint, ()) == 1

        Arblib.flint_set_num_threads(2)
        @test ccall(Arblib.@libflint(flint_get_num_threads), Cint, ()) == 2

        Arblib.flint_set_num_threads(4)
        @test ccall(Arblib.@libflint(flint_get_num_threads), Cint, ()) == 4

        Arblib.flint_set_num_threads(1)
        @test ccall(Arblib.@libflint(flint_get_num_threads), Cint, ()) == 1
    end

    @testset "Disabled threading" begin
        # This manually unsets the threading flag and checks that we
        # get an error when trying to set the number of threads
        Arblib.__isthreaded[] = false

        @test_throws ErrorException Arblib.flint_set_num_threads(2)

        Arblib.__isthreaded[] = true
    end

    @testset "Matrix multiplication" begin
        A = ArbMatrix(reshape(1:10000, 100, 100))

        Arblib.flint_set_num_threads(1)
        B1 = Arblib.mul_threaded!(similar(A), A, A)

        Arblib.flint_set_num_threads(2)
        B2 = Arblib.mul_threaded!(similar(A), A, A)

        @test B1 == B2 == A * A

        Arblib.flint_set_num_threads(2)
    end
end


text:
#Main definitions
export GaussianProcess, SEGP, MSEGP, query, train!, update!

type GaussianProcess{N}
    m::MeanFunc
    K::CovFunc
    sig_n2::Float64
    Xtrn::Matrix{Float64}
    ftrn::Vector{Float64}

    #These variables are computed during training
    L::Matrix{Float64}
    alpha::Vector{Float64}
    #For debug printing
    printlevel::Int

    function GaussianProcess{N}(m, K, sn2) where N
        new(m, K, sn2, zeros(N,0), zeros(0), zeros(0,0), zeros(0), 0)
    end
end
SEGP(K,sn2,Nd) = GaussianProcess{Nd}(ConstantMean(0.0), K, sn2)
SEGP(K,sn2,X,y) = train!(SEGP(K,sn2,size(X,1)), X, y)

MSEGP(m,K,sn2,Nd) = GaussianProcess{Nd}(m, K, sn2)
MSEGP(m,K,sn2,X,y) = train!(MSEGP(m,K,sn2,size(X,1)), X, y)

Base.show{N}(io::IO, gp::GaussianProcess{N}) =
    print(io, "GP{$N}(Ndata = $(length(gp.ftrn)), K = $(gp.K))")
Base.size{D}(gp::GaussianProcess{D}) = D,length(gp)
Base.length(gp::GaussianProcess) = size(gp.Xtrn, 2)
getx{D}(gp::GaussianProcess{D}, idx::Int) = gp.Xtrn[:,idx]
getf{D}(gp::GaussianProcess{D}, idx::Int) = gp.ftrn[idx]

#Do calculations that depend on training data only
function train!{N}(gp::GaussianProcess{N}, xt, ft)
    gp.Xtrn = copy(xt)
    gp.ftrn = copy(ft)

    fmean = Float64[gp.m(gp.Xtrn[:,i]) for i=1:length(ft)]

    Kt = gp.K(gp.Xtrn, gp.Xtrn)
    Kn = Kt + gp.sig_n2*eye(length(ft)) #n x n
    gp.L = try
        chol(Kn)'
    catch ex
        println("training K info: ", gp.K)
        error("chol failure")
    end
    gp.alpha = gp.L' \ (gp.L \ (gp.ftrn - fmean)) #n x 1
    gp
end

#An incremental update given a single new point. Just as expensize as a full
#batch training.
function update!{N}(gp::GaussianProcess{N}, xn::Vector, fn::Real)
    Xtrn = length(gp) == 0 ?
        reshape(xn,length(xn),1) :
        hcat(gp.Xtrn, xn)
    train!(gp, Xtrn, [gp.ftrn; fn])
end

#Incremental update with multiple new points
function update!{N}(gp::GaussianProcess{N}, Xn::Matrix, fn::Vector)
    Xtrn = hcat(gp.Xtrn, Xn)
    ftrn = vcat(gp.ftrn, fn)
    train!(gp, Xtrn, ftrn)
end

#Evaluate the mean and variance of a single test point
function query{N}(x::Vector, gp::GaussianProcess{N})
    ktrn = vec(gp.K(gp.Xtrn, x)) #n x 1
    mt = dot(ktrn, gp.alpha) + gp.m(x)
    v = gp.L \ ktrn
    vt = gp.K(x,x) - dot(v,v)
    mt, vt
end

#Marginal log-likelihood
function loglikelihood{N}(gp::GaussianProcess{N})
    -0.5*dot(gp.ftrn,gp.alpha) -
        sum(log(diag(gp.L))) -
        length(gp.ftrn)*log(2pi)/2
end

#Condition a set of test points given a trained GP
function condition{N}(x, gp::GaussianProcess{N})
    @assert length(gp.ftrn) != 0 "The process must be trained to condition"
    K_x_xt  = gp.K(x, gp.Xtrn)
    K_xt_x  = K_x_xt'
    K_x_x   = gp.K(x, x)
    Kinv = inv(gp.K(gp.Xtrn, gp.Xtrn))

    mcond = K_x_xt * Kinv * gp.ftrn
    Kcond = K_x_x - K_x_xt * Kinv * K_xt_x
    mcond, Kcond
end

function gradient!{N}(grad::Vector, x::Vector, gp::GaussianProcess{N})
    grad[:] = covgrad(x, gp.K, gp.Xtrn)'*gp.alpha
    grad
end
gradient{N}(x::Vector, gp::GaussianProcess{N}) = gradient!(zeros(N), x, gp)

function checkgrad{N}(x::Vector, gp::GaussianProcess{N}, dx)
    grad = zeros(N)
    for i=1:N
        dv = zeros(N)
        dv[i] = dx
        xt = x + dv
        m1,sig1 = query(x, gp)
        m2,sig2 = query(xt,gp)
        grad[i] = (m2-m1)/dx
    end
    grad
end


text:
#if !@isdefined  Distributions
#  using Distributions
#end

struct MultiUniform <: Distribution
  a::Float64
  b::Float64
  num_param::Integer

  MultiUniform(a::Real, b::Real) = new(a, b, 1)
  MultiUniform(a::Real, b::Real, num_param::Integer) = new(a, b, num_param)
  MultiUniform(num_param::Integer) = new(0.0, 1.0, num_param)
  MultiUniform() = new(0.0, 1.0, 1)
end

mean(d::MultiUniform) = middle(d.a, d.b)
median(d::MultiUniform) = mean(d)

pdf(d::MultiUniform, x::Any) = 1.0 / (d.b - d.a)

function rand(d::MultiUniform)
  result_arr = Vector{Float64}(d.num_param)
  for i in 1:d.num_param
    result_arr[i] = d.a + (d.b - d.a) * rand()
  end
  return result_arr
end


text:
working_dir = nothing;

function set_working_dir(dir::AbstractString)
    if !isdir(dir)
        error("$(dir) not found!!!")
    end
    global working_dir = dir;
end

function get_working_dir()
    if working_dir == nothing
        error("You must first defined the working_dir!!!")
    elseif !isdir(working_dir)
        error("$(working_dir) not found!!!")
    end
    return working_dir;
end


text:
using SwitchTimeOpt
using MathOptInterface
const MOI = MathOptInterface

using Plots
using Ipopt

# Define Solver options
maxiter = 100
maxtime = 1000.0
verbose = 5
tolerance = 1e-06

#Define solver
solver = Ipopt.Optimizer()
MOI.set(solver, MOI.RawOptimizerAttribute("tol"), tolerance)
MOI.set(solver, MOI.RawOptimizerAttribute("print_level"), verbose)
MOI.set(solver, MOI.RawOptimizerAttribute("max_cpu_time"), maxtime)
MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), maxiter)

### Define system parameters
# Size of the state space
nx = 2;

# Number of switches
N = 5;

# Initial state
x0 = [1.0; 1.0]

# Define initial and final time
t0 = 0.0
tf = 1.0

### Define system dynamics
A = zeros(nx, nx, N+1)

A[:, :, 1] = [-1 0;
              1  2]
A[:, :, 2] = [1 1;
              1 -2]
for i = 3:N+1
  A[:, :, i] = A[:, :, mod(i+1, 2)+1]
end

### Define and solve switching time optimization problem
# m = stoproblem(x0, A, solver=solver)
#
m = stoproblem(x0, A, t0=t0, tf=tf, ngrid=100)

# Solve problem
solve!(m)

# Obtain values
tauopt = gettau(m)
Jopt = getobjval(m)
soltime = getsoltime(m)

# Simulate linear system
xsim, xpts, Jsim, t = simulate(m)
println("Objective: ", Jsim)

uvec = zeros(2,N+1)
for i=1:N+1
  uvec[mod(i+1,2)+1,i] = 1
end

plotSolution(m, uvec, "Linear example", "Linear example", true, false)


text:
@testset "Generate PSF method" begin
    # large aperture is delta peak
    out = zeros((5, 5))
    out[1,1] = 1

    @test out ≈ simple_psf((5, 5), 100)
    out[1,1] = 0
    out[3,3] = 1
    @test out ≈ simple_psf((5, 5), 100, shift=true)

    # pinhole aperture
    out = ones((10, 10))
    out ./= sum(out)
    @test out ≈ simple_psf((10, 10), 0.01)

    # normalized
    @test 1 ≈ sum(simple_psf((100, 100), 10))
end

@testset "Generate jinc PSF method" begin
    # large aperture is delta peak
    out = zeros((5, 5))
    out[1,1] = 1

    # large aperture is roughly a delta
    out[1,1] = 0
    out[3,3] = 1
    @test out ≈ jinc_psf((5, 5), 1, 1, shift=true)

    # pinhole aperture
    out = ones((10, 10))
    out ./= sum(out)
    @test out ≈ jinc_psf((10, 10), 1, 1e-12)

    # normalized
    @test 1 ≈ sum(jinc_psf((100, 100), 10f0, 2f0, 1; λ=1f0, shift=true))
    @test 1 ≈ sum(jinc_psf((100, 100), 10f0, 2f0, 1e-20; λ=1f0, shift=true))

    # check consistency with Rayleigh resolution
    d = 1
    f = 1
    r = rayleigh_criterion(f, d)

    psf = jinc_psf((12, 12), 2 * r, d, f)
    @test psf[1, 6] > psf[1, 7] && psf[1, 7] < psf[1, 8]

    psf = jinc_psf((1024, 1024), 2 * r, d, f)
    @test 1 .+ psf[513, 513] ≈ 1
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
import Term: RenderableText, Spacer, vLine, hLine

@testset "\e[31mlayout - spacer" begin
    sizes = [(22, 1), (44, 123), (21, 1), (4334, 232)]
    for (w, h) in sizes
        spacer = Spacer(w, h)
        @test spacer.measure.w == w
        @test spacer.measure.h == h
    end
end

@testset "\e[31mlayout - vLine " begin
    for h in [1, 22, 55, 11]
        line = vLine(h)
        @test length(line.segments) == h
        @test line.measure.h == h
        @test line.height == h
    end

    lines = [(22, "bold"), (55, "red on_green")]
    for (h, style) in lines
        @test vLine(h; style = style).measure.h == h
    end

    for box in (:MINIMAL_DOUBLE_HEAD, :DOUBLE, :ASCII, :DOUBLE_EDGE)
        @test vLine(22; box = box).height == 22
    end

    @test vLine().height == displaysize(stdout)[1]

end

@testset "\e[31mlayout - hLine " begin
    for w in [1, 342, 433, 11, 22]
        line = hLine(w)
        @test line.width == w
        @test length(line.segments) == 1
        @test line.measure.w == w
    end

    for box in (:MINIMAL_DOUBLE_HEAD, :DOUBLE, :ASCII, :DOUBLE_EDGE)
        @test hLine(22; box = box).width == 22
        @test hLine(22, "title"; box = box).width == 22
    end

    for style in ("bold", "red on_green", "blue")
        @test hLine(11; style = style).width == 11
        @test hLine(11, "ttl"; style = style).width == 11
    end

    @test hLine().width == displaysize(stdout)[2]
end

@testset "\e[31mlayout - vstack" begin
    r1 = RenderableText("."^100; width = 25)
    r2 = RenderableText("."^100; width = 50)

    r = r1 / r2
    @test r.measure.w == 50
    @test r.measure.h == 6

    s1 = "."^50
    s2 = ".\n"^25

    r = s1 / s2
    @test r.measure.w == 50
    @test r.measure.h == 27

    # stack other renderables
    h1 = hLine(22)
    h2 = hLine(33)
    @test (h1 / h2).measure.w == 33
    @test (h1 / h2).measure.h == 2

end

@testset "\e[31mlayout - hstack" begin
    r1 = RenderableText("."^100; width = 25)
    r2 = RenderableText("."^100; width = 50)

    r = r1 * r2
    @test r.measure.w == 75
    @test r.measure.h == 4

    # stack other renderables
    h1 = vLine(22)
    h2 = vLine(33)
    @test (h1 * h2).measure.w == 2
    @test (h1 * h2).measure.h == 33

end