text:
module CellFieldsTests

using Test
using BlockArrays
using FillArrays
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.CellData
using Gridap.FESpaces
using Random
using StaticArrays

domain = (0,1,0,1)
cells = (3,3)
model = CartesianDiscreteModel(domain,cells)

trian = Triangulation(model)
trian_N =BoundaryTriangulation(model)
trian_D =BoundaryTriangulation(model,tags="tag_8")
trian_S =SkeletonTriangulation(model)
trian_0 =Triangulation(trian_D,Int[])

ϕ = GenericCellField(get_cell_map(trian),trian,ReferenceDomain())
ϕinv = GenericCellField(lazy_map(inverse_map,get_cell_map(trian)),trian,PhysicalDomain())
xref = get_cell_points(trian)
xphy = CellPoint(get_array(xref),trian,PhysicalDomain())
test_array(ϕ(xref),collect1d(ϕ(xphy)),(i,j)->all(map(≈,i,j)))
test_array(ϕinv(xref),collect1d(ϕinv(xphy)),(i,j)->all(map(≈,i,j)))

x = get_cell_points(trian)
@test DomainStyle(x) == ReferenceDomain()
@test get_array(x) == get_cell_coordinates(trian)
@test get_data(x) == get_cell_ref_coordinates(trian)

px = get_physical_coordinate(trian)
test_array(px(x),collect1d(get_array(x)))

_x = change_domain(x,PhysicalDomain())
@test DomainStyle(_x) == PhysicalDomain()
@test get_array(_x) == get_cell_coordinates(trian)
@test get_data(_x) == get_cell_coordinates(trian)

_x = change_domain(x,ReferenceDomain())
@test DomainStyle(_x) == ReferenceDomain()
@test get_array(x) == get_cell_coordinates(trian)
@test get_data(x) == get_cell_ref_coordinates(trian)

ffun(x) = 2*x[1]
f = CellField(ffun,trian)
fx = f(x)
r = map(xs->ffun.(xs),get_array(x))
r = reshape(r,length(r))
test_array(fx,r,≈)

x_0 = get_cell_points(trian_0)
fx_0 = f(x_0)
test_array(fx_0,collect(fx_0))

n_S = get_normal_vector(trian_S)
x_S = get_cell_points(trian_S)

nf_S = n_S⋅∇(f)

jnf_S = jump(n_S⋅∇(f))
jnfx_S = jnf_S(x_S)
test_array(jnfx_S,0*collect(jnfx_S))

h = CellField(rand(num_cells(trian_S)),trian_S)*jump(∇(f))
hx_S = h(x_S)
test_array(hx_S,collect(hx_S))

h = 3*mean(f)⋅jump(n_S⋅∇(f))
hx_S = h(x_S)
test_array(hx_S,0*collect(hx_S))

aa(f) = 4*f
test_array((aa∘f)(x),4*r,≈)

f1 = f
f2 = 2*f
b(f1,f2) = f1+f2
test_array((b∘(f1,f2))(x),3*r,≈)

f2 = ones(num_cells(trian))
test_array((b∘(f1,f2))(x),map(i->i.+1,r),≈)

v = GenericCellField(get_cell_shapefuns(trian),trian,ReferenceDomain())
vx = v(x)
test_array(vx,collect(vx))

u = GenericCellField(lazy_map(transpose,get_data(v)),v.trian,v.domain_style)
m = v*u
test_array(m(x),collect(m(x)))
m = ∇(v)⋅∇(u)
test_array(m(x),collect(m(x)))

∇vx = ∇(v)(x)
test_array(∇vx,collect(∇vx))

∇fx = ∇(f)(x)
test_array(∇fx,collect(∇fx))

k = VectorValue(1.0,2.0)
∇kfx = ((∇+k)(f))(x)
test_array(∇kfx,collect(∇kfx))

∇kvx = ((∇+k)(v))(x)
test_array(∇kvx,collect(∇kvx))

β(x) = 2*x[1]
α = CellField(x->2*x,trian)
ax = ((∇+k)(β*α))(x)
test_array(ax,collect(ax))

ν = CellField(x->2*x,trian)
ax =((∇-k)⋅ν)(x)
test_array(ax,collect(ax))

ax =((∇-k)×ν)(x)
test_array(ax,collect(ax))

ax =((∇-k)⊗ν)(x)
test_array(ax,collect(ax))

ax =(∇.*ν)(x)
test_array(ax,collect(ax))

ax =(ν.*ν)(x)
test_array(ax,collect(ax))

ax =((∇-k).*ν)(x)
test_array(ax,collect(ax))

ax =(ν⊗(∇-k))(x)
test_array(ax,collect(ax))

σ(x) = diagonal_tensor(VectorValue(1*x[1],2*x[2]))
Fields.gradient(::typeof(σ)) = x-> ThirdOrderTensorValue{2,2,2,Float64}(1,0,0,0,0,0,0,2)
ax = ((∇+k)(σ⋅α))(x)
test_array(ax,collect(ax))

h = Operation(*)(2,f)
hx = h(x)
test_array(hx,2*fx)

a = fill(2,num_cells(trian))
h = Operation(*)(a,f)
hx = h(x)
test_array(hx,2*fx)

fx = evaluate(ffun,x)
test_array(fx,r,≈)

f_N = CellField(ffun,trian_N)
x_N = get_cell_points(trian_N)
fx_N = f_N(x_N)
test_array(fx_N,collect(fx_N))

n_N = get_normal_vector(trian_N)
nx_N = n_N(x_N)
test_array(nx_N,collect(nx_N))

h = f*n_N
hx = h(x_N)
test_array(hx,collect(hx))

gfun(x) = 3*x
g = CellField(gfun,trian)

h = Operation(*)(f,g)
gx = g(x)
hx = h(x)
r = map((i,j)->broadcast(*,i,j),fx,gx)
test_array(hx,r)

h_N = Operation(*)(f_N,g)
gx_N = g(x_N)
r = map((i,j)->broadcast(*,i,j),fx_N,gx_N)
hx_N = h_N(x_N)
test_array(hx_N,r)

g_D = CellField(gfun,trian_D)

cell_h = rand(num_cells(trian))
h = CellField(cell_h,trian)
test_array(h(x),collect(h(x)))
test_array(h(x_N),collect(h(x_N)))

h_N = (2*f_N+g)⋅g
hx_N = h_N(x_N)
test_array(hx_N,collect(hx_N))

# Test function evaluation

# Set reproducible random number seed
Random.seed!(0)
@testset "evaluating functions" for D in 1:3
    xmin = 0
    xmax = 1
    domain = repeat([xmin, xmax], D)
    ncells = 3
    partition = repeat([ncells], D)
    base_model = CartesianDiscreteModel(domain, partition)
    order = 2
    reffe = ReferenceFE(lagrangian, Float64, order)

    for M ∈ [:hypercubes, :simplices]
        model = (M == :hypercubes) ? base_model : simplexify(base_model)

        V = FESpace(model, reffe)

        coeff0 = rand(Float64)
        coeffs = rand(SVector{D,Float64})
        f(x) = coeffs ⋅ SVector(Tuple(x)) + coeff0
        # TODO: use this mechanism instead to project
        # Francesc Verdugo @fverdugo 13:11
        # a(u,v) = ∫( u*v )dΩ
        # l(v) = a(f,v)
        # Solve a fe problem with this weak form
        # See also tutorial 10, "Isotropic damage model", section "L2
        # projection", function "project"
        fh = interpolate_everywhere(f, V)
        fhcache = return_cache(fh, VectorValue(zeros(D)...))

        # Test Random points
        xs = [VectorValue(rand(D)...) for i in 1:10]
        for x in xs
            x = VectorValue(rand(D)...)
            fx = f(x)
            fhx = evaluate!(fhcache, fh, x)
            @test fhx ≈ fx
        end
        fhxs = fh(xs)
        @test fhxs ≈ f.(xs)

        nv = num_vertices(model) # Number of vertices
        nf = num_faces(model,D-1) # Number of faces
        trian = Triangulation(model)
        topo = GridTopology(model)

        pts = get_vertex_coordinates(topo) # Vertex coordinates
        face_nodes = get_face_nodes(model, D-1) # face-to-node numbering
        face_coords = lazy_map(Broadcasting(Reindex(pts)), face_nodes) # Get LazyArray of coordinates of face

        # Test a random vertex from the triangulation
        pt = pts[rand(1:nv)]
        fhpt = evaluate!(fhcache, fh, pt)
        @test fhpt .≈ f(pt)

        # Test a random point lying on the face of the polytope
        face_coord = face_coords[rand(1:nf)]
        λ = rand(length(face_coord));
        λ = (D > 1) ? λ./sum(λ) : λ
        pt = face_coord ⋅ λ # Point on the face
        fhpt = evaluate!(fhcache, fh, pt)
        @test fhpt .≈ f(pt)

        # Test with CellPoint
        # Build cell_to_fxs manually
        cache1,cache2 = fhcache
        ncells = num_cells(model)
        x_to_cell(x) = CellData._point_to_cell!(cache1, x)
        point_to_cell = map(x_to_cell, xs)
        cell_to_points, point_to_lpoint = make_inverse_table(point_to_cell, ncells)
        cell_to_xs = lazy_map(Broadcasting(Reindex(xs)), cell_to_points)
        cell_to_f = get_array(fh)
        cell_to_fxs = lazy_map(evaluate, cell_to_f, cell_to_xs)

        # Now build CellPoint with xs instead of building cell_to_xs
        cell_point_xs = compute_cell_points_from_vector_of_points(xs, trian, PhysicalDomain())
        cell_point_fxs = evaluate(fh, cell_point_xs)
        @test cell_point_fxs ≈ cell_to_fxs

    end
end

p = QUAD
D = num_dims(QUAD)
et = Float64
source_model = CartesianDiscreteModel((0,1,0,1),(2,2))

@testset "Test interpolation Lagrangian" begin
  # Lagrangian space -> Lagrangian space
  f(x) = x[1] + x[2]
  reffe = LagrangianRefFE(et, p, 1)
  V₁ = FESpace(source_model, reffe, conformity=:H1)
  fh = interpolate_everywhere(f, V₁)
  # Target Lagrangian Space
  reffe = LagrangianRefFE(et, p, 2)
  model = CartesianDiscreteModel((0,1,0,1),(4,4))
  V₂ = FESpace(model, reffe, conformity=:H1)

  ifh = Interpolable(fh)
  try
    interpolate_everywhere(fh, V₂)
  catch
    gh = interpolate_everywhere(ifh, V₂)
    pts = [VectorValue(rand(2)) for i=1:10]
    for pt in pts
      @test gh(pt) ≈ fh(pt)
    end
  end

  # VectorValued Lagrangian
  fᵥ(x) = VectorValue([x[1], x[1]+x[2]])
  reffe = ReferenceFE(lagrangian, VectorValue{2,et}, 1)
  V₁ = FESpace(source_model, reffe, conformity=:H1)
  fh = interpolate_everywhere(fᵥ, V₁)
  # Target
  reffe = ReferenceFE(lagrangian, VectorValue{2,et},  2)
  V₂ = FESpace(model, reffe, conformity=:H1)

  ifh = Interpolable(fh);
  gh = interpolate_everywhere(ifh, V₂)
  pts = [VectorValue(rand(2)) for i=1:10]
  for pt in pts
    @test gh(pt) ≈ fh(pt)
  end
end

@testset "Test interpolation RT" begin
  # RT Space -> RT Space
  f(x) = VectorValue([x[1], x[2]])
  reffe = RaviartThomasRefFE(et, p, 0)
  V₁ = FESpace(source_model, reffe, conformity=:HDiv)
  fh = interpolate_everywhere(f, V₁);
  # Target RT Space
  reffe = RaviartThomasRefFE(et, p, 1)
  model = CartesianDiscreteModel((0,1,0,1),(40,40))
  V₂ = FESpace(model, reffe, conformity=:HDiv)

  ifh = Interpolable(fh)
  gh = interpolate_everywhere(ifh, V₂)
  pts = [VectorValue(rand(2)) for i=1:10]
  for pt in pts
    @test gh(pt) ≈ fh(pt)
  end
end

#np = 3
#ndofs = 4
#
#p = Point(1,2)
#x = fill(p,np)
#z = 2.0
#
#v = VectorValue(3.0,1.5)
#w = VectorValue(3.4,3.5)
#a = MockBasis{2}(v,ndofs)
#b = MockBasis{2}(w,ndofs)
#c = fill(1.0,ndofs)
#f = OtherMockBasis{2}(ndofs)
#
#g = MockField{2}(v)
#
#l = 10
#xl = Fill(x,l)
#zl = [ z for  i in 1:l]
#cl = fill(c,l)
#fl = Fill(f,l)
#ϕl = lincomb(fl,cl)
#gl = fill(g,l)
#al = Fill(a,l)
#bl = fill(b,l)
#
#gf = GenericCellField(gl,ϕl,Val(true))
#gf_x = evaluate(gf,xl)
#∇gf_x = evaluate(∇(gf),xl)
#test_cell_field(gf,xl,gf_x,grad=∇gf_x)
#
#af = GenericCellField(al,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
#bf = GenericCellField(bl,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
#zf = convert_to_cell_field(zl,ϕl)
#df = af*zf
#dft = trialize_cell_basis(df)
#
#zfr = reindex(zf,[1,4,3])
#@test length(zfr) == 3
#@test length(get_array(zfr)) == 3
#@test length(get_cell_map(zfr)) == 3
#@test length(get_cell_axes(zfr)) == 3
#
## Check memoization
#df_x1 = evaluate(df,xl)
#df_x2 = evaluate(df,xl)
#@test df_x1 === df_x2
#∇gf1 = ∇(gf)
#∇gf2 = ∇(gf)
#@test ∇gf1 === ∇gf2
#@test evaluate(∇gf1,xl) === evaluate(∇gf2,xl)
#εgf1 = ε(gf)
#εgf2 = ε(gf)
#@test εgf1 === εgf2
#@test ∇×gf === ∇×gf
#@test evaluate(∇×gf,xl) === evaluate(∇×gf,xl)
#
#@test is_test(af)
#@test is_trial(dft)
#@test is_basis(af)
#@test is_basis(dft)
#mf = af⋅dft
#@test get_metasize(mf) == (:,:)
#mf_x = evaluate(mf,xl)
#@test size(mf_x[1]) == (np,ndofs,ndofs)
#
#@test get_cell_axes(mf) == Fill((Base.OneTo(ndofs),Base.OneTo(ndofs)),l)
#
#idsL = [ i*collect(1:ndofs)  for i in 1:l]
#idsR = [ 2*i*collect(1:ndofs)  for i in 1:l]
#axesL = Fill((Base.OneTo(ndofs),),l)
#axesR = Fill((Base.OneTo(ndofs),),l)
#
#idsS = merge_cell_dofs_at_skeleton(idsL,idsR,axesL,axesR)
#@test isa(idsS,VectorOfBlockArrayCoo)
#
#afS = merge_cell_fields_at_skeleton(af,2*af)
#@test isa(afS,SkeletonCellField)
#
#afL_x = evaluate(afS.plus,xl)
#afR_x = evaluate(afS.minus,xl)
#@test isa(afL_x,VectorOfBlockArrayCoo)
#@test isa(afR_x,VectorOfBlockArrayCoo)
#
#@test isa(afS*2,SkeletonCellField)
#@test isa(afS+afS,SkeletonCellField)
#
## Checks associated with trial bases
#df = bf
#dft = trialize_cell_basis(df)
#@test is_trial(dft) == true
#cell_vals = [rand(ndofs) for  i in 1:l]
#cell_field = lincomb(df,cell_vals)
#@test isa(cell_field,CellField)
#@test get_metasize(cell_field) == ()
#cell_field_x = evaluate(cell_field,xl)
#@test isa(cell_field_x[1],AbstractVector)

end # module


text:
module EllipsoidAlgebra

using Plots, LinearAlgebra

"""
E(x, P) = { z : (z-x) P ( z-x) <= 1
"""
struct Ellipsoid{Tx, TP}
    x0::Tx
    P::TP
end

Ellipsoid(x0, P, V) = Ellipsoid(x0, P/V)

function plotEllipse!(e::Ellipsoid; kwargs...)

    @assert length(e.x0) == 2

    Pinv_half = real.(sqrt(inv(e.P)))

    x(θ) = e.x0[1] + first(Pinv_half * [cos(θ), sin(θ)])
    y(θ) = e.x0[2] + last(Pinv_half * [cos(θ), sin(θ)])

    plot!(x, y, 0, 2π; kwargs...)

end

function project(e::Ellipsoid, dims=[1,2])

    N = length(e.x0)
    Pinv_half = sqrt(inv(e.P))

    # construct projection matrix
    T = zeros(2, N)
    T[1,dims[1]] = 1
    T[2, dims[2]] = 1

    # get new shape matrix
    TPinvhalf = T * Pinv_half
    PPinv = TPinvhalf * TPinvhalf'
    Pp = real.(inv(PPinv))

    # get new center
    xp0 = T * e.x0

    return Ellipsoid(xp0, Pp)

end

function sampleRand(e::Ellipsoid, N=1)

    @assert isposdef(e.P)

    Pinv_half = sqrt(inv(e.P))

    # sample random from unit sphere, and map to ellipsoid
     y = 2*rand(length(e.x0), N) .- 1
     for i=1:N
        while norm(y[:,i]) > 1.0
            y[:,i] = 2*rand(length(e.x0), 1) .- 1
        end

    end

    xs = similar(y)
    for i=1:N
        xs[:,i] = e.x0 + Pinv_half * y[:,i]
    end

    return xs
end

function creategrid(d::Integer, n::Integer)

    @assert d >= 1 ("d (number of dimensions) must be a positive integer")
    @assert n >= 2 ("n (number of points) must be a at least 2")

    r = range(-1, 1, length = n)

    iter = Iterators.product((r for _ in 1:d)...)

    return vec([collect(i) for i in iter if dot(i, i) <= 1])
end

const sphere_6_grid = creategrid(6, 5)

function sampleUniformGrid(e::Ellipsoid, n=5)
    N = length(e.x0)
    Pinv_half = sqrt(inv(e.P))
    if N == 6
        pts = [e.x0 + Pinv_half * i for i in sphere_6_grid]
        return pts
    else
        pts = [e.x0 + Pinv_half * i for i in creategrid(N, n)]
        return pts
    end
end

end

ext: jl


text:
#========================================================================================#
#	Laboratory 6
#
# Data visualisation
#
# Author: Niall Palfreyman, 27/01/2022
#========================================================================================#
[
	Activity(
		"""
		In this laboratory we learn how to use Julia's Makie plotting package to create
		publication-quality graphics for data visualisation. Makie is the front-end for
		several different backend graphics packages - for example, Cairo, GL and WebGL. We
		shall focus here on Cairo. First, we add and load the CairoMakie package. You already
		know how to do this using the Julia PackageManager, but like everything in Julia, we
		can also do it ourselves from within a program using method calls:

		using Pkg
		Pkg.add("CairoMakie")
		using CairoMakie

		This will take a while - as will also your first graphics call. Julia uses JiT (Just-
		in-time) compilation, so the first time you call a method, it will always take longer
		than the second time you call it. Enter the following command, then tell me the type of
		your return value fig:

		fig = scatterlines(0:10,(0:10).^2);
		""",
		"Enter the command exactly as I have done here",
		x -> x==Main.Makie.FigureAxisPlot
	),
	Activity(
		"""
		If you enter fig at the Julia prompt, you will see that a Makie figure contains a whole
		bunch of fields necessary for producing a graphic image. We will look at this image by
		saving it to a pdf file, but first we need to study the filesystem underlying the Julia
		REPL. Enter now the following command:

		pwd()

		The returned string tells you the path of your Present Working Directory. It should look
		something like this:

		"C:\\\\Users\\\\hswt136nia\\\\...\\\\Projects\\\\Ingolstadt\\\\src"

		In the chapter "Filesystem" of the Julia user manual you will find many functions for
		exploring the filesystem - we shall just use one or two here. Look at the return vale of
		the following command and tell me its type:

		readdir()
		""",
		"It should be a Vector of filenames in your PWD",
		x -> x <: Vector{String}
	),
	Activity(
		"""
		We want to create a Temp directory to save our files into. First, check whether there
		already exists such a directory in the PWD by entering:

		"Temp" in readdir()

		If the result is false, create the Temp subdirectory by entering:

		mkdir("Temp")

		What happens if you use this command when Temp already exists?
		""",
		"",
		x -> occursin( "error", lowercase(x)) || occursin( "exist", lowercase(x))
	),
	Activity(
		"""
		Now change directory into your new Temp directory, then check that you've changed PWD:

		cd("Temp")
		pwd()

		Now we're ready to save our Makie figure. Enter the following command, go and look at
		the resulting file using Adobe Acrobat, then come back here and tell me the highest
		number on the vertical axis:

		save("myfig.png", fig)
		""",
		"",
		x -> x==100
	),
	Activity(
		"""
		Didn't it look pretty? :) Wouldn't it be nice to view the graphic from within the Julia
		REPL? The REPL is a very basic environment that cannot display graphics, but we can easily
		include this functionality by adding and loading the package ElectronDisplay. When you have
		done this, just enter fig at the Julia prompt and ElectronDisplay will automatically
		display() your image.

		What is the type of the return value from this display?
		""",
		"",
		x -> x <: Main.Makie.FigureAxisPlot
	),
	Activity(
		"""
		Our figure is still a little primitive - let's customise its attributes. First notice
		that fig has a structure that includes a figure, an axis and a plot object:

		fg,ax,plt = fig;

		Check out the attributes field of plt to find out how many attributes it has:
		""",
		"",
		x -> x >= 15
	),
	Activity(
		"""
		We can ask for help on any function in the REPL by using the help() function. Try:

		help(lines)

		and then tell me how many different ways there are of calling the lines() function:
		""",
		"The signature of a function is the pattern of arguments that defines which method to call",
		x -> x==3
	),
	Activity(
		"""
		OK, now let's get fancy. We don't always want the scatter points on our plot, so we'll try
		out the lines() function. Also, it would be nice to plot something a little more exciting,
		like maybe the Hill function. The Hill function describes the activation and inhibition of
		DNA expression in biological cells by a transcription factor (TF):

		function hill( tf, K::Real=2, n::Real=1)
			abs_n = abs(n)
			if abs_n!=1 K = K^abs_n; tf = tf.^abs_n end
			n>=0 ? tf./(K.+tf) : K./(K.+tf)
		end
		trfactor = 0:20
		lines(trfactor,hill(trfactor))

		We can brighten up this figure by adding some interesting line attributes:

		lines( trfactor, hill(trfactor); color=:blue, linewidth=3, linestyle=:dash)

		The ; symbol marks the beginning of keyword arguments of a function. Experiment with this
		plot by changing the keyword arguments of lines(). What argument value displays a line
		consisting of a sequence of two dots and a dash (-..-..-..-)?
		""",
		"Search for keyword arguments on the site: https://makie.juliaplots.org/",
		x -> x==:dashdotdot
	),
	Activity(
		"""
		The lines() function returns a Plot object: the line that is to be plotted. This line is
		positioned within a set of axes: an Axis object. We can design the structure of this Axis
		object by adding an Axis specification to lines():

		lines( trfactor, hill(trfactor,2,1); color=:blue, label="activation",
			linewidth=3, linestyle=:dash,
			axis=(;
				xlabel="TF concentration", ylabel="Expression rate",
				title="Transcription regulation",
				xgridstyle=:dash, ygridstyle=:dash
			)
		)

		Notice carefully how I have used indentation to make it clear where each structure or
		function begins and ends. It's not easy to do this in the REPL, but soon we will be writing
		our own code, so it's good to practise making our code easily readable by other programmers.

		Experiment now with the axis parameters ...
		""",
		"",
		x -> true
	),
	Activity(
		"""
		This Axis object is then displayed within a Figure object. We can design the appearance of
		this figure by adding a Figure specification to lines():

		lines( trfactor, hill(trfactor); color=:blue, label="activation",
			linewidth=3, linestyle=:dash,
			axis=(;
				xlabel="TF concentration", ylabel="Expression rate",
				title="Transcription regulation",
				xgridstyle=:dash, ygridstyle=:dash
			),
			figure=(;
				figure_padding=5, resolution=(600,400), font="sans",
				backgroundcolor=:green, fontsize=20
			)
		)

		Experiment now with the figure parameters ...
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Once we have set up a Plot object inside an Axis, inside a Figure, we can add extra Plot
		objects using bang! methods that change the current state of an object, for example:

		scatterlines!( trfactor, hill(trfactor,2,-1); color=:red, label="repression", linewidth=3)
		axislegend( "Legend"; position=:rc)

		Once we have done this, we can redisplay the current figure like this:

		current_figure()

		Take a look at the point where the curves of the two Hill-functions cross. You can see that
		this TF concentration correspondes to the value of K that we pass to the Hill function.
		Check that this is generally true, and tell me what value of Expression rate corresponds
		in general to the TF concentration value K:
		""",
		"Create a new graphic, and set the new value of K in BOTH calls to hill()",
		x -> x==0.5
	),
	Activity(
		"""
		Let's use our new-found knowledge of plotting to display a bubble-plot:

		data = randn(50,3)
		figr, ax, pltobj = scatter( data[:,1], data[:,2];
			color=data[:,3], label="Bubbles", colormap=:plasma, markersize=15*abs.(data[:,3]),
			axis=(; aspect=DataAspect()),
			figure=(; resolution=(600,400))
		)
		limits!(-3,3,-3,3)
		Legend( fig[1,2], ax, valign=:top)
		Colorbar( fig[1,2], pltobj, height=Relative(3/4))
		figr

		Do you like my pretty bubble-plot?
		""",
		"I'm very insecure: I won't let you go on until you say you like my art!",
		x -> occursin("yes",lowercase(x))
	),
	Activity(
		"""
		As you see, to create a nice graphic, we often need to build it up in several steps, so it
		is a idea to collect these steps together inside a single function. For example, suppose
		we want to display several different 2-d functions as heatmaps - we might do like this:

		function prettyheatmap( f::Function)
			figure = (; resolution=(600,400), font="CMU Serif")
			axis = (; xlabel="x", ylabel="y", aspect=DataAspect())
			xs = range(-2, 2, length = 25)
			ys = range(-2, 2, length = 25)
			zs = [f(x,y) for x in xs, y in ys]

			global fig, ax, pltobj = heatmap( xs, ys, zs, axis=axis, figure=figure)
			Colorbar( fig[1,2], pltobj, label="A colour bar")
			fig
		end

		Now call prettyheatmap() as follows, then tell me how many yellow rings you see:

		prettyheatmap( (x,y) -> cos(5*hypot(x,y)) )
		""",
		"Just count the complete rings - not the circle",
		x -> x==2
	),
	Activity(
		"""
		Now that we have wrapped our plotting code inside the function prettyheatmap(), we can
		easily reuse this encapsulated code to display a completely different function:

		mountains(x,y) =
			3 * (1-x)^2 *				exp(-(x^2+(y+1)^2))
			-(1/3) *					exp(-((x+1)^2 + y^2))
			-10 * (x/5 - x^3 - y^5) *	exp(-(x^2+y^2))
		prettyheatmap( mountains)

		How many mountains can you see?
		""",
		"Use the colour bar to count the number of patches where z > 0",
		x -> x==3
	),
	Activity(
		"""
		You see how useful it is to encapsulate (i.e.: hide) code? As soon as we have hidden the
		complicated graphics code inside the function myheatmap(), we can easily reuse this
		encapsulated code, by simply passing to it the name of the function we wish to plot:

		Encapsulation is THE CENTRAL tool for abstraction in modern computer languages!

		However, encapsulation can get broken! Before we move on to the next laboratory, I want
		you to notice something important. Do you remember our graph of the quadratic function
		y=x^2 that we saved in a file at the beginning of this laboratory? If you remember, we
		stored this graphic in a global variable fig, so we should be able to look at it again
		now, yes?

		OK, so please now display the quadratic graph contained in the global variable fig, and
		then, without doing anything else, move on to the next activity in this laboratory:

		fig
		""",
		"Your result may puzzle you, but please just move on to the next activity",
		x -> true
	),
	Activity(
		"""
		You should have found that the variable fig no longer contains the quadratic graph, but
		instead contains our most recent graph of the mountains() function.

		I want you to think now about this question: Why has the value of fig changed, and what
		single word can you delete in this laboratory to fix this problem? Don't worry if you
		can't immediately work it out: we will return to this important issue of encapsulation
		in laboratory 7.
		""",
		"",
		x -> true
	),
]


text:
export timeslotmat
export adjustLossWeights

"""
    timeslotmat(matrix::AbstractMatrix, timestamp::AbstractVector; dim=2, slotvalue=1.0)

mark `matrix` with `timestamp`, all elements of `timestamp` is ∈ (0,1), standing for time ratio

# Example
    julia> timeslotmat(reshape(1:24,2,12), [0.2, 0.8], dim=2, slotvalue=0)
    2×12 Array{Int64,2}:
    1  3  0  7   9  11  13  15  17  0  21  23
    2  4  0  8  10  12  14  16  18  0  22  24
          ↑                         ↑
          20% position              80% position
"""
function timeslotmat(matrix::AbstractMatrix{S1}, timestamp::AbstractVector{S2}; dim=2, slotvalue=1.0) where {S1,S2}
    x = copy(matrix)
    T = size(x, dim)               # time dimention
    I = ceil.(Int, T .* timestamp) # map time in (0,1) to integer [1,T]
    if dim==2
        for t in I
            x[:,t] .= slotvalue
        end
    elseif dim==1
        for t in I
            x[t,:] .= slotvalue
        end
    else
        @error "dim is 1 or 2, but got $dim"
    end
    return x
end

"""
    adjustLossWeights(x...) -> w

`x` is the loss values of multiple losses at training step i, so the weights for
each loss function at training step (i+1) is `w` , `w` is from:

    yᵢ = 1 / xᵢ
    wᵢ = yᵢ / ∑ yᵢ
"""
function adjustLossWeights(x...)
    n = length(x)
    w = zeros(n)
    for i = 1:n
        w[i] = 1 / x[i]
    end
    return w ./ sum(w)
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
@testset "Testing Inner Product" begin
    # Note: use `abs` to serve as an example of an activation function that cannot be fused
    # with the `innerproduct` kernel.
    activations = [identity, Flux.relu, Flux.sigmoid, abs]
    for σ in activations
        @show σ
        m = Flux.Dense(128, 256, σ)
        M = OneDNN.Dense(m)

        x = randn(Float32, 128, 128)
        y = m(x)

        X = OneDNN.Memory(x)
        @test M.optimized_weights == false
        Y = M(X)
        # Did the automatic weight transformation take place?
        @test M.optimized_weights
        @test isapprox(y, OneDNN.materialize(Y))

        # Ensure inference works correctly
        @inferred M(X)

        # Now - try backprop
        y, back_jl = Zygote._pullback(m, x)
        Y, back_dnnl = Zygote._pullback(M, X)
        @test isa(Y, OneDNN.Memory)
        @test isapprox(y, OneDNN.materialize(Y))
        @inferred Zygote._pullback(M, X)
        out, back = Zygote._pullback(M, X)
        @inferred back(out)

        dy = Float32(0.125) * randn(Float32, size(y))
        DY = OneDNN.Memory(dy)

        grads_jl = back_jl(dy)
        grads_dnnl = back_dnnl(DY)
        @test isa(grads_jl, Tuple)
        @test isa(grads_dnnl, Tuple)

        # N.B.: Notice the difference in the plurality of the `weights` field ...
        # In the OneDNN code, we try to stay with the plural version since that's what the
        # oneDNN source code generally does.
        @test isapprox(
            grads_jl[1].weight, transpose(OneDNN.materialize(grads_dnnl[1].weights))
        )
        @test isapprox(grads_jl[1].bias, OneDNN.materialize(grads_dnnl[1].bias))
        @test isapprox(grads_jl[2], OneDNN.materialize(grads_dnnl[2]))
        @inferred back_dnnl(DY)
    end
end


text:
module CBindingGen
	using Clang
	using Clang: LibClang
	# using Intervals
	using Todo

	export ConverterContext, generate

	struct CodeLocation
		file::String
		line::Int
		col::Int
	end
	CodeLocation(decl::CLCursor) = CodeLocation(decl.cursor)
	CodeLocation(decl::Clang.LibClang.CXCursor) = CodeLocation(location(decl))
	function CodeLocation(loc::Clang.LibClang.CXSourceLocation)
		cxfile = Ref{Clang.LibClang.CXFile}()
		line = Ref{Cuint}()
		col = Ref{Cuint}()
		offset = Ref{Cuint}()
		GC.@preserve loc begin
			Clang.LibClang.clang_getSpellingLocation(loc, cxfile, line, col, offset)
			cxstr = Clang.LibClang.clang_getFileName(cxfile[])
			file = try
				ptr = Clang.LibClang.clang_getCString(cxstr)
				ptr == C_NULL ? "<unknown>" : unsafe_string(ptr)
			finally
				Clang.LibClang.clang_disposeString(cxstr)
			end
		end
		return CodeLocation(file, line[], col[])
	end

	Base.isless(a::CodeLocation, b::CodeLocation) = a.file == b.file && (a.line < b.line || (a.line == b.line && a.col < b.col))
	Base.show(io::IO, cl::CodeLocation) = print(io, "$(cl.file):$(cl.line):$(cl.col)")

	# Intervals.Interval(decl::CLCursor) = Interval(decl.cursor)
	# Intervals.Interval(decl::Clang.LibClang.CXCursor) = Interval(extent(decl))
	# Intervals.Interval(ext::Clang.LibClang.CXSourceRange) = Interval(CodeLocation(Clang.LibClang.clang_getRangeStart(ext)), CodeLocation(Clang.LibClang.clang_getRangeEnd(ext)))

	struct JuliaizedC
		decl::Union{CLCursor, Nothing}
		expr::String
		kind::Symbol  # :atdevelop, :atcompile, :atcompile_typedefs, :atcompile_bindings, :atload

		JuliaizedC(decl::CLCursor, expr::String, kind::Symbol) = new(decl, expr, kind)
	end

	# struct MacroUse
	# 	c::String
	# 	loc::Interval{CodeLocation}
	# end

	_gensym(ctx, str) = join(("", string(str), string(hash(str), base = 16, pad = 16)), '_')
	_export(ctx, str) = !isempty(str) && !startswith(str, '_') && push!(ctx.exports, str)

	struct ConverterContext <: AbstractContext
		ctx::DefaultContext
		quiet::Bool
		filt::Function
		libs::Union{Vector{String}, Nothing}   # the library name (or library path) to dlopen, an array of library paths or names, or `nothing` for dlopening the running Julia process itself
		includes::Vector{String}
		args::Vector{String}
		exports::Vector{String}
		oneofs::Set{String}
		typedefs::Dict{String, String}
		converted::Vector{JuliaizedC}
		# macroUses::Vector{MacroUse}

		function ConverterContext(filter::Function, libs::Union{Vector{String}, Nothing} = String[], includes::Vector{String} = String[], args::Vector{String} = String[]; quiet::Bool = false)
			return new(
				DefaultContext(),

				quiet,

				filter,
				libs,
				includes,
				args,

				String[],
				Set{String}(),
				Dict{String, String}(),
				JuliaizedC[],
				# MacroUse[],
			)
		end
	end

	Clang.parse_header!(ctx::ConverterContext, header::String; kwargs...) = Clang.parse_headers!(ctx, [header]; kwargs...)
	function Clang.parse_headers!(ctx::ConverterContext, headers::Vector{String}; args::Vector{String} = String[], includes::Vector{String} = String[], builtin::Bool = false)
		(o, c)  = builtin ? ('<', '>') : ('"', '"')

		return mktempdir() do dir
			hdr = joinpath(dir, "parse_headers.h")
			open(hdr, "w+") do file
				for header in headers
					print(file,
						"""
						#include $(o)$(header)$(c)
						"""
					)
				end
			end

			tu = parse_header(hdr; args = vcat(ctx.args, args), includes = vcat(Clang.CLANG_INCLUDE, ctx.includes, includes), flags = Clang.LibClang.CXTranslationUnit_DetailedPreprocessingRecord)
			num = LibClang.clang_getNumDiagnostics(tu.ptr)
			for i in 1:num
				diag = LibClang.clang_getDiagnostic(tu.ptr, i-1)
				sev = LibClang.clang_getDiagnosticSeverity(diag)
				sev in (LibClang.CXDiagnostic_Error, LibClang.CXDiagnostic_Fatal) && error("Errors encountered parsing headers, unable to proceed with conversion")
			end

			_convert(ctx, tu)
		end

		return
	end

	function generate(ctx::ConverterContext, where::Union{AbstractString, Nothing} = nothing, prefix::Union{AbstractString, Nothing} = nothing)
		if isnothing(where)
			result = quote end
		else
			prefix = prefix === nothing ? "" : prefix*"-"

			try
				mkpath(where)
			catch
			end
		end

		for kind in (:atdevelop, :atcompile, :atcompile_typedefs, :atcompile_bindings, :atload)
			f = function (io)
				if kind === :atcompile
					for i in 1:10:length(ctx.exports)
						expr = "export $(join(ctx.exports[i:min(length(ctx.exports), i+9)], ", "))"
						if isnothing(where)
							push!(result.args, Meta.parse(expr))
						else
							println(io, expr)
						end
					end
				elseif kind === :atload
					libs = isnothing(ctx.libs) ? ("CBinding.Clibrary()",) : map(lib -> "CBinding.Clibrary($(repr(lib)))", ctx.libs)
					expr = "$(_gensym(ctx, "libraries")) = ($(isempty(libs) ? "" : join(libs, ", ")*","))"
					if isnothing(where)
						push!(result.args, Meta.parse(expr))
					else
						println(io, expr)
					end
				end

				for c in ctx.converted
					if c.kind === kind
						if isnothing(where)
							push!(result.args, Meta.parse(c.expr))
						else
							if !isnothing(c.decl)
								println(io, "# ", CodeLocation(c.decl))
							end
							println(io, c.expr)
						end
					end
				end
			end

			if isnothing(where)
				f(nothing)
			else
				open(f, joinpath(where, "$(prefix)$(kind).jl"), "w+")
			end
		end

		return isnothing(where) ? result : nothing
	end

	function _convert(ctx::ConverterContext, tu::TranslationUnit)
		root = getcursor(tu)

		for decl in children(root)
			if decl isa CLMacroInstantiation
				# macroUse = MacroUse(spelling(decl), Interval(decl))
				# push!(ctx.macroUses, macroUse)
			elseif ctx.filt(decl)
				_convert(ctx, decl)
			end
		end
	end

	_convert(ctx::ConverterContext, decl) = @warn "Not wrapping $(decl)"
	_convert(ctx::ConverterContext, decl::CLLastPreprocessing) = nothing

	include("types.jl")
	include("functions.jl")
	include("variables.jl")
	include("macros.jl")
	include("code.jl")
end


text:
module Spring

end # module


text:
include("../examples/rijke.jl")
include("../src/clvs.jl")
using OrdinaryDiffEq
using PyPlot
	s = [7.0, 0.2]
	#n = 100000
	n = 50000
	d = N
	nRunup = 400000
	u = Rijke(rand(d),s,nRunup)
	du_trj = zeros(d,d,n)
	duT_trj = zeros(d,d,n)
	f_trj = zeros(d,n)
	for i = 1:n
		f!(view(f_trj,:,i), u, s, 1.)
		du_trj[:,:,i] = dRijke(u, s, 1.e-8)
		duT_trj[:,:,n+1-i] = du_trj[:,:,i]'
		u .= Rijke(u, s, 1)
	end
	d_u = 20
	les, tan_clvs = clvs(du_trj, d_u)
	les1, adj_clvs = clvs(duT_trj, d_u)

	tan_angles = zeros(d_u, d_u, n)
	adj_angles = zeros(d_u, d_u, n)

	for i = 1:n
		tan_angles[:,:,i] = tan_clvs[:,:,i]'*tan_clvs[:,:,i]
		adj_angles[:,:,n+1-i] = adj_clvs[:,:,i]'*
								adj_clvs[:,:,i]
	end
	n_spinup = 5000
	tan_adj_angles = zeros(d_u, d_u, n-1)
	for i = 2:n
		tan_adj_angles[:,:,i-1] = tan_clvs[:,:,i]'*
		adj_clvs[:,:,n+2-i]
	end
	tan_adj_angles = tan_adj_angles[:,:,n_spinup:n-1-n_spinup]
	mean_tan_adj_angles = 180*acos.(sum(tan_adj_angles,
										dims=3)/n)/pi
	mean_tan_adj_angles = mean_tan_adj_angles[:,:,1]

	mean_tan_angles = sum(tan_angles,dims=3)/n
	mean_tan_angles = mean_tan_angles[1:6,1:6,1]
	mean_tan_angles = 180*acos.(mean_tan_angles)/pi
	mean_adj_angles = sum(adj_angles,dims=3)/n
	mean_adj_angles = mean_adj_angles[1:6,1:6,1]
	mean_adj_angles = 180*acos.(mean_adj_angles)/pi

	fig, ax = subplots(1,1)
	ax.xaxis.set_tick_params(labelsize=25)
	ax.yaxis.set_tick_params(labelsize=25)
	angle_plot = ax.matshow(mean_tan_angles)
	cbar = fig.colorbar(angle_plot)
	cbar.ax.tick_params(labelsize=25)
	fig2, ax2 = subplots(1,1)
	ax2.xaxis.set_tick_params(labelsize=25)
	ax2.yaxis.set_tick_params(labelsize=25)
	angle_plot2 = ax2.matshow(mean_adj_angles)
	cbar = fig2.colorbar(angle_plot2)
	cbar.ax.tick_params(labelsize=25)

	fig1, ax1 = subplots(1,1)
	ax1.xaxis.set_tick_params(labelsize=25)
	ax1.yaxis.set_tick_params(labelsize=25)
	ax1.set_title(L"Probability distribution of
				  angles between $V^1$ and $V^2$",
				  fontsize=25)
	ax1.hist(acos.(tan_angles[1,2,:])*180/pi,
			 density=true,bins=100)
	# minimum angle between 1 and 2 is 15.6 degrees
	ax1.set_xlabel("Angle in degrees",fontsize=25)
	ax1.grid(true)
	fig3, ax3 = subplots(1,1)
	ax3.xaxis.set_tick_params(labelsize=25)
	ax3.yaxis.set_tick_params(labelsize=25)
	ax3.set_title(L"Probability distribution of
				  angles between $W^1$ and $W^2$",
				  fontsize=25)
	ax3.hist(acos.(adj_angles[1,2,:])*180/pi,
			 density=true,bins=100)
	# minimum angle between 1 and 2 is 4 degrees
	ax3.set_xlabel("Angle in degrees",fontsize=25)
	ax3.grid(true)

	fig4, ax4 = subplots(1,1)
	ax4.xaxis.set_tick_params(labelsize=25)
	ax4.yaxis.set_tick_params(labelsize=25)
	angle_plot4 = ax4.matshow(mean_tan_adj_angles)
	cbar = fig4.colorbar(angle_plot4)
	cbar.ax.tick_params(labelsize=25)

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
##################################################
# Compute the NR orbit-averaged diffusion coefficients in (Jr,L)
##################################################

include("../sources/julia/Main.jl")

Jr = 0.2 # Radial action
L = 0.3 # Angular momentum
m = 0.1 # Mass of field star

println("--------------------")

println("Plummer potential")

println("q  = ", qCalc)
println("m  = ", m)
println("Jr = ", Jr)
println("L  = ", L)

println("--------------------")

println("Coulomb log = ", logCoulomb)

println("--------------------")

@time dJr, dL, dJrJr, dJrL, dLL = avr_action_coefficients(Jr,L,m)

println("D_Jr  =  ", dJr)
println("D_JrJr =  ", dJrJr)
println("D_L  =  ", dL)
println("D_LL =  ", dLL)
println("D_JrL =  ", dJrL)

println("--------------------")

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
#-----------------------------------------------------------------------# Stats
n = 1000
x,  y,  z  = rand(Bool, n), randn(n), rand(1:10, n)
x2, y2, z2 = rand(Bool, n), randn(n), rand(1:10, n)
xs, ys, zs = vcat(x, x2), vcat(y, y2), vcat(z, z2)

p = 5
xmat,  ymat,  zmat  = rand(Bool, n, p), randn(n, p), rand(1:10, n, p)
xmat2, ymat2, zmat2 = rand(Bool, n, p), randn(n, p), rand(1:10, n, p)

function mergestats(a::OnlineStat, y1, y2; nobs_equals_length=true)
    b = copy(a)
    fit!(a, y1)             # fit a on y1
    fit!(b, y2)             # fit b on y2
    merge!(a, b)            # merge b into a
    fit!(b, y1)             # fit b on y1
    if nobs_equals_length
        @test nobs(a) == nobs(b) == length(y1) + length(y2)
    else
        @test nobs(a) == nobs(b)
    end
    a, b
end
mergevals(o1::OnlineStat, y1, y2; kw...) = map(value, mergestats(o1, y1, y2; kw...))

@testset "Testing Stats" begin
#-----------------------------------------------------------------------------# CircBuff
println("  > CircBuff")
@testset "CircBuff" begin
    o = CircBuff(Int, 5)
    fit!(o, 1:2)
    @test value(o) == [1,2]
    fit!(o, 3:11)
    @test o[1] == 7
    @test o[end] == 11
end

#-----------------------------------------------------------------------# Counter
println("  > Counter")
@testset "Counter" begin
    o = fit!(Counter(Int), 1:10)
    @test value(o) == 10
    o2 = fit!(Counter(Int), 1)
    @test value(merge!(o, o2)) == 11
    ==(mergevals(Counter(), y, y2)...)
end

#-----------------------------------------------------------------------# CountMap
println("  > CountMap")
@testset "CountMap" begin
    a = fit!(CountMap(Bool), x)
    @test sort(value(a)) == sort(countmap(x))
    @test O.pdf(a, true) == mean(x)
    @test O.pdf(a, false) == mean(!, x)
    @test O.pdf(a, 2) == 0.0
    @test keys(a) == keys(a.value)
    @test values(a) == values(a.value)
    @test a[true] == a.value[true]
    @test O.nkeys(a) == 2

    b = fit!(CountMap(Int), z)
    @test sort(value(b)) == sort(countmap(z))
    for i in 1:11
        @test O.pdf(b, i) ≈ sum(==(i), z) / n
    end
    @test all(x -> 0 < x < 1, O.probs(b))
    @test O.probs(b, 11:13) == zeros(3)

    @test ==(mergevals(CountMap(Bool), x, x2)...)
    @test ==(mergevals(CountMap(Int), z, z2)...)
    @test ==(mergevals(CountMap{Bool}(), x, x2)...)
    @test ==(mergevals(CountMap{Int}(), z, z2)...)
    @test ==(mergevals(CountMap(Dict{Bool,Int}()), x, x2)...)
    @test ==(mergevals(CountMap(Dict{Int,Int}()), z, z2)...)
    # Pair method
    @test ==(mergevals(CountMap(Bool), Pair.(x,z), Pair.(x2,z2); nobs_equals_length=false)...)
    @test ==(mergevals(CountMap(Int), Pair.(z,z), Pair.(z2,z2); nobs_equals_length=false)...)
end
#-----------------------------------------------------------------------# CountMissing
println("  > CountMissing")
@testset "CountMissing" begin
    data = Vector{Union{Missing,Float64}}(copy(y))
    data2 = Vector{Union{Missing,Float64}}(copy(y2))
    data[x] .= missing
    data2[x2] .= missing
    a, b = mergevals(CountMissing(Mean()), data, data2)
    @test value(a.stat) ≈ value(b.stat)
    @test a.nmissing == b.nmissing
end
#-----------------------------------------------------------------------# CovMatrix
println("  > CovMatrix")
@testset "CovMatrix" begin
    o = fit!(CovMatrix(), eachrow(ymat))
    @test O.nvars(o) == size(ymat, 2)
    @test value(o) ≈ cov(ymat)
    @test cov(o) ≈ cov(ymat)
    @test cor(o) ≈ cor(ymat)
    @test all(x -> ≈(x...), zip(var(o), var(ymat; dims=1)))
    @test all(x -> ≈(x...), zip(std(o), std(ymat; dims=1)))
    @test all(x -> ≈(x...), zip(mean(o), mean(ymat; dims=1)))

    @test ≈(mergevals(CovMatrix(), eachrow(ymat), eachrow(ymat2))...)
    @test ≈(mergevals(CovMatrix(), O.eachcol(ymat'), O.eachcol(ymat2'))...)
    @test ≈(mergevals(CovMatrix(Complex{Float64}), eachrow(ymat * im), eachrow(ymat2))...)
    @test ≈(mergevals(CovMatrix(Complex{Float64}), eachrow(ymat * im), eachrow(ymat2 * im))...)
end
#-----------------------------------------------------------------------# Extrema
println("  > Extrema")
@testset "Extrema" begin
    o = fit!(Extrema(), y)
    @test extrema(o) == extrema(y)
    @test minimum(o) == minimum(y)
    @test maximum(o) == maximum(y)

    f(x) = (
        min = minimum(x),
        max = maximum(x),
        nmin = sum(collect(x) .== minimum(x)),
        nmax = sum(collect(x) .== maximum(x))
    )

    @test value(fit!(Extrema(Bool), x)) == f(x)
    @test value(fit!(Extrema(Int), z)) == f(z)

    @test ==(mergevals(Extrema(), y, y2)...)
    @test ==(mergevals(Extrema(), x, x2)...)

    o = fit!(Extrema(Date), Date(2010):Day(1):Date(2011))
    @test minimum(o) == Date(2010)
    @test maximum(o) == Date(2011)

    @test value(fit!(Extrema(Char), 'a':'z')) == f('a':'z')
    @test value(fit!(Extrema(Char), "abc")) == f("abc")
    @test value(fit!(Extrema(String), ["a", "b"])) == f(["a","b"])

    o = fit!(Extrema(), x)
    @test o.nmin == length(x) - sum(x)
    @test o.nmax == sum(x)
end
#-----------------------------------------------------------------------# Group
println("  > Group")
@testset "Group" begin
    o = fit!(5Mean(), eachrow(ymat))
    @test o[1] == first(o)
    @test o[end] == last(o)
    @test 5Mean() == 5Mean()
    @test collect(map(value, value(o))) ≈ vec(mean(ymat, dims=1))

    o2 = Group(m1=Mean(), m2=Mean(), m3=Mean(), m4=Mean(), m5=Mean())
    fit!(o2, eachrow(ymat))
    @test collect(map(value, value(o2))) ≈ vec(mean(ymat, dims=1))
    @test length(o2) == 5

    a, b = mergevals(
        Group(Mean(), Variance(), Sum(), Moments(), Mean()), eachrow(ymat), eachrow(ymat2)
    )
    for (ai, bi) in zip(a, b)
        @test value(ai) ≈ value(bi)
    end

    @test length(values(a)) == 5

    c = fit!(Group([Mean(), Mean()]), zip(1:10, 1:10))
    for m in c
        @test value(m) ≈ mean(1:10)
    end
end
#-----------------------------------------------------------------------# GroupBy
println("  > GroupBy")
@testset "GroupBy" begin
    @test GroupBy(Bool, Mean()) == GroupBy(Bool, Mean())
    d = value(fit!(GroupBy(Bool, Mean()), zip(x,y)))
    @test value(d[true]) ≈ mean(y[x])
    @test value(d[false]) ≈ mean(y[map(!,x)])
    string(GroupBy(Bool, Mean()))

    a, b = mergevals(GroupBy(Int, Mean()), zip(z,y), zip(z2, y2))
    for (ai,bi) in zip(values(sort(a)), values(sort(b)))
        @test value(ai) ≈ value(bi)
    end
    @test value(a[1]) ≈ value(b[1])
    string(GroupBy(Int, Mean()))
end
#-----------------------------------------------------------------------# Mean
println("  > Mean")
@testset "Mean" begin
    o = fit!(Mean(), y)
    @test value(o) ≈ mean(y)
    @test mean(o) ≈ mean(y)
    @test ≈(mergevals(Mean(), y, y2)...)
end
#-----------------------------------------------------------------------# Moments
println("  > Moments")
@testset "Moments" begin
    o = fit!(Moments(), y)
    @test value(o) ≈ [mean(y), mean(y .^ 2), mean(y .^ 3), mean(y .^ 4)]
    @test mean(o) ≈ mean(y)
    @test var(o) ≈ var(y)
    @test std(o) ≈ std(y)
    @test skewness(o) ≈ skewness(y)
    @test kurtosis(o) ≈ kurtosis(y)
    for (v1,v2) in zip(mergevals(Moments(), y, y2)...)
        @test v1 ≈ v2
    end
end

#-----------------------------------------------------------------------# Series/FTSeries
println("  > Series/FTSeries")
@testset "Series/FTSeries" begin
    @testset "Series" begin
        a, b = mergevals(Series(Mean(), Variance()), y, y2)
        @test a[1] ≈ b[1]
        @test a[2] ≈ b[2]

        a, b = mergevals(Series(m=Mean(), v=Variance()), y, y2)
        @test a.m ≈ b.m
        @test a.v ≈ b.v
    end

    @testset "FTSeries" begin
        o = fit!(FTSeries(Mean(); transform=abs), y)
        @test value(o)[1] ≈ mean(abs, y)

        data = vcat(y, fill(missing, 20))
        o = fit!(FTSeries(Union{Number,Missing}, Mean(); transform=abs, filter=!ismissing), data)
        @test value(o)[1] ≈ mean(abs, y)
        @test o.nfiltered == 20

        o2 = fit!(FTSeries(Union{Missing,Float64}, Mean(); transform=abs, filter=!ismissing), data)
        @test value(o2)[1] ≈ mean(abs, y)
        @test o2.nfiltered == 20

        a, b = mergestats(FTSeries(Union{Missing,Number}, Mean(), Variance(); transform=abs, filter=!ismissing), y, y2)
        @test a.nfiltered == b.nfiltered
    end
end

#-----------------------------------------------------------------------------# SkipMissing
println("  > SkipMissing")
@testset "SkipMissing" begin
    data = [rand() > .5 ? missing : rand() for i in 1:1000]
    o = fit!(skipmissing(Mean()), data)
    @test nobs(o) == sum(!ismissing, data)
    @test value(o) ≈ mean(skipmissing(data))
end

#-----------------------------------------------------------------------# Sum
println("  > Sum")
@testset "Sum" begin
    @test sum(fit!(Sum(Int), x)) == sum(x)
    @test value(fit!(Sum(), y)) ≈ sum(y)
    @test value(fit!(Sum(Int), z)) == sum(z)

    @test ==(mergevals(Sum(Int), x, x2)...)
    @test ≈(mergevals(Sum(), y, y2)...)
    @test ==(mergevals(Sum(Int), z, z2)...)
end
#-----------------------------------------------------------------------# Variance
println("  > Variance")
@testset "Variance" begin
    o = fit!(Variance(), y)
    @test mean(o) ≈ mean(y)
    @test var(o) ≈ var(y)
    @test std(o) ≈ std(y)

    @test ≈(mergevals(Variance(), x, x2)...)
    @test ≈(mergevals(Variance(), y, y2)...)
    @test ≈(mergevals(Variance(Float32), Float32.(y), Float32.(y2))...; atol=.001)
    @test ≈(mergevals(Variance(), z, z2)...)
    # Issue 116
    @test std(Variance()) == 1
    @test std(fit!(Variance(), 1)) == 1
    @test std(fit!(Variance(), [1, 2])) == sqrt(.5)
    # https://github.com/joshday/OnlineStats.jl/issues/217
    @test value(fit!(Variance(Float32), randn(Float32, 10))) isa Float32
end

end # end "Test Stats"


text:
let
  i = Index(2)
  A = randomITensor(i, i', i'')
  for T in (Float64, ComplexF64)
    A = randomITensor(T, i, i')
    B = randomITensor(T, i', i'')
    C = A * B
    U, S, V = svd(A, i)
    A, B = factorize(A, i)
  end
end


text:
using Documenter, Mongoc

makedocs(
    format = :html,
    sitename = "Mongoc.jl",
    modules = [ Mongoc ],
    pages = [ "Home" => "index.md",
              "Tutorial" => "tutorial.md",
              "CRUD Operations" => "crud.md",
              "Authentication" => "authentication.md",
              "Transactions" => "transaction.md",
              "API Reference" => "api.md" ]
)

deploydocs(
    repo = "github.com/felipenoris/Mongoc.jl.git",
    target = "build",
    julia  = "1.0",
    deps   = nothing,
    make   = nothing
)


text:
macro mustimplement(sig)
    fname = sig.args[1]
    arg1 = sig.args[2]
    if isa(arg1,Expr)
        arg1 = arg1.args[1]
    end
    :($(esc(sig)) = error(typeof($(esc(arg1))),
                          " must implement ", $(Expr(:quote,fname))))
end
export @mustimplement

macro returntype(expr)
  esc(:(eltype([$expr for i in []])))
end
export @returntype

eltypes{T1}(::Type{Tuple{T1}}) = (T1,)
eltypes{T1,T2}(::Type{Tuple{T1,T2}}) = (T1, T2)
eltypes{T1,T2,T3}(::Type{Tuple{T1,T2,T3}}) = (T1, T2, T3)
eltypes{T1,T2,T3,T4}(::Type{Tuple{T1,T2,T3,T4}}) = (T1, T2, T3, T4)

function flip_tuple_array{T<:Tuple}(a::Array{T})
  types = eltypes(T)
  N = length(a)
  b = [Array{t}(N) for t in types]
  for i=1:length(types)
    for j=1:N
      b[i][j] = a[j][i]
    end
  end
  return (b...)
end
export flip_tuple_array


text:
"""
    growGenerator(gpExperimentInfo::GEPInfo, rng::Random.AbstractRNG)

documentation
"""
function growGenerator(gpExperimentInfo::GEPInfo, rng::Random.AbstractRNG)
    functionSet = gpExperimentInfo._functionSet
    terminalSet = deepcopy(gpExperimentInfo._terminalSet)
    headSize = gpExperimentInfo._head
    tailSize = gpExperimentInfo._tail
    repSize = headSize + tailSize

    representation = Array{Node}(undef, repSize)
    prob = length(terminalSet)/(length(terminalSet) + length(functionSet))

    for i=1:headSize
        if rand(rng) < prob
            representation[i] = rand(rng, terminalSet)
        else
            representation[i] = rand(rng, functionSet)
        end
    end

    for i=headSize+1:repSize
        representation[i] = rand(rng, terminalSet)
    end

    return GEPGenotype(representation)
end # function

"""
    fullGenerator(gpExperimentInfo::GEPInfo, rng::Random.AbstractRNG)

documentation
"""
function fullGenerator(gpExperimentInfo::GEPInfo, rng::Random.AbstractRNG)
    functionSet = gpExperimentInfo._functionSet
    terminalSet = deepcopy(gpExperimentInfo._terminalSet)
    headSize = gpExperimentInfo._head
    tailSize = gpExperimentInfo._tail
    repSize = headSize + tailSize

    representation = Array{Node}(undef, repSize)

    for i=1:headSize
        representation[i] = rand(rng, functionSet)
    end

    for i=headSize+1:repSize
        representation[i] = rand(rng, terminalSet)
    end

    return GEPGenotype(representation)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# Autogenerated wrapper script for PMIx_jll for aarch64-apple-darwin
export libpmix, pmix_info

using libevent_jll
using Hwloc_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("PMIx")
JLLWrappers.@declare_library_product(libpmix, "@rpath/libpmix.2.dylib")
JLLWrappers.@declare_executable_product(pmix_info)
function __init__()
    JLLWrappers.@generate_init_header(libevent_jll, Hwloc_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libpmix,
        "lib/libpmix.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        pmix_info,
        "bin/pmix_info",
    )

    JLLWrappers.@generate_init_footer()
    ENV["PMIX_PREFIX"] = artifact_dir

end  # __init__()


text:
function equals(obj::ComplexField, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_instance(::Type{ComplexField})
    return jcall(ComplexField, "getInstance", ComplexField, ())
end

function get_one(obj::ComplexField)
    return jcall(obj, "getOne", JComplex, ())
end

function get_runtime_class(obj::ComplexField)
    return jcall(obj, "getRuntimeClass", Class, ())
end

function get_zero(obj::ComplexField)
    return jcall(obj, "getZero", JComplex, ())
end

function hash_code(obj::ComplexField)
    return jcall(obj, "hashCode", jint, ())
end


text:
using .BenchmarkTools: @belapsed
using .ForwardDiff

"""
$(SIGNATURES)

Default chunk sizes to try for benchmarking. Fewer than `M`, always contains `1` and `N`.
"""
function heuristic_chunks(N, M = 20)
    step = max(N ÷ M, 1)
    Ns = 1:step:N
    if N ∉ Ns
        Ns = vcat(Ns, N)
    end
    Ns
end

"""
$(SIGNATURES)

Benchmark a log density problem with various chunk sizes using ForwardDiff.

`chunks`, which defaults to all possible chunk sizes, determines the chunks that are tried.

The function returns `chunk => time` pairs, where `time` is the benchmarked runtime in
seconds, as determined by `BenchmarkTools.@belapsed`. The gradient is evaluated at `x`
(defaults to zeros).

*Runtime may be long* because of tuned benchmarks, so when `markprogress == true` (the
default), dots are printed to mark progress.

This function is not exported, but part of the API.

*It is loaded conditionally when both `ForwardDiff` and `BenchmarkTools` are loaded.*
"""
function benchmark_ForwardDiff_chunks(ℓ;
                                      chunks = heuristic_chunks(dimension(ℓ), 20),
                                      markprogress = true,
                                      x = zeros(dimension(ℓ)))
    map(chunks) do chunk
        ∇ℓ = ADgradient(Val(:ForwardDiff), ℓ; chunk = ForwardDiff.Chunk(chunk))
        markprogress && print(".")
        chunk => @belapsed logdensity_and_gradient($(∇ℓ), $(x))
    end
end


text:
#=
* This is a very minimal implementation along the lines of properties from this:
  http://bioportal.bioontology.org/ontologies/FMA. Instead of defining a complex
  type hierarchy we need to use traits because this will invetiably evolve.

* It might be good adopt something like FilePathsBase.jl and parts of each "path"
  could be subject to these methods/traits.
=#
"""
    is_left(x) -> Bool

Returns `true` if `x` represent the left position.
"""
is_left(x::Symbol) = (x === :left) || (x === :L)

"""
    is_right(x) -> Bool

Returns `true` if `x` represent the right position.
"""
is_right(x::Symbol) = (x === :right) || (x === :R)

"""
    is_anterior(x) -> Bool

Returns `true` if `x` represent the anterior position.
"""
is_anterior(x::Symbol) = (x === :anterior) || (x === :A)

"""
    is_posterior(x) -> Bool

Returns `true` if `x` represent the posterior position.
"""
is_posterior(x::Symbol) = (x === :posterior) || (x === :P)

"""
    is_inferior(x) -> Bool

Returns `true` if `x` represent the inferior position.
"""
is_inferior(x::Symbol) = (x === :inferior) || (x === :I)

"""
    is_superior(x) -> Bool

Returns `true` if `x` represent the superior position.
"""
is_superior(x::Symbol) = (x === :superior) || (x === :S)

"""
    is_sagittal(x) -> Bool

Returns `true` if `x` represent the sagittal orientation.
"""
is_sagittal(x::Symbol) = is_left(x) || is_right(x) || (x === :sagittal)

"""
    is_coronal(x) -> Bool

Returns `true` if `x` represent the coronal orientation.
"""
is_coronal(x::Symbol) = is_anterior(x) || is_posterior(x) || (x === :coronal)

"""
    is_axial(x) -> Bool

Returns `true` if `x` represent the axial orientation.
"""
is_axial(x::Symbol) = is_superior(x) || is_inferior(x) || (x === :axial)

"""
    is_cortical(::T) -> Bool

Returns `true` if `T` represents a cortical region.
"""
is_cortical(x::Symbol) = x === :cortical

"""
    is_csf(::T) -> Bool

Returns `true` if `T` represents a region of corticospinal fluid (CSF).
"""
is_csf(x::Symbol) = x === :csf

"""
    is_white_matter(::T) -> Bool

Returns `true` if `T` represents white matter.
"""
is_white_matter(x::Symbol) = x === :white_matter

"""
    is_gyrus(::T) -> Bool

Returns `true` if `T` represents a gyrus.
"""
is_gyrus(x::Symbol) = x === :gyrus

"""
    is_sulcus(::T) -> Bool

Returns `true` if `T` represents a gyrus.
"""
is_sulcus(x::Symbol) = x === :sulcus

"""
    is_channel(::T) -> Bool

Returns `true` if `T` represents a gyrus.
"""
is_channel(x::Symbol) = x === :channel


text:
struct QueueDispatch
    queues::Dictionary{QueueFlag,Vector{Created{Queue,DeviceQueueInfo2}}}
    present_queue::Optional{Created{Queue,DeviceQueueInfo2}}
    """
    Build a `QueueDispatch` structure from a given device and configuration.

    !!! warning
        `device` must have been created with a consistent number of queues as requested in the provided queue configuration.
        It is highly recommended to have created the device with the result of `queue_infos(QueueDispatch, physical_device, config)`.
    """
    function QueueDispatch(device, config; surface = nothing)
        physical_device = handle(device).physical_device
        families = find_queue_family.(physical_device, collect(keys(config)))
        queues = dictionary(map(zip(keys(config), families)) do (property, family)
            property => map(0:config[property] - 1) do count
                info = DeviceQueueInfo2.(family, count)
                Created(get_device_queue_2(device, info), info)
            end
        end)
        present_queue = if !isnothing(surface)
            idx = findfirst(families) do family
                unwrap(get_physical_device_surface_support_khr(physical_device, family, surface))
            end
            if isnothing(idx)
                error("Could not find a queue that supports presentation on the provided surface.")
            else
                first(collect(queues)[idx])
            end
        else
            nothing
        end
        new(queues, present_queue)
    end
end

function queue_infos(::Type{QueueDispatch}, physical_device::PhysicalDevice, config)
    all(==(1), config) || error("Only one queue per property is currently supported")
    families = find_queue_family.(physical_device, collect(keys(config)))
    DeviceQueueCreateInfo.(families, ones.(collect(config)))
end

function submit(dispatch::QueueDispatch, properties::QueueFlag, submit_infos; fence = C_NULL)
    q = queue(dispatch, properties)
    unwrap(queue_submit_2_khr(q, submit_infos; fence))
    q
end

function queue(dispatch::QueueDispatch, properties::QueueFlag)
    if properties in keys(dispatch.queues)
        first(dispatch.queues[properties])
    else
        for props in keys(dispatch.queues)
            if properties in props
                return first(dispatch.queues[props])
            end
        end
        error("Could not find a queue matching with the required properties $properties.")
    end
end

function present(dispatch::QueueDispatch, present_info::PresentInfoKHR)
    queue = dispatch.present_queue
    if isnothing(queue)
        error("No presentation queue was specified for $dispatch")
    else
        queue_present_khr(queue, present_info)
    end
end

function queue_family_indices(dispatch::QueueDispatch; include_present = true)
    indices = map(dispatch.queues) do queues
        map(queues) do queue
            info(queue).queue_family_index
        end
    end
    indices = reduce(vcat, indices)
    if include_present && !isnothing(dispatch.present_queue)
        push!(indices, info(dispatch.present_queue).queue_family_index)
    end
    sort!(unique!(indices))
end


text:
using Turing

@model function two_model(x)
    # Hyper-parameters
    μ0 = 0.0
    σ0 = 1.0

    # Draw weights.
    π1 ~ Beta(1, 1)
    π2 = 1 - π1

    # Draw locations of the components.
    μ1 ~ Normal(μ0, σ0)
    μ2 ~ Normal(μ0, σ0)

    # Draw latent assignment.
    z ~ Categorical([π1, π2])

    # Draw observation from selected component.
    if z == 1
        x ~ Normal(μ1, 1.0)
    else
        x ~ Normal(μ2, 1.0)
    end
end

using Turing.RandomMeasures

# Concentration parameter.
α = 10.0

# Random measure, e.g. Dirichlet process.
rpm = DirichletProcess(α)

# Cluster assignments for each observation.
z = Vector{Int}()

# Maximum number of observations we observe.
Nmax = 500

for i in 1:Nmax
    # Number of observations per cluster.
    K = isempty(z) ? 0 : maximum(z)
    nk = Vector{Int}(map(k -> sum(z .== k), 1:K))

    # Draw new assignment.
    push!(z, rand(ChineseRestaurantProcess(rpm, nk)))
end

using Plots

# Plot the cluster assignments over time
@gif for i in 1:Nmax
    scatter(
        collect(1:i),
        z[1:i];
        markersize=2,
        xlabel="observation (i)",
        ylabel="cluster (k)",
        legend=false,
    )
end

@model function infiniteGMM(x)
    # Hyper-parameters, i.e. concentration parameter and parameters of H.
    α = 1.0
    μ0 = 0.0
    σ0 = 1.0

    # Define random measure, e.g. Dirichlet process.
    rpm = DirichletProcess(α)

    # Define the base distribution, i.e. expected value of the Dirichlet process.
    H = Normal(μ0, σ0)

    # Latent assignment.
    z = tzeros(Int, length(x))

    # Locations of the infinitely many clusters.
    μ = tzeros(Float64, 0)

    for i in 1:length(x)

        # Number of clusters.
        K = maximum(z)
        nk = Vector{Int}(map(k -> sum(z .== k), 1:K))

        # Draw the latent assignment.
        z[i] ~ ChineseRestaurantProcess(rpm, nk)

        # Create a new cluster?
        if z[i] > K
            push!(μ, 0.0)

            # Draw location of new cluster.
            μ[z[i]] ~ H
        end

        # Draw observation.
        x[i] ~ Normal(μ[z[i]], 1.0)
    end
end

using Plots, Random

# Generate some test data.
Random.seed!(1)
data = vcat(randn(10), randn(10) .- 5, randn(10) .+ 10)
data .-= mean(data)
data /= std(data);

# MCMC sampling
Random.seed!(2)
iterations = 1000
model_fun = infiniteGMM(data);
chain = sample(model_fun, SMC(), iterations);

# Extract the number of clusters for each sample of the Markov chain.
k = map(
    t -> length(unique(vec(chain[t, MCMCChains.namesingroup(chain, :z), :].value))),
    1:iterations,
);

# Visualize the number of clusters.
plot(k; xlabel="Iteration", ylabel="Number of clusters", label="Chain 1")

histogram(k; xlabel="Number of clusters", legend=false)

if isdefined(Main, :TuringTutorials)
    Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
end

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
# Classifies MNIST digits with a convolutional network.
# Writes out saved model to the file "mnist_conv.bson".
# Demonstrates basic model construction, training, saving,
# conditional early-exit, and learning rate scheduling.
#
# This model, while simple, should hit around 99% test
# accuracy after training for approximately 20 epochs.

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy
using Base.Iterators: repeated, partition
using Printf, BSON
using CUDAapi
if has_cuda()
    @info "CUDA is on"
    import CuArrays
    CuArrays.allowscalar(false)
end

# Load labels and images from Flux.Data.MNIST
@info("Loading data set")
train_labels = MNIST.labels()
train_imgs = MNIST.images()

# Bundle images together with labels and group into minibatchess
function make_minibatch(X, Y, idxs)
    X_batch = Array{Float32}(undef, size(X[1])..., 1, length(idxs))
    for i in 1:length(idxs)
        X_batch[:, :, :, i] = Float32.(X[idxs[i]])
    end
    Y_batch = onehotbatch(Y[idxs], 0:9)
    return (X_batch, Y_batch)
end
batch_size = 128
mb_idxs = partition(1:length(train_imgs), batch_size)
train_set = [make_minibatch(train_imgs, train_labels, i) for i in mb_idxs]

# Prepare test set as one giant minibatch:
test_imgs = MNIST.images(:test)
test_labels = MNIST.labels(:test)
test_set = make_minibatch(test_imgs, test_labels, 1:length(test_imgs))

# Define our model.  We will use a simple convolutional architecture with
# three iterations of Conv -> ReLU -> MaxPool, followed by a final Dense
# layer that feeds into a softmax probability output.
@info("Constructing model...")
model = Chain(
    # First convolution, operating upon a 28x28 image
    Conv((3, 3), 1=>16, pad=(1,1), relu),
    MaxPool((2,2)),

    # Second convolution, operating upon a 14x14 image
    Conv((3, 3), 16=>32, pad=(1,1), relu),
    MaxPool((2,2)),

    # Third convolution, operating upon a 7x7 image
    Conv((3, 3), 32=>32, pad=(1,1), relu),
    MaxPool((2,2)),

    # Reshape 3d tensor into a 2d one, at this point it should be (3, 3, 32, N)
    # which is where we get the 288 in the `Dense` layer below:
    x -> reshape(x, :, size(x, 4)),
    Dense(288, 10),

    # Finally, softmax to get nice probabilities
    softmax,
)

# Load model and datasets onto GPU, if enabled
train_set = gpu.(train_set)
test_set = gpu.(test_set)
model = gpu(model)

# Make sure our model is nicely precompiled before starting our training loop
model(train_set[1][1])

# We augment `x` a little bit here, adding in random noise.
augment(x) = x .+ gpu(0.1f0*randn(eltype(x), size(x)))

paramvec(m) = vcat(map(p->reshape(p, :), params(m))...)
anynan(x) = any(isnan.(x))

# `loss()` calculates the crossentropy loss between our prediction `y_hat`
# (calculated from `model(x)`) and the ground truth `y`.  We augment the data
# a bit, adding gaussian random noise to our image to make it more robust.
function loss(x, y)
    x̂ = augment(x)
    ŷ = model(x̂)
    return crossentropy(ŷ, y)
end
accuracy(x, y) = mean(onecold(cpu(model(x))) .== onecold(cpu(y)))

# Train our model with the given training set using the ADAM optimizer and
# printing out performance against the test set as we go.
opt = ADAM(0.001)

@info("Beginning training loop...")
best_acc = 0.0
last_improvement = 0
for epoch_idx in 1:100
    global best_acc, last_improvement
    # Train for a single epoch
    Flux.train!(loss, params(model), train_set, opt)

    if anynan(paramvec(model))
        @error "NaN params"
        break
    end

    # Calculate accuracy:
    acc = accuracy(test_set...)
    @info(@sprintf("[%d]: Test accuracy: %.4f", epoch_idx, acc))

    # If our accuracy is good enough, quit out.
    if acc >= 0.999
        @info(" -> Early-exiting: We reached our target accuracy of 99.9%")
        break
    end

    # If this is the best accuracy we've seen so far, save the model out
    if acc >= best_acc
        @info(" -> New best accuracy! Saving model out to mnist_conv.bson")
        BSON.@save joinpath(dirname(@__FILE__), "mnist_conv.bson") params=cpu.(params(model)) epoch_idx acc
        best_acc = acc
        last_improvement = epoch_idx
    end

    # If we haven't seen improvement in 5 epochs, drop our learning rate:
    if epoch_idx - last_improvement >= 5 && opt.eta > 1e-6
        opt.eta /= 10.0
        @warn(" -> Haven't improved in a while, dropping learning rate to $(opt.eta)!")

        # After dropping learning rate, give it a few epochs to improve
        last_improvement = epoch_idx
    end

    if epoch_idx - last_improvement >= 10
        @warn(" -> We're calling this converged.")
        break
    end
end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
using Shapefile, Plots
using Navigation
# using DDR2import

gr()

# function initialize_plot()
#     plot(xaxis=((-0.75,0.75), font(6)), yaxis=((-0.375,0.375), font(6)), legend = false,
#     grid = false, showaxis = true, ticks = false, xticks = -180.0:10.0:180.0,
#     yticks = -90.0:0.5:90.0, framestyle=:box)
# end

# function initialize_plot()
    # plot(xaxis=((-180.0,180.0), font(6)), yaxis=((-90.0,90.0), font(6)), legend = false,
    # grid = false, showaxis = false, ticks = false, xticks = -180.0:10.0:180.0,
    # yticks = -90.0:5.0:90.0, framestyle=:box)
# end

function initialize_plot()
    plot(legend = false, grid = false, showaxis = false, ticks = false,
    framestyle=:box)
    # plot(legend = false, grid = false, showaxis = false, ticks = false,
    # framestyle=:box, xticks = -180.0:10.0:180.0, yticks = -85.0:2.5:85.0)
    # plot!(showaxis = true, ticks = true, xticks = -180.0:10.0:180.0,
    # yticks = -85.0:2.5:85.0, xlim=(-40.0,40.0), ylim=(-20.0, 20.0))
end

function shapes2lonlat(shapes::Array{Union{Missing, Shapefile.Polygon},1})
    lon = Vector{Float64}()
    lat = Vector{Float64}()
    for shape in shapes
        for point in shape.points
            append!(lon, point.x)
            append!(lat, point.y)
        end
        append!(lon, NaN)
        append!(lat, NaN)
    end
    lon, lat
end

function sectors2lonlat(sectors::Array{Array{Array{Float64, 2},1},1})
    lon = Vector{Float64}()
    lat = Vector{Float64}()
    for sector in sectors
        for points in sector
                append!(lon, points[:,2])
                append!(lat, points[:,1])
        end
        append!(lon, NaN)
        append!(lat, NaN)
    end
    lon, lat
end

function grid2lonlat(grid_lon_deg::StepRangeLen, grid_lat_deg::StepRangeLen)
    lon = Vector{Float64}()
    lat = Vector{Float64}()
    lon_prev = NaN
    lat_prev = NaN
    fractions = 2
    function grid_loop(grid_lon::Float64, grid_lat::Float64,
        lon_prev::Float64, lat_prev::Float64, lon::Vector{Float64},
        lat::Vector{Float64}, fractions::Int64)
        if !isnan(lon_prev)
            for part in 1:(fractions - 1)
                fraction = 1.0 / fractions * part
                append!(lon, lon_prev + (grid_lon - lon_prev) * fraction)
                append!(lat, lat_prev + (grid_lat - lat_prev) * fraction)
            end
        end
        append!(lon, grid_lon)
        append!(lat, grid_lat)
        lon_prev = grid_lon
        lat_prev = grid_lat
        lon, lat, lon_prev, lat_prev
    end

    for grid_lat in grid_lat_deg
        for grid_lon in grid_lon_deg
            lon, lat, lon_prev, lat_prev = grid_loop(grid_lon, grid_lat,
            lon_prev, lat_prev, lon, lat, fractions)
        end
        append!(lon, NaN)
        append!(lat, NaN)
        lon_prev = NaN
        lat_prev = NaN
    end
    for grid_lon in grid_lon_deg
        for grid_lat in grid_lat_deg
            lon, lat, lon_prev, lat_prev = grid_loop(grid_lon, grid_lat,
            lon_prev, lat_prev, lon, lat, fractions)
        end
        append!(lon, NaN)
        append!(lat, NaN)
        lon_prev = NaN
        lat_prev = NaN
    end
    lon, lat
end

@recipe function f(shapes::Array{Union{Missing, Shapefile.Polygon},1}, projparam)
    lon, lat = shapes2lonlat(shapes)
    projection_short(ln::Float64, lt::Float64) = projection(ln, lt, projparam)
    points = vcat(projection_short.(lon, lat)...)
    points[:,1], points[:,2]
end

@recipe function f(grid_lon_deg::StepRangeLen, grid_lat_deg::StepRangeLen,
    projparam)
    lon, lat = grid2lonlat(grid_lon_deg, grid_lat_deg)
    projection_short(ln::Float64, lt::Float64) = projection(ln, lt, projparam)
    points = vcat(projection_short.(lon, lat)...)
    points[:,1], points[:,2]
end

@recipe function f(sectors::Vector{Vector{Array{Float64,2}}}, projparam)
    lon, lat = sectors2lonlat(sectors)
    projection_short(ln::Float64, lt::Float64) = projection(ln, lt, projparam)
    points = vcat(projection_short.(lon, lat)...)
    points[:,1], points[:,2]
end

function get_shapefile(shp::String)
    return Shapefile.shapes(Shapefile.Table(shp))
end

# function add_polygon!(shp::String, projparam=Equirectangular(); kwargs...)
#     plot!(get_shapefile(shp), projparam; kwargs...)
# end

function add_polygon!(shapes::Array{Union{Missing, Shapefile.Polygon},1}, projparam=Equirectangular(); kwargs...)
    plot!(shapes, projparam; kwargs...)
end

#TODO Fix are file reading
function add_polygon!(sectors::Vector{Vector{Array{Float64,2}}}, projparam=Equirectangular(); kwargs...)
    plot!(sectors, projparam; kwargs...)
end


text:
@testset "Utilities" begin
  @testset "spheredir" begin
    @test spheredir(90, 0) ≈ [1,0,0]
    @test spheredir(90,90) ≈ [0,1,0]
    @test spheredir(0,  0) ≈ [0,0,1]
  end

  @testset "aniso2distance" begin
    x, y = rand(2), rand(2)
    d = aniso2distance([1.0,1.0], [0.0])
    @test evaluate(d, [0,0], [1,0]) == 1
    @test evaluate(d, [0,0], [0,1]) == 1
    @test evaluate(d, x, y) ≈ euclidean(x, y)
  end
end


text:
export CuArray, CuVector, CuMatrix, CuVecOrMat, cu

mutable struct CuArray{T,N,P} <: AbstractGPUArray{T,N}
  ptr::CuPtr{T}
  dims::Dims{N}

  parent::P       # parent array (for, e.g., contiguous views keeping their parent alive)
  pooled::Bool    # is this memory backed by the memory pool?

  # for early freeing outside of the GC
  refcount::Int
  freed::Bool

  ctx::CuContext

  # constrain P
  CuArray{T,N,P}(ptr, dims, parent::Union{Nothing,CuArray}, pooled, ctx) where {T,N,P} =
    new(ptr, dims, parent, pooled, 0, false, ctx)
end

# primary array
function CuArray{T,N}(ptr::CuPtr{T}, dims::Dims{N}, pooled::Bool=true;
                      ctx=context()) where {T,N}
  self = CuArray{T,N,Nothing}(ptr, dims, nothing, pooled, ctx)
  retain(self)
  finalizer(unsafe_free!, self)
  return self
end

# derived array (e.g. view, reinterpret, ...)
function CuArray{T,N}(ptr::CuPtr{T}, dims::Dims{N}, parent::P) where {T,N,P<:CuArray}
  self = CuArray{T,N,P}(ptr, dims, parent, parent.pooled, parent.ctx)
  retain(self)
  retain(parent)
  finalizer(unsafe_free!, self)
  return self
end

function unsafe_free!(xs::CuArray)
  # this call should only have an effect once, becuase both the user and the GC can call it
  xs.freed && return
  _unsafe_free!(xs)
  xs.freed = true
  return
end

function _unsafe_free!(xs::CuArray)
  @assert xs.refcount >= 0
  if release(xs)
    if xs.parent === nothing
      # primary array with all references gone
      if xs.pooled && isvalid(xs.ctx)
        free(convert(CuPtr{Nothing}, pointer(xs)))
      end
    else
      # derived object
      _unsafe_free!(xs.parent)
    end

    # the object is dead, so we can also wipe the pointer
    xs.ptr = CU_NULL
  end

  return
end

@inline function retain(a::CuArray)
  a.refcount += 1
  return
end

@inline function release(a::CuArray)
  a.refcount -= 1
  return a.refcount == 0
end

Base.parent(A::CuArray{<:Any,<:Any,Nothing})     = A
Base.parent(A::CuArray{<:Any,<:Any,P}) where {P} = A.parent

Base.dataids(A::CuArray{<:Any,<:Any,Nothing})     = (UInt(pointer(A)),)
Base.dataids(A::CuArray{<:Any,<:Any,P}) where {P} = (Base.dataids(A.parent)..., UInt(pointer(A)),)

# TODO: implement array alias detection from https://github.com/JuliaLang/julia/pull/25890

Base.unaliascopy(A::CuArray{<:Any,<:Any,Nothing}) = copy(A)
function Base.unaliascopy(A::CuArray{<:Any,<:Any,P}) where {P}
  offset = pointer(A) - pointer(A.parent)
  new_parent = Base.unaliascopy(A.parent)
  typeof(A)(pointer(new_parent) + offset, A.dims, new_parent, A.pooled, A.ctx)
end

## convenience constructors

# discard the P typevar
#
# P is just used to specialize the parent field, and does not actually affect the object,
# so we can safely discard this information when creating similar objects (`typeof(A)(...)`)
CuArray{T,N,P}(args...) where {T,N,P} = CuArray{T,N}(args...)

CuVector{T} = CuArray{T,1}
CuMatrix{T} = CuArray{T,2}
CuVecOrMat{T} = Union{CuVector{T},CuMatrix{T}}

# type and dimensionality specified, accepting dims as tuples of Ints
function CuArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N}
  Base.isbitsunion(T) && error("CuArray does not yet support union bits types")
  Base.isbitstype(T)  || error("CuArray only supports bits types") # allocatedinline on 1.3+
  ptr = alloc(prod(dims) * sizeof(T))
  CuArray{T,N}(convert(CuPtr{T}, ptr), dims)
end

# type and dimensionality specified, accepting dims as series of Ints
CuArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = CuArray{T,N}(undef, dims)

# type but not dimensionality specified
CuArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = CuArray{T,N}(undef, dims)
CuArray{T}(::UndefInitializer, dims::Integer...) where {T} =
  CuArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# empty vector constructor
CuArray{T,1}() where {T} = CuArray{T,1}(undef, 0)

# do-block constructors
for (ctor, tvars) in (:CuArray => (), :(CuArray{T}) => (:T,), :(CuArray{T,N}) => (:T, :N))
  @eval begin
    function $ctor(f::Function, args...) where {$(tvars...)}
      xs = $ctor(args...)
      try
        f(xs)
      finally
        unsafe_free!(xs)
      end
    end
  end
end

Base.similar(a::CuArray{T,N}) where {T,N} = CuArray{T,N}(undef, size(a))
Base.similar(a::CuArray{T}, dims::Base.Dims{N}) where {T,N} = CuArray{T,N}(undef, dims)
Base.similar(a::CuArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = CuArray{T,N}(undef, dims)

"""
  unsafe_wrap(::CuArray, ptr::CuPtr{T}, dims; own=false, ctx=context())

Wrap a `CuArray` object around the data at the address given by `ptr`. The pointer
element type `T` determines the array element type. `dims` is either an integer (for a 1d
array) or a tuple of the array dimensions. `own` optionally specified whether Julia should
take ownership of the memory, calling `cudaFree` when the array is no longer referenced. The
`ctx` argument determines the CUDA context where the data is allocated in.
"""
function Base.unsafe_wrap(::Union{Type{CuArray},Type{CuArray{T}},Type{CuArray{T,N}}},
                          p::CuPtr{T}, dims::NTuple{N,Int};
                          own::Bool=false, ctx::CuContext=context()) where {T,N}
  xs = CuArray{T, length(dims)}(p, dims, false; ctx=ctx)
  if own
    base = convert(CuPtr{Cvoid}, p)
    buf = Mem.DeviceBuffer(base, prod(dims) * sizeof(T), ctx)
    finalizer(xs) do obj
      if isvalid(obj.ctx)
        Mem.free(buf)
      end
    end
  end
  return xs
end

function Base.unsafe_wrap(Atype::Union{Type{CuArray},Type{CuArray{T}},Type{CuArray{T,1}}},
                          p::CuPtr{T}, dim::Integer;
                          own::Bool=false, ctx::CuContext=context()) where {T}
  unsafe_wrap(Atype, p, (dim,); own=own, ctx=ctx)
end

Base.unsafe_wrap(T::Type{<:CuArray}, ::Ptr, dims::NTuple{N,Int}; kwargs...) where {N} =
  throw(ArgumentError("cannot wrap a CPU pointer with a $T"))

## array interface

Base.elsize(::Type{<:CuArray{T}}) where {T} = sizeof(T)

Base.size(x::CuArray) = x.dims
Base.sizeof(x::CuArray) = Base.elsize(x) * length(x)

Base.pointer(x::CuArray) = x.ptr
Base.pointer(x::CuArray, i::Integer) = x.ptr + (i-1) * Base.elsize(x)

## interop with other arrays

@inline function CuArray{T,N}(xs::AbstractArray{T,N}) where {T,N}
  A = CuArray{T,N}(undef, size(xs))
  copyto!(A, xs)
  return A
end

# FIXME: `map(T, xs)`, https://github.com/FluxML/Flux.jl/issues/958
CuArray{T,N}(xs::AbstractArray{S,N}) where {T,N,S} = CuArray{T,N}((x -> T(x)).(xs))

# underspecified constructors
CuArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = CuArray{T,N}(xs)
(::Type{CuArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = CuArray{S,N}(x)
CuArray(A::AbstractArray{T,N}) where {T,N} = CuArray{T,N}(A)

# idempotency
CuArray{T,N}(xs::CuArray{T,N}) where {T,N} = xs

## conversions

Base.convert(::Type{T}, x::T) where T <: CuArray = x

function Base._reshape(parent::CuArray, dims::Dims)
  n = length(parent)
  prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
  return CuArray{eltype(parent),length(dims)}(pointer(parent), dims, parent)
end
function Base._reshape(parent::CuArray{T,1}, dims::Tuple{Int}) where T
  n = length(parent)
  prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
  return parent
end

## interop with C libraries

Base.unsafe_convert(::Type{Ptr{T}}, x::CuArray{T}) where {T} = throw(ArgumentError("cannot take the CPU address of a $(typeof(x))"))
Base.unsafe_convert(::Type{Ptr{S}}, x::CuArray{T}) where {S,T} = throw(ArgumentError("cannot take the CPU address of a $(typeof(x))"))

Base.unsafe_convert(::Type{CuPtr{T}}, x::CuArray{T}) where {T} = pointer(x)
Base.unsafe_convert(::Type{CuPtr{S}}, x::CuArray{T}) where {S,T} = convert(CuPtr{S}, Base.unsafe_convert(CuPtr{T}, x))

## interop with device arrays

function Base.convert(::Type{CuDeviceArray{T,N,AS.Global}}, a::CuArray{T,N}) where {T,N}
  CuDeviceArray{T,N,AS.Global}(a.dims, DevicePtr{T,AS.Global}(pointer(a)))
end

Adapt.adapt_storage(::Adaptor, xs::CuArray{T,N}) where {T,N} =
  convert(CuDeviceArray{T,N,AS.Global}, xs)

## interop with CPU arrays

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{CuArray}, xs::AbstractArray) =
  isbits(xs) ? xs : convert(CuArray, xs)

# aggressively convert arrays of floats to float32
Adapt.adapt_storage(::Type{CuArray}, xs::AbstractArray{<:AbstractFloat}) =
  isbits(xs) ? xs : convert(CuArray{Float32}, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:CuArray{T}}, xs::AbstractArray) where {T} =
  isbits(xs) ? xs : convert(CuArray{T}, xs)

Adapt.adapt_storage(::Type{Array}, xs::CuArray) = convert(Array, xs)

Base.collect(x::CuArray{T,N}) where {T,N} = copyto!(Array{T,N}(undef, size(x)), x)

function Base.copyto!(dest::CuArray{T}, doffs::Integer, src::Array{T}, soffs::Integer,
                      n::Integer) where T
  @assert !dest.freed "Use of freed memory"
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

function Base.copyto!(dest::Array{T}, doffs::Integer, src::CuArray{T}, soffs::Integer,
                      n::Integer) where T
  @assert !src.freed "Use of freed memory"
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

function Base.copyto!(dest::CuArray{T}, doffs::Integer, src::CuArray{T}, soffs::Integer,
                      n::Integer) where T
  @assert !dest.freed && !src.freed "Use of freed memory"
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

function Base.unsafe_copyto!(dest::CuArray{T}, doffs, src::Array{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.unsafe_copyto!(dest::Array{T}, doffs, src::CuArray{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.unsafe_copyto!(dest::CuArray{T}, doffs, src::CuArray{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.deepcopy_internal(x::CuArray, dict::IdDict)
  haskey(dict, x) && return dict[x]::typeof(x)
  return dict[x] = copy(x)
end

## utilities

cu(xs) = adapt(CuArray, xs)
Base.getindex(::typeof(cu), xs...) = CuArray([xs...])

zeros(T::Type, dims...) = fill!(CuArray{T}(undef, dims...), 0)
ones(T::Type, dims...) = fill!(CuArray{T}(undef, dims...), 1)
zeros(dims...) = zeros(Float32, dims...)
ones(dims...) = ones(Float32, dims...)
fill(v, dims...) = fill!(CuArray{typeof(v)}(undef, dims...), v)
fill(v, dims::Dims) = fill!(CuArray{typeof(v)}(undef, dims...), v)

# optimized implementation of `fill!` for types that are directly supported by memset
const MemsetTypes = Dict(1=>UInt8, 2=>UInt16, 4=>UInt32)
const MemsetCompatTypes = Union{UInt8, Int8,
                                UInt16, Int16, Float16,
                                UInt32, Int32, Float32}
function Base.fill!(A::CuArray{T}, x) where T <: MemsetCompatTypes
  U = MemsetTypes[sizeof(T)]
  y = reinterpret(U, convert(T, x))
  Mem.set!(convert(CuPtr{U}, pointer(A)), y, length(A))
  A
end

## reversing

# the kernel works by treating the array as 1d. after reversing by dimension x an element at
# pos [i1, i2, i3, ... , i{x},            ..., i{n}] will be at
# pos [i1, i2, i3, ... , d{x} - i{x} + 1, ..., i{n}] where d{x} is the size of dimension x

# out-of-place version, copying a single value per thread from input to output
function _reverse(input::CuArray{T, N}, output::CuArray{T, N}; dims::Integer=1) where {T, N}
    @assert size(input) == size(output)
    shape = [size(input)...]
    numelemsinprevdims = prod(shape[1:dims-1])
    numelemsincurrdim = shape[dims]

    function kernel(input::CuDeviceArray{T, N}, output::CuDeviceArray{T, N}) where {T, N}
        offset_in = blockDim().x * (blockIdx().x - 1)

        index_in = offset_in + threadIdx().x

        if index_in <= length(input)
            element = @inbounds input[index_in]

            # the index of an element in the original array along dimension that we will flip
            #assume(numelemsinprevdims > 0)
            #assume(numelemsincurrdim > 0)
            ik = ((cld(index_in, numelemsinprevdims) - 1) % numelemsincurrdim) + 1

            index_out = index_in + (numelemsincurrdim - 2ik + 1) * numelemsinprevdims

            @inbounds output[index_out] = element
        end

        return
    end

    nthreads = 256
    nblocks = cld(prod(shape), nthreads)
    shmem = nthreads * sizeof(T)

    @cuda threads=nthreads blocks=nblocks kernel(input, output)
end

# in-place version, swapping two elements on half the number of threads
function _reverse(data::CuArray{T, N}; dims::Integer=1) where {T, N}
    shape = [size(data)...]
    numelemsinprevdims = prod(shape[1:dims-1])
    numelemsincurrdim = shape[dims]

    function kernel(data::CuDeviceArray{T, N}) where {T, N}
        offset_in = blockDim().x * (blockIdx().x - 1)

        index_in = offset_in + threadIdx().x

        # the index of an element in the original array along dimension that we will flip
        #assume(numelemsinprevdims > 0)
        #assume(numelemsincurrdim > 0)
        ik = ((cld(index_in, numelemsinprevdims) - 1) % numelemsincurrdim) + 1

        index_out = index_in + (numelemsincurrdim - 2ik + 1) * numelemsinprevdims

        if index_in <= length(data) && index_in < index_out
            @inbounds begin
                temp = data[index_out]
                data[index_out] = data[index_in]
                data[index_in] = temp
            end
        end

        return
    end

    # NOTE: we launch twice the number of threads, which is wasteful, but the ND index
    #       calculations don't allow using only the first half of the threads
    #       (e.g. [1 2 3; 4 5 6] where threads 1 and 2 swap respectively (1,2) and (2,1)).

    nthreads = 256
    nblocks = cld(prod(shape), nthreads)
    shmem = nthreads * sizeof(T)

    @cuda threads=nthreads blocks=nblocks kernel(data)
end

# n-dimensional API

# in-place
function Base.reverse!(data::CuArray{T, N}; dims::Integer) where {T, N}
    if !(1 ≤ dims ≤ length(size(data)))
      ArgumentError("dimension $dims is not 1 ≤ $dims ≤ $length(size(input))")
    end

    _reverse(data; dims=dims)

    return data
end

# out-of-place
function Base.reverse(input::CuArray{T, N}; dims::Integer) where {T, N}
    if !(1 ≤ dims ≤ length(size(input)))
      ArgumentError("dimension $dims is not 1 ≤ $dims ≤ $length(size(input))")
    end

    output = similar(input)
    _reverse(input, output; dims=dims)

    return output
end

# 1-dimensional API

# in-place
function Base.reverse!(data::CuVector{T}, start=1, stop=length(data)) where {T}
    _reverse(view(data, start:stop))
    return data
end

# out-of-place
function Base.reverse(input::CuVector{T}, start=1, stop=length(input)) where {T}
    output = similar(input)

    start > 1 && copyto!(output, 1, input, 1, start-1)
    _reverse(view(input, start:stop), view(output, start:stop))
    stop < length(input) && copyto!(output, stop+1, input, stop+1)

    return output
end

## resizing

"""
  resize!(a::CuVector, n::Int)

Resize `a` to contain `n` elements. If `n` is smaller than the current collection length,
the first `n` elements will be retained. If `n` is larger, the new elements are not
guaranteed to be initialized.

Several restrictions apply to which types of `CuArray`s can be resized:

- the array should be backed by the memory pool, and not have been constructed with `unsafe_wrap`
- the array cannot be derived (view, reshape) from another array
- the array cannot have any derived arrays itself

"""
function Base.resize!(A::CuVector{T}, n::Int) where T
  A.parent === nothing || error("cannot resize derived CuArray")
  A.refcount == 1 || error("cannot resize shared CuArray")
  A.pooled || error("cannot resize wrapped CuArray")

  ptr = convert(CuPtr{T}, alloc(n * sizeof(T)))
  m = Base.min(length(A), n)
  unsafe_copyto!(ptr, pointer(A), m)

  free(convert(CuPtr{Nothing}, pointer(A)))
  A.dims = (n,)
  A.ptr = ptr

  A
end


text:
using Documenter, TableShowUtils

makedocs(
	modules = [TableShowUtils],
	sitename = "TableShowUtils.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/TableShowUtils.jl.git"
)


text:
using NGWP, Plots, LightGraphs, Random
using PrettyTables

nsim = 500

function generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q; edge_length = 1)
    Random.seed!(1234)
    ρ = zeros(nsim)
    for i = 1:nsim
        p = rand(N); p ./= norm(p, 1)
        q = rand(N); q ./= norm(q, 1)
        W1 = ROT_Distance(p, q, Q; edge_length = edge_length)
        K = K_functional(p, q, 𝚽, ∇𝚽, 𝛌; length = edge_length)[1]
        ρ[i] = K / W1
    end
    return ρ
end

function display_basic_stats(ρs)
    header = ["min" "max" "mean" "std"]
    basic_stats = zeros(0, 4)
    for ρ in ρs
        basic_stats = vcat(
            basic_stats,
            round.([minimum(ρ) maximum(ρ) mean(ρ) std(ρ)]; digits = 4)
            )
    end
    pretty_table(basic_stats, header)
end

function ROT_TSD_ratio_histogram(ρ)
    plt = histogram(ρ, grid = false, legend = false, c = :teal,
              xlims = [minimum(ρ) - std(ρ), maximum(ρ) + std(ρ)],
              xlab = "ρ")
    return plt
end


text:
# Script to run a test case for understanding how the probability of success scales with active learning #

# Add paths
loc_gm_package = "../src/"
if !(loc_gm_package in LOAD_PATH)
    push!(LOAD_PATH, "../src/")
end

# Required packages
using GraphicalModelLearning
using LightGraphs
using Cairo
using Junet

# Function definitions
include("beta_scaling.jl")

# parameters
FLAG_create_struct_gm = true
FLAG_regular_random_gm = false
FLAG_lattice_gm = true
FLAG_weak_impurity = false
sampling_regime = M_regime()
learning_algo = CD(1e-12)

if FLAG_lattice_gm
    # text file to save adjacency matrix in
    file_adj_matrix_gm = "adj_matrix_ferro_lattice_gm_A_ME.txt"
    # name of picture to save graphical model in
    file_plot_gm = "ferro_lattice_gm_A_ME.eps"
    # File to save final results
    file_M_opt_gm = "M_opt_Ferro_Lattice_A_ME.txt"
elseif FLAG_regular_random_gm
    # text file to save adjacency matrix in
    file_adj_matrix_gm = "adj_matrix_spin_glass_gm.txt"
    # name of picture to save graphical model in
    file_plot_gm = "spin_glass_gm3.eps"
    # File to save final results
    file_M_opt_gm = "M_opt_M_Regime_Regular_Random_B.txt"
end

# Simple test
N = 16
d = 3
α = 0.4
β_array = [0.4+0.1*i for i=0:46]

# Create and plot the initial graphical graphical model
β = copy(β_array[4])
if FLAG_create_struct_gm
    if FLAG_regular_random_gm
        adj_matrix, struct_adj_matrix = spin_glass_random_regular_model(N,d,α,β)
        gm = Junet.Graph(struct_adj_matrix,directed=false)

        h = plot(gm, size=(300,300), layout=layout_fruchterman_reingold(gm), node_color="white", node_border_color="black",
            node_border_width=2, edge_color=[edge_color_gm(i, adj_matrix, α, β) for i = 1:edgecount(gm)],
            edge_width=2, format=:eps, file=file_plot_gm)
    elseif FLAG_lattice_gm
        adj_matrix, struct_adj_matrix = ferro_lattice(N,α,β,FLAG_weak_impurity)
        m = n = Int(sqrt(N))
        draw_periodic_lattice(adj_matrix,m,n,α,β,file_plot_gm)
        open(file_adj_matrix_gm, "w") do io
            writedlm(io, adj_matrix)
        end
    end
else
    adj_matrix = readdlm(file_adj_matrix_gm)
end

## Start the complexity studies
τ = α/2
L = 45
M_guess = 1000
M_factor = 0.05

M_opt = Array{Int64,1}(length(β_array))

if FLAG_lattice_gm
    for i = 1:length(β_array)
        # define β
        β = copy(β_array[i])
        println(β)

        # Create the adjacency matrix
        adj_matrix, struct_adj_matrix = ferro_lattice(N,α,β,FLAG_weak_impurity)

        # Get the optimum number of samples
        M_opt[i] = get_M_opt_glauber_dynamics(adj_matrix, RISE(), learning_algo, sampling_regime, τ, L, M_guess, M_factor)

        M_guess = Int(floor(0.8*copy(M_opt[i])))
    end
end

writedlm(file_M_opt_gm,M_opt)


text:
using DifferentialEquations
using Plots

# State Vector
#  X = [u,v,w,p,q,r,phi,theta,psi,x,y,z]

# time span
t_begin = 0
t_end = 10.0

tspan = (t_begin,t_end)

##########################################
# initial conditions
##########################################
u0 = 100
v0 = 0
w0 = 3

p0 = 0
q0 = 0
r0 = 0

phi0 = 0
theta0 = deg2rad(2)
psi0 = 0

x0 = 0
y0 = 0
z0 = -500

###############
# Let's get it running with CVLM in the loop - then, we will need to come up with a trim routine here for steady level flight. I've done this before, so I can add it.
##############

# inital condition vector [tranlational velocities | angular velocities | euler angles | cartesian position ]
y0 = [u0;v0;w0;p0;q0;r0;phi0;theta0;psi0;x0;y0;z0]

#############################################
# parameter definitions -
## **** Rohit - Can you make these into a Tuple? I tried unsuccessfully *******
#############################################
m = 10
g = 9.81
rho = 1.225
T = 0.2 # thrust (will need to determine from induced drag at trim?)
# Inertia constants for calculating GAMMA, C1-C9
Ix = 1
Iy = 1
Iz = 1
Ixz = 0.2
# constants to deal with inertial coupling
GAMMA = Ix*Iz-Ixz^2
c1 = [(Iy-Iz)*Iz-Ixz^2]/GAMMA
c2 = [(Ix-Iy+Iz)*Ixz]/GAMMA
c3 =  Iz/GAMMA
c4 = Ixz/GAMMA
c5 = (Iz-Ix)/Iy
c6 = Ixz/Iy
c7 = 1/Iy
c8 = [(Ix-Iy)*Ix-Ixz^2]/GAMMA
c9 = Ix/GAMMA
# *** Rohit -  add matrices from CVLM here
# all we need to do is define Px, Py, Pz, Qx, Qy, Qz as parameter to ODE solver
Px = zeros(6,6)
Py = zeros(6,6)
Pz = zeros(6,6)
Qx = zeros(6,6)
Qy = zeros(6,6)
Qz = zeros(6,6)

# rotating mass along longitudinal axis (spinning propellor or turbomachinery) - these will probably always be zero for using
Ip = 0  # rotating inertia
OmegaP = 0  # rpm of rotating machinery

##################################################
# Define parameter vector (would be more convienient as a tuple)
##################################################
p = (m,g,rho,c1,c2,c3,c4,c5,c6,c7,c8,c9,T,Ip,OmegaP)

## ODE Function - For some reason I can't save this as a separate file.... probably Julia syntax
function calc_dx(dx,x,p,t)

    #X = [u,v,w,p,q,r,phi,theta,psi,x,y,z]

    # parameters:
    m,g,rho,c1,c2,c3,c4,c5,c6,c7,c8,T,Ip,OmegaP = p
    # current state values:
    u,v,w,p,q,r,phi,theta,psi,x,y,z = x

    # Aerodynamic Force/Moment Calculation (without controls {for now})
    X = [u;v;w;p;q;r]

    Fx = rho*X'*Px*X
    Fy = rho*X'*Py*X
    Fz = rho*X'*Pz*X

    Mx = rho*X'*Qx*X
    Mx = rho*X'*Qy*X
    Mz = rho*X'*Qz*X

    # state derivative calculations

    dx[1] = r*v-q*w+Fx/m-g*sin(theta)+T/m
    dx[2] = p*w-r*u+Fy/m+g*cos(theta)*sin(phi)  # Force Equations
    dx[3] = q*u-p*v+Fz/m+g*cos(theta)*cos(phi)

    dx[4] = (c1*r+c2*p-c4*Ip*OmegaP)*q.+c3*Mx.+c4*Mz
    dx[5] = (c5*p+c7*Ip*OmegaP)*r-c6*(p^2-r^2)+c7*My   # moment equations
    dx[6] = (c8*p-c2*r-c9*Ip*OmegaP)*q+c9*Mz+c4*Mx

    dx[7] = p+tan(theta)*(q*sin(phi)+r*cos(phi))
    dx[8] = q*cos(phi)-r*sin(phi)
    dx[9] = (q*sin(phi)+r*cos(phi))/cos(theta)

    dx[10:12] = [1 0 0; 0 cos(phi) sin(phi);0 -sin(phi) cos(phi)]*[cos(theta) 0 -sin(theta);0 1 0; sin(theta) 0 cos(theta)]*[cos(phi) sin(phi) 0;-sin(phi) cos(phi) 0; 0 0 1]*x[1:3]

end

prob = ODEProblem(calc_dx,y0,tspan,p)
sol = solve(prob)

plot(sol)

max_stars_repo_name: rohitvuppala/CVLM

max_issues_repo_name: rohitvuppala/CVLM

max_forks_repo_name: rohitvuppala/CVLM


text:
function factorize_wrapper!(ips::InteriorPointSolver)
    @trace(ips.logger,"Factorization started.")
    build_kkt!(ips.kkt)
    ips.cnt.linear_solver_time += @elapsed factorize!(ips.linear_solver)
end

function solve_refine_wrapper!(
    ips::InteriorPointSolver,
    x::AbstractKKTVector,
    b::AbstractKKTVector,
)
    cnt = ips.cnt
    @trace(ips.logger,"Iterative solution started.")
    fixed_variable_treatment_vec!(full(b), ips.ind_fixed)

    cnt.linear_solver_time += @elapsed begin
        result = solve_refine!(x, ips.iterator, b)
    end

    if result == :Solved
        solve_status =  true
    else
        if improve!(ips.linear_solver)
            cnt.linear_solver_time += @elapsed begin
                factorize!(ips.linear_solver)
                ret = solve_refine!(x, ips.iterator, b)
                solve_status = (ret == :Solved)
            end
        else
            solve_status = false
        end
    end
    fixed_variable_treatment_vec!(full(x), ips.ind_fixed)
    return solve_status
end

function solve_refine_wrapper!(
    ips::InteriorPointSolver{<:DenseCondensedKKTSystem},
    x::AbstractKKTVector,
    b::AbstractKKTVector,
)
    cnt = ips.cnt
    @trace(ips.logger,"Iterative solution started.")
    fixed_variable_treatment_vec!(full(b), ips.ind_fixed)

    kkt = ips.kkt

    n = num_variables(kkt)
    n_eq, ns = kkt.n_eq, kkt.n_ineq
    n_condensed = n + n_eq

    # load buffers
    b_c = view(full(ips._w1), 1:n_condensed)
    x_c = view(full(ips._w2), 1:n_condensed)
    jv_x = view(full(ips._w3), 1:ns) # for jprod
    jv_t = primal(ips._w4)             # for jtprod
    v_c = dual(ips._w4)

    Σs = get_slack_regularization(kkt)
    α = get_scaling_inequalities(kkt)

    # Decompose right hand side
    bx = view(full(b), 1:n)
    bs = view(full(b), n+1:n+ns)
    by = view(full(b), kkt.ind_eq_shifted)
    bz = view(full(b), kkt.ind_ineq_shifted)

    # Decompose results
    xx = view(full(x), 1:n)
    xs = view(full(x), n+1:n+ns)
    xy = view(full(x), kkt.ind_eq_shifted)
    xz = view(full(x), kkt.ind_ineq_shifted)

    v_c .= 0.0
    v_c[kkt.ind_ineq] .= (Σs .* bz .+ α .* bs) ./ α.^2
    jtprod!(jv_t, kkt, v_c)
    # init right-hand-side
    b_c[1:n] .= bx .+ jv_t[1:n]
    b_c[1+n:n+n_eq] .= by

    cnt.linear_solver_time += @elapsed (result = solve_refine!(x_c, ips.iterator, b_c))
    solve_status = (result == :Solved)

    # Expand solution
    xx .= x_c[1:n]
    xy .= x_c[1+n:end]
    jprod_ineq!(jv_x, kkt, xx)
    xz .= sqrt.(Σs) ./ α .* jv_x .- Σs .* bz ./ α.^2 .- bs ./ α
    xs .= (bs .+ α .* xz) ./ Σs

    fixed_variable_treatment_vec!(full(x), ips.ind_fixed)
    return solve_status
end