text:
Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
[@LorenaABarba](https://twitter.com/LorenaABarba)

12 steps to Navier–Stokes
=====
***

In the previous step, we solved the [2D Burgers' equation](./10_Step_8.ipynb): an important equation in the study of fluid mechanics because it contains the full convective nonlinearity of the flow equations. With that exercise, we also build the experience to incrementatlly code a Navier–Stokes solver.

In the next two steps, we will solve Laplace and then Poisson equation. We will then put it all together!

Step 9: 2D Laplace Equation
----
***

Here is Laplace's equation in 2D:

$$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0$$

We know how to discretize a 2nd order derivative. But think about this for a minute — Laplace's equation has the features typical of diffusion phenomena. For this reason, it has to be discretized with *central differences*, so that the discretization is consistent with the physics we want to simulate.

The discretized equation is:

$$\frac{p_{i+1, j}^n - 2p_{i,j}^n + p_{i-1,j}^n}{\Delta x^2} + \frac{p_{i,j+1}^n - 2p_{i,j}^n + p_{i, j-1}^n}{\Delta y^2} = 0$$

Notice that the Laplace Equation does not have a time dependence — there is no $p^{n+1}$.  Instead of tracking a wave through time (like in the previous steps), the Laplace equation calculates the equilibrium state of a system under the supplied boundary conditions.

If you have taken coursework in Heat Transfer, you will recognize the Laplace Equation as the steady-state heat equation.

Instead of calculating where the system will be at some time $t$, we will iteratively solve for $p_{i,j}^n$ until it meets a condition that we specify.  The system will reach equilibrium only as the number of iterations tends to $\infty$, but we can approximate the equilibrium state by iterating until the change between one iteration and the next is *very* small.

Let's rearrange the discretized equation, solving for $p_{i,j}^n$:

$$p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}$$

Using second-order central-difference schemes in both directions is the most widely applied method for the Laplace operator. It is also known as the **five-point difference operator**, alluding to its stencil.

We are going to solve Laplace's equation numerically by assuming an initial state of $p=0$ everywhere. Then we add boundary conditions as follows:

$p=0$ at $x=0$

$p=y$ at $x=2$

$\frac{\partial p}{\partial y}=0$ at $y=0, \ 1$

Under these conditions, there is an analytical solution for Laplace's equation:

$$p(x,y)=\frac{x}{4}-4\sum_{n=1,odd}^{\infty}\frac{1}{(n\pi)^2\sinh2n\pi}\sinh n\pi x\cos n\pi y$$

##### Exercise

Write your own code to solve Poisson's equation using loops, in the style of coding used in our first lessons. Then, consider the demonstration of how to write it using functions (below) and modify your code in that style. Can you think of reasons to abandon the old style and adopt modular coding?

Other tips:

+ Visualize each step of the iterative process
+ Think about what the boundary conditions are doing
+ Think about what the PDE is doing

### Using functions

Remember the lesson on writing [functions with Python](./11_Defining_Function_in_Python.ipynb)? We will use that style of code in this exercise.

We're going to define two functions: one that plots our data in a 3D projection plot and the other that iterates to solve for $p$ until the change in the [L1 Norm][1] of $p$ is less than a specified value.

[1]: http://en.wikipedia.org/wiki/Norm_(mathematics)#Taxicab_norm_or_Manhattan_norm

```python
import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
```

```python
def plot2D(x, y, p):
    fig = pyplot.figure(figsize=(11, 7), dpi=100)
    ax = fig.gca(projection='3d')
    X, Y = numpy.meshgrid(x, y)
    surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,
            linewidth=0, antialiased=False)
    ax.set_xlim(0, 2)
    ax.set_ylim(0, 1)
    ax.view_init(30, 225)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

```

The function `plot2D` takes three arguments, an x-vector, a y-vector and our p matrix.  Given these three values, it produces a 3D projection plot, sets the plot limits and gives us a nice viewing angle.

$$p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}$$

```python
def laplace2d(p, y, dx, dy, l1norm_target):
    l1norm = 1
    pn = numpy.empty_like(p)

    while l1norm > l1norm_target:
        pn = p.copy()
        p[1:-1, 1:-1] = ((dy**2 * (pn[1:-1, 2:] + pn[1:-1, 0:-2]) +
                         dx**2 * (pn[2:, 1:-1] + pn[0:-2, 1:-1])) /
                        (2 * (dx**2 + dy**2)))

        p[:, 0] = 0  # p = 0 @ x = 0
        p[:, -1] = y  # p = y @ x = 2
        p[0, :] = p[1, :]  # dp/dy = 0 @ y = 0
        p[-1, :] = p[-2, :]  # dp/dy = 0 @ y = 1
        l1norm = (numpy.sum(numpy.abs(p[:]) - numpy.abs(pn[:])) /
                numpy.sum(numpy.abs(pn[:])))

    return p
```

`laplace2d` takes five arguments, the `p` matrix, the `y`-vector, `dx`, `dy` and the value `l1norm_target`.  This last value defines how close the `p` matrix should be in two consecutive iterations before the loop breaks and returns the calculated `p` value.

Note that when executing the cells above in your own notebook, there will be no output.  You have *defined* the function but you have not yet *called* the function.  It is now available for you to use, the same as `numpy.linspace` or any other function in our namespace.

```python
##variable declarations
nx = 31
ny = 31
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)

##initial conditions
p = numpy.zeros((ny, nx))  # create a XxY vector of 0's

##plotting aids
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 1, ny)

##boundary conditions
p[:, 0] = 0  # p = 0 @ x = 0
p[:, -1] = y  # p = y @ x = 2
p[0, :] = p[1, :]  # dp/dy = 0 @ y = 0
p[-1, :] = p[-2, :]  # dp/dy = 0 @ y = 1
```

Now let's try using our `plot2D` function to look at our initial conditions.  If the function has been correctly defined, you should be able to begin typing `plot2D` and hit the **tab** key for auto-complete options.

```python
plot2D(x, y, p)
```

![png](output_19_0.png)

It worked!  This is the initial state of our problem, where the value of `p` is zero everywhere except for along $x=2$ where $p=y$.  Now let's try to run our `laplace2d` function with a specified L1 target of .01

[Hint: if you are having trouble remembering the order in which variables are sent to a function, you can just type `laplace2d(` and the iPython Notebook will put up a little popup box to remind you]

```python
p = laplace2d(p, y, dx, dy, 1e-4)
```

Now try plotting this new value of `p` with our plot function.

```python
plot2D(x, y, p)
```

![png](output_23_0.png)

***

## Learn More

The [next step](./13_Step_10.ipynb) will be to solve Poisson's equation. Watch **Video Lesson 11** on You Tube to understand why we need Poisson's equation in CFD.

```python
from IPython.display import YouTubeVideo
YouTubeVideo('ZjfxA3qq2Lg')
```

And for a detailed walk-through of the discretization of Laplace and Poisson equations (steps 9 and 10), watch **Video Lesson 12** on You Tube:

```python
from IPython.display import YouTubeVideo
YouTubeVideo('iwL8ashXhWU')
```

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

> (The cell above executes the style for this notebook.)

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
# Julia 1.0 compat
if VERSION < v"1.1"
    @eval isnothing(::Any) = false
    @eval isnothing(::Nothing) = true
end

"""
    has_card(deck, card)

Return true if `deck` contains `card`.
"""
has_card(deck, card) = card ∈ deck

"""
    find_card(deck, card)

Return the index of `card` in `deck`.
"""
find_card(deck, card) = findfirst(==(card), deck)

"""
    all_cards_even(deck)

Return true if all cards in `deck` are even.
"""
all_cards_even(deck) = all(iseven, deck)

"""
    any_odd_cards(deck)

Return true if any cards in `deck` are odd.
"""
any_odd_cards(deck) = any(isodd, deck)

"""
    first_even_card_idx(deck)

Return the index of the first even card in `deck`.
"""
first_even_card_idx(deck) = findfirst(iseven, deck)

"""
    first_odd_card(deck)

Return the first odd card in `deck`, and `nothing` if the deck does not contain any odd cards.
"""
function first_odd_card(deck)
    idx = findfirst(isodd, deck)
    isnothing(idx) ? nothing : deck[idx]
end

ext: jl


text:
# function to perform tests of the rotation functions in the Rotations module

#=
########################
# Define helper methods
########################

import Rotations: numel
numel{T <: Rotations.RotationTypes}(X::T) = Rotations.numel(Rotations.strip_eltype(T))
numel{T}(X::T) = (T <: FixedSizeArrays.FixedArray) || (T <: FixedSizeArrays.AbstractArray)  ? length(X) : error("numel undefined for input of type $(T)")

# a macro to test if the contents of two containers are approximatley equal
macro contents_approx_eq(a, b)
    quote
        n = numel($(esc(a)))
        @test n == numel($(esc(b)))
        for i = 1:n
            ai, bi = getindex($(esc(a)), i), getindex($(esc(b)), i)
            @test typeof(ai) == typeof(bi)
            @test_approx_eq ai bi
        end
    end
end
macro contents_approx_eq_eps(a, b, eps)
    quote
        n = numel($(esc(a)))
        @test n == numel($(esc(b)))
        for i = 1:n
            ai, bi = getindex($(esc(a)), i), getindex($(esc(b)), i)
            @test typeof(ai) == typeof(bi)
            @test_approx_eq_eps ai bi $(esc(eps))
        end
    end
end

# a macro to test if two types are approximatey equal
macro types_approx_eq(a, b)
    quote
        @test typeof($(esc(a))) == typeof($(esc(b)))
        @contents_approx_eq($(esc(a)), $(esc(b)))
    end
end
macro types_approx_eq_eps(a, b, eps)
    quote
        @test typeof($(esc(a))) == typeof($(esc(b)))
        @contents_approx_eq_eps($(esc(a)), $(esc(b)), $(esc(eps)))
    end
end

# a macro to test if the contents of two containers are approximatley equal, without examining the element types
macro contents_approx_eq_notype(a, b)
    quote
        n = numel($(esc(a)))
        @test n == numel($(esc(b)))
        for i = 1:n
            ai, bi = getindex($(esc(a)), i), getindex($(esc(b)), i)
            @test_approx_eq ai bi
        end
    end
end
=#

#####################################################################################
# build a full list of rotation types including the different angle ordering schemas
#####################################################################################

rot_types = (RotMatrix{3}, Quat, SPQuat, AngleAxis, RodriguesVec,
			 UnitQuaternion, RodriguesParam, MRP,
             RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
             RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ)

one_types = (RotX, RotY, RotZ)
two_types = (RotXY, RotYZ, RotZX, RotXZ, RotYX, RotZY)
taitbyran_types = (RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX)
all_types = (RotMatrix{3}, Quat, SPQuat, AngleAxis, RodriguesVec,
			 UnitQuaternion, RodriguesParam, MRP,
             RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
             RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ,
             RotX, RotY, RotZ,
             RotXY, RotYZ, RotZX, RotXZ, RotYX, RotZY)

quat_types = (Quat, UnitQuaternion)

###############################
# Start testing
###############################

@testset "Rotations Tests" begin

    ###############################
    # Check fixed relationships
    ###############################

    @testset "Identity rotation checks" begin
        I = one(SMatrix{3,3,Float64})
        I32 = one(SMatrix{3,3,Float32})
        @testset "$(R)" for R in all_types
            @test one(R)::R == I
            @test one(R{Float32})::R{Float32} == I32
        end
    end

    ################################
    # check on the inverse function
    ################################

    @testset "Testing inverse()" begin
        repeats = 100
        I = one(RotMatrix{3,Float64})
        @testset "$(R)" for R in all_types
            Random.seed!(0)
            for i = 1:repeats
                r = rand(R)
                @test inv(r) == adjoint(r)
                @test inv(r) == transpose(r)
                @test inv(r)*r ≈ I
                @test r*inv(r) ≈ I
            end
        end
    end

    #########################################################################
    # Rotate some stuff
    #########################################################################

    # a random rotation of a random point
    @testset "Rotate Points" begin
        repeats = 100
        @testset "$(R)" for R in all_types
            I = one(R)::R

            Random.seed!(0)
            for i = 1:repeats
                r = rand(R)
                m = SMatrix(r)
                v1 = randn(SVector{3})
                v2 = randn(3)

                @test I*v1 ≈ v1
                @test I*v2 ≈ v2
                @test r*v1 ≈ m*v1
                @test r*v2 ≈ m*v2
            end
        end
    end

    @testset "Quaternion double cover" begin
        repeats = 100
		@testset "$Q double cover" for Q in quat_types
	        for i = 1 : repeats
	            q = rand(Quat)

	            q2 = Quat(-q.w, -q.x, -q.y, -q.z) # normalize: need a tolerance
	            @test SVector(q2.w, q2.x, q2.y, q2.z) ≈ SVector(-q.w, -q.x, -q.y, -q.z) atol = 100 * eps()
	            @test q ≈ q2 atol = 100 * eps()

	            q3 = Quat(-q.w, -q.x, -q.y, -q.z, false) # don't normalize: everything is exact
	            @test (q3.w, q3.x, q3.y, q3.z) == (-q.w, -q.x, -q.y, -q.z)
	            @test q == q3

	            Δq = q \ q3
	            @test Δq ≈ one(Quat) atol = 100 * eps()
	        end
		end
    end

    # compose two random rotations
    @testset "Compose rotations" begin
        repeats = 100
        @testset "$(R1) * $(R2)" for R1 in all_types, R2 in all_types
            Random.seed!(0)
            for i = 1:repeats
                r1 = rand(R1)
                m1 = SMatrix(r1)

                r2 = rand(R2)
                m2 = SMatrix(r2)

                @test r1*r2 ≈ m1*m2
            end
        end
    end

    #########################################################################
    # Test conversions between rotation types
    #########################################################################
    @testset "Convert rotations" begin
        repeats = 100
        @testset "convert $(R1) -> $(R2)" for R1 in all_types, R2 in rot_types
            Random.seed!(0)
            for i = 1:repeats
                r1 = rand(R1)
                m1 = SMatrix(r1)

                r2 = R2(r1)

                @test r2 ≈ m1
            end
        end
    end

    #########################################################################
    # Test robustness of DCM to Quat function
    #########################################################################
    function nearest_orthonormal(not_orthonormal)
        u,s,v = svd(not_orthonormal)
        return u * Diagonal([1, 1, sign(det(u * transpose(v)))]) * transpose(v)
    end

    @testset "DCM to Quat" begin
        pert = 1e-3
		for quat in quat_types
	        for type_rot in all_types
	            for _ = 1:100
	                not_orthonormal = rand(type_rot) + rand(3, 3) * pert
	                @test norm(quat(not_orthonormal) - nearest_orthonormal(not_orthonormal)) < 10pert
	            end
	        end
		end
    end

    #########################################################################
    # Check angle and axis and inv work as expected
    #########################################################################

    @testset "Testing angle / axis extraction" begin
        repeats = 100
        @testset "$(R)" for R in rot_types
            Random.seed!(0)
            for i = 1:repeats
                r1 = rand(AngleAxis)

                angle = rotation_angle(r1)
                axis = rotation_axis(r1)

                r2 = R(r1)

                @test rotation_angle(r2) ≈ angle
                @test rotation_axis(r2) ≈ axis
            end
        end
		@test norm(rotation_axis(Quat(1.0, 0.0, 0.0, 0.0))) ≈ 1.0

        # TODO RotX, RotXY?
    end

    #########################################################################
    # Check construction of Quat given two vectors
    #########################################################################

    @testset "Testing construction of Quat given two vectors" begin
        angle_axis_test(from, to, rot, atol) = isapprox(rot * from * norm(to) / norm(from), to; atol = atol)

        for i = 1 : 10000
            from = randn(SVector{3, Float64})
            to = rand(SVector{3, Float64})
            rot = rotation_between(from, to)
            @test angle_axis_test(from, to, rot, 1e-10)
        end

        # degenerate cases
        for i = 1 : 10000
            from = randn(SVector{3, Float64})
            to = randn() * from # either from and to are aligned, or in opposite directions
            rot = rotation_between(from, to)
            # @show rot
            @test angle_axis_test(from, to, rot, 1e-7)
        end
        for direction = 1 : 3
            for i = 1 : 10000
                from = @SVector [ifelse(i == direction, 1., 0.) for i = 1 : 3] # unit vector in direction 'direction'
                to = randn() * from
                rot = rotation_between(from, to)
                @test angle_axis_test(from, to, rot, 1e-7)
            end
        end
        @test_throws ArgumentError rotation_between(zero(SVector{3}), rand(SVector{3}))
        @test_throws ArgumentError rotation_between(rand(SVector{3}), zero(SVector{3}))
    end

    #########################################################################
    # Check roll, pitch, yaw constructors
    #########################################################################

    @testset "Testing roll / pitch / yaw constructors" begin
        repeats = 100
        s = 1e-4
        @testset "$(R)" for R in taitbyran_types
            Random.seed!(0)
            for i = 1:repeats
                roll = s*(rand()-0.5)
                pitch = s*(rand()-0.5)
                yaw = s*(rand()-0.5)

                # This tests whether the rotations are the same to first order
                # in roll, pitch and yaw. Second-order terms are small enough
                # to pass the isapprox() test, but not first-order terms.
                r1 = RotXYZ(roll, pitch, yaw)
                r2 = R(roll=roll, pitch=pitch, yaw=yaw)

                @test r1 ≈ r2
            end
        end
    end

    #########################################################################
    # Check that isrotation works
    #########################################################################
    @testset "Testing isrotation" begin
      a=[1.0 0.0 0.0
         0.0 0.0 -1.0
         0.0 1.0 0.0]
      @test isrotation(a)
      foreach(rot_types) do rot_type
        foreach(1:20) do idx
          @test isrotation(rand(rot_type))
        end
      end
      a=[40.0 0.0 0.0
         0.0 0.0 1.0
         0.0 1.0 0.0]
      @test !isrotation(a)
      a=[1.0 0.0 0.0
         0.0 0.0 1.0
         0.0 1.0 0.0]
      @test !isrotation(a)

      # isrotation should work for integer (or boolean) matrices (issue #94)
      @test isrotation([0 1 0; -1 0 0; 0 0 1])
      @test isrotation(Matrix(I,3,3))
    end

    @testset "Testing type aliases" begin
        @test one(RotMatrix{2, Float64}) isa RotMatrix2{Float64}
        @test one(RotMatrix{3, Float64}) isa RotMatrix3{Float64}
    end

    @testset "Testing normalization" begin
        θ, x, y, z = 1., 2., 3., 4.
        aa = AngleAxis(θ, x, y, z)
        @test norm([aa.axis_x, aa.axis_y, aa.axis_z]) ≈ 1.
        aa = AngleAxis(θ, x, y, z, false)
        @test norm([aa.axis_x, aa.axis_y, aa.axis_z]) ≈ norm([x, y, z])

        w, x, y, z = 1., 2., 3., 4.
        quat = Quat(w, x, y, z)
        @test norm([quat.w, quat.x, quat.y, quat.z]) ≈ 1.
        quat = Quat(w, x, y, z, false)
        @test norm([quat.w, quat.x, quat.y, quat.z]) ≈ norm([w, x, y, z])

        w, x, y, z = 1., 2., 3., 4.
        quat = UnitQuaternion(w, x, y, z)
        @test norm([quat.w, quat.x, quat.y, quat.z]) ≈ 1.
        quat = UnitQuaternion(w, x, y, z, false)
        @test norm([quat.w, quat.x, quat.y, quat.z]) ≈ norm([w, x, y, z])
    end

    @testset "Testing RotMatrix conversion to Tuple" begin
        rot = one(RotMatrix{3, Float64})
        @inferred Tuple(rot)
    end

    @testset "Testing show" begin
        io = IOBuffer()
        r = rand(RotMatrix{2})
        show(io, MIME("text/plain"), r)
        str = String(take!(io))
        @test startswith(str, "2×2 RotMatrix{2,Float64,4}")

        rxyz = RotXYZ(1.0, 2.0, 3.0)
        show(io, MIME("text/plain"), rxyz)
        str = String(take!(io))
        @test startswith(str, "3×3 RotXYZ{Float64}") && occursin("(1.0, 2.0, 3.0):", str)
    end
end


text:
using Test
using Documenter

@test isdefined(Documenter.Writers, :enable_backend)
@test isdefined(Documenter.Writers, :backends_enabled)

@test Documenter.Writers.backends_enabled[:latex] === false

using DocumenterLaTeX

@test Documenter.Writers.backends_enabled[:latex] === true


text:
#!/usr/bin/env julia

##########################################################################
# run_convergence-randomlargescale.jl
#
# Convergence as a function of number of iterations.
##########################################################################

using CoordinatedPrecoding
using Compat, JLD

##########################################################################
# Set up logging
Lumberjack.configure(Lumberjack._lumber_mill.timber_trucks["console"]; mode = "warn")
Lumberjack.add_truck(Lumberjack.LumberjackTruck("debug.log", "debug"), "debug")

##########################################################################
# General settings
srand(83196723)
start_time = Libc.strftime("%Y%m%dT%H%M%S", time())

##########################################################################
# RandomLargeScale network
simulation_params = @compat Dict(
    "simulation_name" => "convergence_$(start_time)-randomlargescale",
    "I" => 10, "Kc" => 1, "N" => 2, "M" => 2, "d" => 1,
    "Ndrops" => 10, "Nsim" => 10,
    "precoding_methods" => [
        Shi2011_WMMSE,
        Gomadam2008_MaxSINR,
        Komulainen2013_WMMSE,
        # Razaviyayn2013_MinMaxWMMSE,
        Eigenprecoding
    ],
    "aux_precoding_params" => Dict(
        "initial_precoders" => "eigendirection",
        "stop_crit" => 0.,
        "max_iters" => 20,
    ),
    "aux_independent_variables" => [
        (set_transmit_powers_dBm!, [10, 30]),
    ],
)
network =
    setup_random_large_scale_network(simulation_params["I"],
        simulation_params["Kc"], simulation_params["N"], simulation_params["M"],
        num_streams=simulation_params["d"])
raw_results = simulate_precoding_convergence(network, simulation_params)

println("-- Saving $(simulation_params["simulation_name"]) results")
save("$(simulation_params["simulation_name"]).jld",
     "simulation_params", clean_simulation_params_for_jld(simulation_params),
     "raw_results", raw_results)


text:
# https://adventofcode.com/2021/day/7
using AdventOfCode

input = parse.(readlines("2021/data/day_7.txt"))[begin]

function part_1(input)
    med=Int(median(input))
    sum(abs∘Base.Fix2(-,med),input)
end
function part_2(input)
    μ=mean(input)
    minimum(sum((x->x^2+abs(x))∘Base.Fix2(-,m),input)÷2 for m in (floor(Int,μ),ceil(Int,µ)))
end


text:
import Pkg; Pkg.activate(".")
include("../required.jl")

## Test landcover variables
# Define coordinates range
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)

# Test loading variables
lc_vars = SimpleSDMPredictor(Copernicus, LandCover, 1:10; resolution = 10.0, coords...)

# Load worldclim variables to compare
wc_vars = SimpleSDMPredictor(WorldClim, BioClim, 1:19; resolution = 10.0, coords...)

## Load with geotiff instead

# Check which layer is crops
glossary = CSV.read(joinpath("data", "proc", "glossary.csv"), DataFrame) # lc2

# Load crops with geotiff
lc2_geo_full = geotiff(SimpleSDMPredictor, joinpath("assets", "landcover", "lc_crops_10m.tif"))
lc2_geo = geotiff(SimpleSDMPredictor, joinpath("assets", "landcover", "lc_crops_10m.tif"); coords...)
lc2_geo = convert(Float32, lc2_geo)
# Get corresponding layer
lc2_csv = lc_vars[2]

# Tests
isequal(lc2_geo, lc2_csv)
isequal(lc2_geo.grid, lc2_csv.grid)
isapprox(lc2_geo.grid, lc2_csv.grid)
isequal(collect(lc2_geo), collect(lc2_csv))
isequal(length(lc2_geo), length(lc2_csv)) # not same number of elements...

testdf = DataFrame([lc2_geo, lc2_csv]) # not same top coordinates

lc2_geo
lc2_csv
isapprox(stride(lc2_geo, 1), stride(lc2_csv, 1)) # not same stride...

size(lc2_geo_full) # 840 x 2159, but should be 840 x 2160...

## Check details from landcover function

# Arguments
layers = collect(1:10)
resolution = 10.0
path = joinpath("assets", "landcover/")

# Function content
# function landcover(layers::Vector{Int64}; resolution::AbstractFloat=10.0, path::AbstractString=joinpath("assets", "landcover/"))
    ## Get file paths
    # List files in path
    lc_files = readdir("$(path)")
    # Filter for selected resolution
    filter!(x -> occursin.("$(Int(resolution))m.tif", x), lc_files)
    # Create paths for selected layers only
    paths = joinpath.(path, lc_files)[layers]

    ## Load data
    # Read raster dataset
    datasets = [ArchGDAL.read(p) for p in paths]
    datasets[2] # 2160 x 240
    # Read raster values from band 1 (only band in this case)
    values = [ArchGDAL.read(d, 1) for d in datasets]
    values[2] # 2160 x 840
    # HOW???
    # Convert from UInt8 to Float64
    values_int = [convert.(Float64, v) for v in values]

    ## Preliminary manipulations
    # Reverse rows and permute dimensions
    landcover_mat = [permutedims(v[:,end:-1:1]) for v in values_int]
    # Replace 255 (default no data values) by nothing
    landcover_mat = convert.(Array{Union{Float32, Nothing}}, landcover_mat)
    [replace!(l, 255 => nothing) for l in landcover_mat]

    # Fill missing latitudes with nothing (latitude extent is only (-60,80) for landcover data)
    nlat, nlon = size(landcover_mat[1])
    slim, nlim = abs(-60), 80
    res = Int64(nlat/(nlim+slim))
    south_nothings = fill(nothing, ((90-slim)*res, nlon))
    north_nothings = fill(nothing, ((90-nlim)*res, nlon))
    landcover_grids = [vcat(south_nothings, l, north_nothings) for l in landcover_mat]

    # Convert to SimpleSDMLayers
    landcover_layers = SimpleSDMPredictor.(landcover_grids, -180.0, 180.0, -90.0, 90.0)

    # return landcover_layers
# end

[geotiff(SimpleSDMPredictor, p) for p in paths] # always returning 840 x 2159...

coords_landcover = (left = -180.0, right = 180.0, bottom = -60.0, top = 80.0)
geotiff(SimpleSDMPredictor, paths[2]; coords_landcover...) # still the same

## Check details from geotiff
file = paths[2]
bandnumber = 1
left = -180.0
right = 180.0
bottom = -90.0
top = 90.0

# function geotiff(
#     ::Type{LT},
#     file::AbstractString,
#     bandnumber::Integer=1;
#     left = -180.0,
#     right = 180.0,
#     bottom = -90.0,
#     top = 90.0
# ) where {LT<:SimpleSDMLayer}

    # This next block is reading the geotiff file, but also making sure that we
    # clip the file correctly to avoid reading more than we need.
    # This next block is reading the geotiff file, but also making sure that we
    # clip the file correctly to avoid reading more than we need.
    # layer = ArchGDAL.read(file) do dataset
    dataset = ArchGDAL.read(file)

        transform = ArchGDAL.getgeotransform(dataset)
        wkt = ArchGDAL.getproj(dataset)

        # The data we need is pretty much always going to be stored in the first
        # band, but this is not the case for the future WorldClim data.
        band = ArchGDAL.getband(dataset, bandnumber)
        T = ArchGDAL.pixeltype(band)

        # The nodata is not always correclty identified, so if it is not found, we assumed it is the smallest value in the band
        nodata = isnothing(ArchGDAL.getnodatavalue(band)) ? convert(T, ArchGDAL.minimum(band)) : convert(T, ArchGDAL.getnodatavalue(band))

        # Get the correct latitudes
        minlon = transform[1]
        maxlat = transform[4]
        maxlon = minlon + size(band,1)*transform[2]
        minlat = maxlat - abs(size(band,2)*transform[6])

        left = isnothing(left) ? minlon : max(left, minlon)
        right = isnothing(right) ? maxlon : min(right, maxlon)
        bottom = isnothing(bottom) ? minlat : max(bottom, minlat)
        top = isnothing(top) ? maxlat : min(top, maxlat)

        lon_stride, lat_stride = transform[2], transform[6]

        width = ArchGDAL.width(dataset)
        height = ArchGDAL.height(dataset)

        #global lon_stride, lat_stride
        #global left_pos, right_pos
        #global bottom_pos, top_pos

        lon_stride, left_pos, min_width = SimpleSDMLayers._find_span(width, minlon, maxlon, left)
        _, right_pos, max_width = SimpleSDMLayers._find_span(width, minlon, maxlon, right) # here's the problem!!
        lat_stride, top_pos, max_height = SimpleSDMLayers._find_span(height, minlat, maxlat, top)
        _, bottom_pos, min_height = SimpleSDMLayers._find_span(height, minlat, maxlat, bottom)

        max_height, min_height = height .- (min_height, max_height) .+ 1

        # We are now ready to initialize a matrix of the correct type.
        buffer = Matrix{T}(undef, length(min_width:max_width), length(min_height:max_height))
        ArchGDAL.read!(dataset, buffer, bandnumber, min_height:max_height, min_width:max_width)
        buffer = convert(Matrix{Union{Nothing,eltype(buffer)}}, rotl90(buffer))
        replace!(buffer, nodata => nothing)
        # LT(buffer, left_pos-0.5lon_stride, right_pos+0.5lon_stride, bottom_pos-0.5lat_stride, top_pos+0.5lat_stride)
        SimpleSDMPredictor(buffer, left_pos-0.5lon_stride, right_pos+0.5lon_stride, bottom_pos-0.5lat_stride, top_pos+0.5lat_stride)
    end

    # return layer

# end

geotiff(SimpleSDMPredictor, paths[2]; left=nothing, right=nothing, bottom=nothing, top=nothing)
lc2_geo_full

# Extract problem
_, right_pos, max_width = SimpleSDMLayers._find_span(width, minlon, maxlon, right) # here's the problem!!

# Check details of _find_span
n = width
m = minlon
M = maxlon
pos = right
# function _find_span(n, m, M, pos)
    pos > M && return nothing
    pos < m && return nothing
    stride2 = (M - m) / n
    centers = (m + 0.5stride2):stride2:(M-0.5stride2)
    span_pos = last(findmin(abs.(pos .- centers)))
    return (stride, centers[span_pos], span_pos)
# end

lc2_geo_fix = similar(lc2_csv)
lc2_geo_fix.grid = copy(lc2_geo.grid)
lc2_geo_fix = convert(SimpleSDMPredictor, lc2_geo_fix)

compdf = DataFrame([lc2_geo_fix, lc2_csv])
rename!(compdf, :x1 => :geo_fix, :x2 => :csv)
allowmissing!(compdf)
for col in [:geo_fix, :csv]
    replace!(compdf[!, col], nothing => missing)
end
compdf
dropmissing(compdf, :geo_fix)

plot(lc2_geo_fix, c = :viridis)
plot(lc2_csv, c = :viridis) # one cell misplacement

## Attempt to load full 100m resolution layer
megafilepath = "./assets/landcover/test/landcover_copernicus_global_100m_v2.0.2_moss.tif"
megafilepath = "./assets/landcover/landcover_copernicus_global_100m_v2.0.2_urban.tif"
isfile(megafilepath)
test = geotiff(SimpleSDMPredictor, megafilepath; coords...)
ArchGDAL.read(megafilepath)
dataset = ArchGDAL.read(megafilepath)
ArchGDAL.getgeotransform(dataset)
ArchGDAL.read("./assets/landcover/lc_bare_10m.tif")
ArchGDAL.getgeotransform(ArchGDAL.read("./assets/landcover/lc_bare_10m.tif"))

wc1_path = "/home/gdansereau/.data/SimpleSDMLayers.jl/assets/WorldClim/BioClim/10/wc2.1_10m_bio_1.tif"
geotiff(SimpleSDMPredictor, wc1_path)
wc1_dataset = ArchGDAL.read(wc1_path)
ArchGDAL.getgeotransform(wc1_dataset)

## In bash
# -tap option
gdalwarp -tr 0.166667 0.166667 -tap -r average --config GDAL_CACHEMAX 500 -wm 500 -multi assets/landcover/lc_moss_10m.tif -overwrite assets/landcover/landcover_test.tif

# -te option
gdalwarp -tr 0.166667 0.166667 -te -180.0 -60.0 180.0 80.0 -r average --config GDAL_CACHEMAX 500 -wm 500 -multi assets/landcover/lc_moss_10m.tif -overwrite assets/landcover/landcover_test.tif

# -te option and same stride as WorldClim
gdalwarp -tr 0.1666666666666666574 0.1666666666666666574 -te -180.0 -60.0 180.0 80.0 -r average --config GDAL_CACHEMAX 500 -wm 500 -multi assets/landcover/lc_moss_10m.tif -overwrite assets/landcover/landcover_test.tif

# -te option and same stride as WorldClim on full layer
gdalwarp -tr 0.1666666666666666574 0.1666666666666666574 -te -180.0 -60.0 180.0 80.0 -r average --config GDAL_CACHEMAX 500 -wm 500 -multi assets/landcover/landcover_copernicus_global_100m_v2.0.2_moss.tif -overwrite assets/landcover/landcover_fulltest.tif

# full spatial extent
gdalwarp -tr 0.1666666666666666574 0.1666666666666666574 -te -180.0 -90.0 180.0 90.0 -r average --config GDAL_CACHEMAX 500 -wm 500 -multi assets/landcover/landcover_copernicus_global_100m_v2.0.2_moss.tif -overwrite assets/landcover/landcover_fulltest.tif

# Back in Julia
test_path = "./assets/landcover/landcover_test.tif"
test_dataset = ArchGDAL.read(test_path)
ArchGDAL.getgeotransform(test_dataset)

moss_path = "./assets/landcover/lc_moss_10m.tif"
moss_dataset = ArchGDAL.read(moss_path)
ArchGDAL.getgeotransform(moss_dataset)

fulltest_path = "./assets/landcover/landcover_fulltest.tif"
fulltest_dataset = ArchGDAL.read(fulltest_path)
ArchGDAL.getgeotransform(fulltest_dataset)

geotiff(SimpleSDMPredictor, test_path)
geotiff(SimpleSDMPredictor, fulltest_path)
test_layer = convert(Float32, geotiff(SimpleSDMPredictor, test_path; coords...))
fulltest_layer = convert(Float32, geotiff(SimpleSDMPredictor, fulltest_path; coords...))
wctest_layer = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
moss_indx = readdir("assets/landcover/") |> x ->
    filter(f -> occursin.("10m.tif", f), x) |> x ->
    findfirst(occursin.("moss", x))
lc_moss = lc_vars[moss_indx]

isequal(latitudes(fulltest_layer), latitudes(wctest_layer))
isequal(longitudes(fulltest_layer), longitudes(wctest_layer))
SimpleSDMLayers._layers_are_compatible(fulltest_layer, wctest_layer)

plot(test_layer)
plot(lc_moss)
plot(fulltest_layer)

isequal(lc_moss, fulltest_layer)
isequal(lc_moss.grid, fulltest_layer.grid)
testdf = DataFrame([lc_moss, fulltest_layer]) # not exactly same coordinates

isapprox(lc_moss.bottom, fulltest_layer.bottom)

moss_new = similar(lc_moss)
moss_new.grid = copy(fulltest_layer.grid)
moss_new = convert(SimpleSDMPredictor, moss_new)

testdf = DataFrame([lc_moss, moss_new])
isequal(testdf.x1, testdf.x2)
isapprox(testdf.x1, testdf.x2)
isequal(collect(lc_moss), collect(moss_new))
isapprox(collect(lc_moss), collect(moss_new))

writedlm("test.csv", collect(lc_moss))
writedlm("test.csv", collect(moss_new))

allowmissing!(testdf)
for col in [:x1, :x2]
    replace!(testdf[!, col], nothing => missing)
end
testdf
filter(x -> ismissing(x.x1) && !ismissing(x.x2), testdf)
filter(x -> !ismissing(x.x1) && ismissing(x.x2), testdf)
isequal(nrow(dropmissing(testdf)), length(lc_moss))

dropmissing!(testdf)
testdiff = filter(x -> x.x1 != x.x2, testdf) # 800 sites with difference
insertcols!(testdiff, :diff => testdiff.x2 .- testdiff.x1)
unique(testdiff.diff)
sort(countmap(testdiff.diff))

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
#This file contains functions that are pervasive throughout this package as well as structures to define common objects like cells of input variables

####################################################
# Helper Functions
####################################################

⊗(A,B) = kron(A,B) #just makes things look nice

#Kronecker delta function
#δ(x::T, y::T) where {T<:Number} = isequal(x,y) ? one(T) : zero(T) #this typing is excessive, but good practice
δ(x, y) = isequal(x,y) ? 1 : 0

#The following functions help generate adjacency matrices for the underlying network. Circulant arrays are used for periodic boundaries and off-diagonal arrays are used for non-periodic graphs
#See https://stackoverflow.com/a/45958661

#Creates a circulant array. Rows are shifted over one down the array
function circulant!(A,n)

    vec = spzeros(Int, n)
    vec[[2,n]] .= 1 #put values into the 2nd and last entry

    for i in 1:n
        A[i,:] = vec
        circshift!(vec, vec, 1) #shift the row over by 1
    end

    return A
end

#Fills the diagonals above and below main diagonal
function offDiags!(A)
    A[diagind(A,1)] .= 1
    A[diagind(A,-1)] .= 1
    return A
end

#Generates a square (n,n) sparse adjacency matrix to convert into a network
function genAdj(n::Int, isPeriodic::Bool)
    A = spzeros(Int, n,n)

    if isPeriodic #Does the graph wrap around?
        circulant!(A,n)
    else
        offDiags!(A)
    end

    return A ⊗ I(n) + I(n) ⊗ A
end

#Generates a rectangular (m,n) sparse adjacency matrix to convert into a network
function genAdj(m::Int, n::Int, isPeriodic::Bool)
    Am = spzeros(Int, m, m)
    An = spzeros(Int, n, n)

    if isPeriodic #Does the graph wrap around?
        circulant!(Am,m)
        circulant!(An,n)
    else
        offDiags!(Am)
        offDiags!(An)
    end

    return An ⊗ I(m) + I(n) ⊗ Am
end

#Generates a 3D (l,m,n) sparse adjacency matrix to convert into a network
function genAdj(n::Int, m::Int, l::Int, isPeriodic::Bool)
    Al = spzeros(Int, l,l)
    Am = spzeros(Int, m,m)
    An = spzeros(Int, n,n)

    if isPeriodic #Does the graph wrap around?
        circulant!(Al,l)
        circulant!(Am,m)
        circulant!(An,n)
    else
        offDiags!(Al)
        offDiags!(Am)
        offDiags!(An)
    end

    return I(l) ⊗ (Am ⊗ I(n) + I(m) ⊗ An) + Al ⊗ I(m*n) #This took a while to figure out...
end

####################################################
# Hamiltonian Penalties
####################################################

"Define abstract type for different penalties"
abstract type Hamiltonian end

mutable struct AdhesionPenalty <: Hamiltonian
    J::OffsetMatrix{Int, Matrix{Int}} #J[n,m] gives the adhesion penality for cells with types n and m

    function AdhesionPenalty(J::Matrix{Int})
        issymmetric(J) ? nothing : error("J needs to be symmetric")

        Joff = OffsetArray(J, CartesianIndex(0, 0):CartesianIndex(size(J).-1))
        return new(Joff)
    end
end

mutable struct VolumePenalty <: Hamiltonian
    desiredVolumes::OffsetVector{Int,Vector{Int64}}
    λᵥ::OffsetVector{Int,Vector{Int64}}

    function VolumePenalty(desiredVolumes::Vector{Int}, λᵥ::Vector{Int})
        desiredVolumesOff = OffsetVector([0; desiredVolumes], 0:length(desiredVolumes))
        λᵥOff = OffsetVector([0; λᵥ], 0:length(λᵥ))
        return new(desiredVolumesOff, λᵥOff)
    end
end

####################################################
# Structure for user to provide information
####################################################

Base.@kwdef struct ModelParameters{N} #N is graph dimension (makes types stable)
    graphDimension::NTuple{N, Int} = (100,100) #What are the dimensions of the space?
    isPeriodic::Bool = true #Do you want the space to wrap around?
    cellTypes::Vector{String} = ["TCell"] #What kinds of cells do you want?
    cellCounts::Vector{Int} = [10] #How many of each cell do you want?
    cellVolumes::Vector{Int} = [200] #Desired size for each cell? (need same type cells with different sizes?)
    penalties::Vector{Hamiltonian} = [AdhesionPenalty([0 50; 50 100]), VolumePenalty(fill(200,10),[5])] #What penalties are added to the model?
    temperature::Float64 = 20.0 #Temperature to weight site flips (higher increases chance of accepting energy increasing flips)
end

####################################################
# Structure for graph information
####################################################

mutable struct NamedGraph
    network::SimpleGraph{Int} #a network of nodes and edges equivalent to the grid
    σ::Vector{Int} #Cell ID for each node
    τ::Vector{String} #Each node has a cell type
    isArticulation::BitVector #Check for whether changing a node's ID will dissconnent the cell

    #Two inner contructors, the first assume you have all the individual fields
    NamedGraph(network::SimpleGraph{Int}, σ::Vector{Int}, τ::Vector{String}, isArticulation::BitVector) = new(network, σ, τ, isArticulation)

    function NamedGraph(gr::SimpleGraph{Int}, cellMembership::Array{Int}, M::ModelParameters{N}) where N

        #Cell ID for each node
        σ = vec(cellMembership)

        #Initialize a list of cell types proportional to the number of that type
        #Shuffle is to make use all the same times are not next to each other
        cellAssignments = shuffle(inverse_rle(M.cellTypes, M.cellCounts)) #inverse_rle(["a","b"], [2,3]) = ["a","a","b","b","b"]

        #Each node has a cell type
        τ = [σᵢ == 0 ? "Medium" : cellAssignments[σᵢ] for σᵢ in σ]

        #Additionally, each node can be a articulation point (removing will disconnect the cell)
        isArticulation = falses(length(σ))

        graph = new(gr, σ, τ, isArticulation)

        #Update the articulation points
        UpdateConnections!(graph; checkConnect=true)

        return graph
    end
end

####################################################
# Structure to hold cell level attributes
####################################################

mutable struct CellAttributes
    #Vectors are offset to include medium (medium gets index 0, cell 1 gets index 1, etc.)
    ids::OffsetVector{Int,Vector{Int}} #vector of cell IDs
    volumes::OffsetVector{Int,Vector{Int}} #vector of cell volumes
    desiredVolumes::OffsetVector{Int,Vector{Int}} #vector of desired cell volumes
    types::OffsetVector{String,Vector{String}} #given a cell index, output it's type
    typeMap::Dict{String, Int} #mapping cell types (e.g. "Medium") to an index (e.g. 0)

    function CellAttributes(graph::NamedGraph, M::ModelParameters)
        #How many cells are there?
        totalCells = sum(M.cellCounts)

        ids = OffsetVector(collect(0:totalCells), 0:totalCells)

        volumes = OffsetVector([sum(x->x==i, graph.σ) for i in 0:totalCells], 0:totalCells)

        types = OffsetVector(["Medium"; inverse_rle(M.cellTypes, M.cellCounts)], 0:totalCells)

        #Look for the VolumePenalty and use it to update the desired volumes
        volIdx = findfirst(p -> isa(p,VolumePenalty), M.penalties)
        desiredVolumes = M.penalties[volIdx].desiredVolumes #This is type unstable but it only happens once

        typeMap = Dict( M.cellTypes .=> 1:length(M.cellTypes) )
        typeMap["Medium"] = 0

        return new(ids, volumes, desiredVolumes, types, typeMap)
    end
end

####################################################
# Variables for Markov Step
####################################################

Base.@kwdef mutable struct MHStepInfo{T<:Int}
    sourceNode::T = 1 #Index of node choosen
    sourceNodeNeighbors::Vector{T} = [1] # Indicies for the neighboring nodes
    possibleCellTargets::Vector{T} = [1] # Unique cell IDs of the neighboring nodes
    sourceCell::T = 1 #ID of sourceNode
    targetCell::T = 1 #ID of chosen cell target
end

####################################################
# Model Structure
####################################################

mutable struct CellPotts{N}
    M::ModelParameters{N} #input parameters from the user
    cell::CellAttributes #summary of cell attributes
    graph::NamedGraph #node properties: cell membership, cell type, bridge status
    energy::Int #Total penality energy across graph
    visual::Array{Int,N} #An array of cell memberships for plotting
    stepCounter::Int #counts the number of MHSteps performed
    stepInfo::MHStepInfo #Information about nodes choosen for flip attempts

    function CellPotts(M::ModelParameters{N}) where N

        #Generate an adjacency matrix
        adjMat = genAdj(M.graphDimension..., M.isPeriodic)

        #create a graph based on the adjacency matrix
        gr = SimpleGraph(adjMat)

        #Choose a cell initialization method
        cellMembership = GrowCells(gr, M)

        #Add attributes to the nodes
        graph = NamedGraph(gr, cellMembership, M)

        #create a summary of the cells
        cell = CellAttributes(graph, M)

        #Create an instance of the model with zero energy
        CPM = new{N}(M, cell, graph, 0, cellMembership, 0, MHStepInfo())

        #Calculate the energy from the given penalties
        #Update the energy with the given penalties
        CPM.energy = sum([f(CPM) for f in M.penalties]) #😮

        return CPM
    end
end

####################################################
# Penalty Functions
####################################################

#=
Need to define two methods for each penality:
    - calculate penality on the entire grid
    - calculate change in penality after markov step
=#

#----Calculate on entire Model (for initialization)----

#redirect to method below
(AP::AdhesionPenalty)(CPM::CellPotts) = AP(CPM.graph, CPM.cell.typeMap)

#This code gets repeated a lot so I made it a separate function
function (AP::AdhesionPenalty)(graph::NamedGraph, typeMap::Dict{String, Int64})

    #Initialize energy
    H = 0

    #Loop through all edges in network and calculate adhesion adjacencies
    for edge in edges(graph.network)

        #Given a node index, get the cell id and type
        (σᵢ, σⱼ) = graph.σ[ [edge.src, edge.dst] ]
        (typeᵢ, typeⱼ) = graph.τ[ [edge.src, edge.dst] ]

        #Convert the type (string) to an index for J
        (τᵢ, τⱼ) = ( typeMap[typeᵢ], typeMap[typeⱼ] )

        H += AP.J[τᵢ, τⱼ] * (1-δ(σᵢ, σⱼ))
    end

    return H
end

function (VP::VolumePenalty)(CPM::CellPotts)
    #Create slots to add up penalties for each cell type
    ΔV = similar(VP.λᵥ)
    ΔV .= 0 #Making an Offset zeros vector

    #Loop through cells, see how far they are from a desired volume, and put into appropriate slot
    for id in CPM.cell.ids
        cellTypeIdx = CPM.cell.typeMap[ CPM.cell.types[id] ] # cell id (21) → cell type ("Medium") → cell type index (0)
        ΔV[cellTypeIdx] += VP.λᵥ[cellTypeIdx] * (CPM.cell.volumes[id] - CPM.cell.desiredVolumes[id])^2
    end

    #Calculate the penalty
    return sum(ΔV)
end

#----Calculate after markov step----

function (AP::AdhesionPenalty)(CPM::CellPotts, stepInfo::MHStepInfo)

    #Just want to look at the local neighborhood
    subGraphIdx = vcat(stepInfo.sourceNode, stepInfo.sourceNodeNeighbors)

    subgraph = CPM.graph[subGraphIdx] #uses the custom getindex to subset all fields

    sourceH = AP(subgraph, CPM.cell.typeMap)

    #Change the subgraph to calculate H for the target node
    subgraph.σ[1] = stepInfo.targetCell
    subgraph.τ[1] = CPM.cell.types[stepInfo.targetCell]

    #Loop through all edges in network and calculate adhesion adjacencies
    targetH = AP(subgraph, CPM.cell.typeMap)

    #Calculate the adhesion difference between target and source
    return targetH - sourceH
end

function (VP::VolumePenalty)(CPM::CellPotts, stepInfo::MHStepInfo)

    #Create a vector of source and target memberships
    sourceTarget = [stepInfo.sourceCell, stepInfo.targetCell]

    #Extract cell types
    cellTypeIdx = [CPM.cell.typeMap[type] for type in CPM.cell.types[sourceTarget]]

    #Extract the desired volumes
    desiredVolumes = CPM.cell.desiredVolumes[sourceTarget]

    #Extract the current volumes
    currentVolumes = CPM.cell.volumes[sourceTarget]

    #Calculate the target volumes
    targetVolumes = currentVolumes .+ [-1, 1]

    #Return the delta (Target - Source)
    return sum(@. VP.λᵥ[cellTypeIdx]*(targetVolumes - desiredVolumes)^2) - sum(@. VP.λᵥ[cellTypeIdx]*(currentVolumes - desiredVolumes)^2)
end

####################################################
# Helper Functions that need structs defined
####################################################

function UpdateConnections!(graph::NamedGraph; checkConnect::Bool=false)

    #Reset the articulation points (is this needed?)
    graph.isArticulation .= falses(size(graph.isArticulation))

    #Loop through cells to find articulation points
    for σᵢ in unique(graph.σ)

        #Get the subgraph for a given cell ID
        cellIdx = findall(isequal(σᵢ), graph.σ)
        subgraph = graph.network[cellIdx]

        if checkConnect
            if !is_connected(subgraph) #if not connected
                error("some cells are disconnected, try rerunning or use a different cell initialization")
            end
        end

        #Update the articulation points
        graph.isArticulation[cellIdx[articulation(subgraph)]] .= true
    end

    return nothing
end

#If not the source or target cell, then the articulations points will not change...right?
function UpdateConnections!(graph::NamedGraph, stepInfo::MHStepInfo)

    #Loop through cells that were changed
    for σᵢ in [stepInfo.sourceCell, stepInfo.targetCell]

        #Get the subgraph for a given cell ID
        cellIdx = findall(isequal(σᵢ), graph.σ)
        subgraph = graph.network[cellIdx] #<------------- THIS IS SLOW

        #Update the articulation points
        graph.isArticulation[cellIdx] .= false
        graph.isArticulation[cellIdx[articulation(subgraph)]] .= true
    end

    return nothing
end

#Overload getindex to get subset of Namedgraph
Base.getindex(graph::NamedGraph, i) = NamedGraph(graph.network[i], graph.σ[i], graph.τ[i], graph.isArticulation[i])

####################################################
# Override Base.show for each struct
####################################################

function Base.show(io::IO, CPM::CellPotts)
    println("Cell Potts Model:")
    #Grid
    dim = length(CPM.M.graphDimension)
    if dim == 1
        println("$(CPM.M.graphDimension[1])×$(CPM.M.graphDimension[1]) Grid")
    elseif dim == 2
        println("$(CPM.M.graphDimension[1])×$(CPM.M.graphDimension[2]) Grid")
    else
        println("$(CPM.M.graphDimension[1])×$(CPM.M.graphDimension[2])×$(CPM.M.graphDimension[3]) Grid")
    end

    #Cells and types
    print("Cell Counts:")
    for (key, value) in countmap(CPM.cell.types) #remove medium
        if key ≠ "Medium"
        print(" [$(key) → $(value)]")
        end
    end

    if length(CPM.M.cellTypes) > 1
        println(" [Total → $(length(CPM.cell.ids))]")
    else
        print("\n")
    end

    print("Model Penalties:")
    for p in typeof.(CPM.M.penalties)
        print(" $(p)")
    end
    print("\n")
    println("Current Energy: ", CPM.energy)
    println("Grid Temperature: ", CPM.M.temperature)
    println("Steps: ", CPM.stepCounter)
end


text:
#=
SandiaExample:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-12
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/Sandia.jl")
dim = 10000
matrix = GBMatrix(sparse(rand(0:1, dim, dim)))
@time sandia(matrix)


text:
@testset "MCIntegrationVI" begin
    C = 10.0
    i = MCIntegrationVI(; clipping=C)
    @test repr(i) == "Numerical Inference by Monte Carlo Integration"
    l = GaussianLikelihood()
end


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

## SubString, RevString and Cstring tests ##

## SubString tests ##
u8str = "∀ ε > 0, ∃ δ > 0: |x-y| < δ ⇒ |f(x)-f(y)| < ε"
u8str2 = u8str^2
len_u8str = length(u8str)
slen_u8str = length(u8str)
len_u8str2 = length(u8str2)
slen_u8str2 = length(u8str2)

@test len_u8str2 == 2 * len_u8str
@test slen_u8str2 == 2 * slen_u8str

u8str2plain = String(u8str2)

for i1 = 1:length(u8str2)
    if !isvalid(u8str2, i1); continue; end
    for i2 = i1:length(u8str2)
        if !isvalid(u8str2, i2); continue; end
        @test length(u8str2[i1:i2]) == length(u8str2plain[i1:i2])
        @test length(u8str2[i1:i2]) == length(u8str2plain[i1:i2])
        @test u8str2[i1:i2] == u8str2plain[i1:i2]
    end
end

str="tempus fugit"              #length(str)==12
ss=SubString(str,1,length(str)) #match source string
@test length(ss)==length(str)

ss=SubString(str,1,0)    #empty SubString
@test length(ss)==0

ss=SubString(str,14,20)  #start indexed beyond source string length
@test length(ss)==0

ss=SubString(str,10,16)  #end indexed beyond source string length
@test length(ss)==3

str2=""
ss=SubString(str2,1,4)  #empty source string
@test length(ss)==0

ss=SubString(str2,1,1)  #empty source string, identical start and end index
@test length(ss)==0

@test SubString("foobar",big(1),big(3)) == "foo"

str = "aa\u2200\u2222bb"
u = SubString(str, 3, 6)
@test length(u)==2
b = IOBuffer()
write(b, u)
@test String(take!(b)) == "\u2200\u2222"

@test_throws ArgumentError SubString(str, 4, 5)
@test_throws BoundsError next(u, 0)
@test_throws BoundsError next(u, 7)
@test_throws BoundsError getindex(u, 0)
@test_throws BoundsError getindex(u, 7)
@test_throws BoundsError getindex(u, 0:1)
@test_throws BoundsError getindex(u, 7:7)
@test reverseind(u, 1) == 4
@test typeof(Base.cconvert(Ptr{Int8},u)) == SubString{String}
@test Base.cconvert(Ptr{Int8},u) == u

str = "føøbar"
u = SubString(str, 4, 3)
@test length(u)==0
b = IOBuffer()
write(b, u)
@test String(take!(b)) == ""

str = "føøbar"
u = SubString(str, 10, 10)
@test length(u)==0
b = IOBuffer()
write(b, u)
@test String(take!(b)) == ""

# search and SubString (issue #5679)
str = "Hello, world!"
u = SubString(str, 1, 5)
@test rsearch(u, "World") == 0:-1
@test rsearch(u, 'z') == 0
@test rsearch(u, "ll") == 3:4

# sizeof
@test sizeof(SubString("abc\u2222def",4,4)) == 3

# issue #3710
@test prevind(SubString("{var}",2,4),4) == 3

# issue #4183
@test split(SubString("x", 2, 0), "y") == AbstractString[""]

# issue #6772
@test float(SubString("10",1,1)) === 1.0
@test float(SubString("1 0",1,1)) === 1.0
@test parse(Float32,SubString("10",1,1)) === 1.0f0

# issue #5870
@test !ismatch(Regex("aa"), SubString("",1,0))
@test ismatch(Regex(""), SubString("",1,0))

# isvalid(), chr2ind() and ind2chr() for SubString{DirectIndexString}
let s="lorem ipsum",
    sdict=Dict(SubString(s,1,11)=>s,
               SubString(s,1,6)=>"lorem ",
               SubString(s,1,0)=>"",
               SubString(s,2,4)=>"ore",
               SubString(s,2,16)=>"orem ipsum",
               SubString(s,12,14)=>""
               )
    for (ss,s) in sdict
        local ss
        for i in -1:12
            @test isvalid(ss,i)==isvalid(s,i)
        end
    end
    for (ss,s) in sdict
        local ss
        for i in 1:length(ss)
            @test ind2chr(ss,i)==ind2chr(s,i)
        end
    end
    for (ss,s) in sdict
        local ss
        for i in 1:length(ss)
            @test chr2ind(ss,i)==chr2ind(s,i)
        end
    end
end #let

#for isvalid(SubString{String})
let s = "Σx + βz - 2"
  for i in -1:length(s)+2
      ss=SubString(s,1,i)
      @test isvalid(ss,i)==isvalid(s,i)
  end
end

ss=SubString("hello",1,5)
@test_throws BoundsError ind2chr(ss, -1)
@test_throws BoundsError chr2ind(ss, -1)
@test_throws BoundsError chr2ind(ss, 10)
@test_throws BoundsError ind2chr(ss, 10)

# length(SubString{String}) performance specialization
let s = "|η(α)-ϕ(κ)| < ε"
    @test length(SubString(s,1,0))==length(s[1:0])
    @test length(SubString(s,4,4))==length(s[4:4])
    @test length(SubString(s,1,7))==length(s[1:7])
    @test length(SubString(s,4,11))==length(s[4:11])
end

## Reverse strings ##

rs = RevString("foobar")
@test length(rs) == 6
@test sizeof(rs) == 6
@test isascii(rs)

# issue #4586
@test rsplit(RevString("ailuj"),'l') == ["ju","ia"]
@test parse(Float64,RevString("64")) === 46.0

# reverseind
for T in (String, GenericString)
    for prefix in ("", "abcd", "\U0001d6a4\U0001d4c1", "\U0001d6a4\U0001d4c1c", " \U0001d6a4\U0001d4c1")
        for suffix in ("", "abcde", "\U0001d4c1β\U0001d6a4", "\U0001d4c1β\U0001d6a4c", " \U0001d4c1β\U0001d6a4")
            for c in ('X', 'δ', '\U0001d6a5')
                s = convert(T, string(prefix, c, suffix))
                r = reverse(s)
                ri = search(r, c)
                @test r == RevString(s)
                @test c == s[reverseind(s, ri)] == r[ri]
                s = RevString(s)
                r = reverse(s)
                ri = search(r, c)
                @test c == s[reverseind(s, ri)] == r[ri]
                s = convert(T, string(prefix, prefix, c, suffix, suffix))
                pre = convert(T, prefix)
                sb = SubString(s, nextind(pre, endof(pre)), endof(convert(T, string(prefix, prefix, c, suffix))))
                r = reverse(sb)
                ri = search(r, c)
                @test c == sb[reverseind(sb, ri)] == r[ri]
            end
        end
    end
end

## Cstring tests ##

# issue #13974: comparison against pointers

str = String("foobar")
ptr = pointer(str)
cstring = Cstring(ptr)
@test ptr == cstring
@test cstring == ptr

# convenient NULL string creation from Ptr{Void}
nullstr = Cstring(C_NULL)

# Comparisons against NULL strings
@test ptr != nullstr
@test nullstr != ptr

# Short-hand comparison against C_NULL
@test nullstr == C_NULL
@test C_NULL == nullstr
@test cstring != C_NULL
@test C_NULL != cstring

max_stars_repo_name: baggepinnen/julia

max_issues_repo_name: baggepinnen/julia

max_forks_repo_name: afniedermayer/julia


text:
using CHIP8
using Base.Test

include("constructor.jl")
include("io.jl")
include("emulate.jl")


text:
module TestTSet

using Base.Test
using Transit

a = Transit.TSet([1,2])
a1 = Transit.TSet([1,2])

b = Transit.TSet([3,4])
b1 = Transit.TSet([4,3])

c = Transit.TSet(Any[a, b])
c1 = Transit.TSet(Any[a,b])

s0 = Transit.TSet([Transit.TSymbol("aaa"), Transit.TSymbol("bbb"),  Transit.TSymbol("ccc")])
s1 = Transit.TSet([Transit.TSymbol("aaa"), Transit.TSymbol("bbb"),  Transit.TSymbol("ccc")])

s2 = Transit.TSet([Transit.TSymbol("aaa"), Transit.TSymbol("ccc"),  Transit.TSymbol("bbb")])
s3 = Transit.TSet([Transit.TSymbol("ccc"), Transit.TSymbol("aaa"),  Transit.TSymbol("bbb")])

@test s0 == s0
@test s0 == s1
@test s1 == s2
@test s2 == s1
@test s3 == s2
@test s3 == s3

@test a == a
@test a == a1
@test b == b
@test b == b1
@test a != b
@test b != a

@test c == c
@test c == c1

end


text:
@doc raw"""
    layerBackProp(
        cLayer::FCLayer,
        model::Model,
        actFun::SoS,
        Ao::AbstractArray,
        labels::AbstractArray,
    ) where {SoS<:Union{Type{softmax},Type{σ}}}

For output `FCLayer` layers with softmax and sigmoid activation functions

"""
function layerBackProp(
    cLayer::FCLayer,
    model::Model,
    actFun::SoS,
    Ao::AbstractArray,
    labels::AbstractArray,
) where {SoS<:Union{Type{softmax},Type{σ}}}

    lossFun = model.lossFun
    dlossFun = Symbol("d", lossFun)
    Y = labels
    dZ = eval(:($dlossFun($Ao, $Y)))

    # dAi = cLayer.W'dZ
    return dZ #, dAi
end #softmax or σ layerBackProp

@doc raw"""
    layerBackProp(
        cLayer::Activation,
        model::Model,
        actFun::SoS,
        Ao::AbstractArray,
        labels::AbstractArray,
    ) where {SoS<:Union{Type{softmax},Type{σ}}}

For output `Activation` layers with softmax and sigmoid activation functions
"""
function layerBackProp(
    cLayer::Activation,
    model::Model,
    actFun::SoS,
    Ao::AbstractArray,
    labels::AbstractArray,
) where {SoS<:Union{Type{softmax},Type{σ}}}

    lossFun = model.lossFun
    dlossFun = Symbol("d", lossFun)
    Y = labels
    dZ = eval(:($dlossFun($Ao, $Y)))

    return dZ
end #softmax or σ layerBackProp

@doc raw"""
    layerBackProp(
        cLayer::FCLayer,
        model::Model,
        FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        dAo::AbstractArray{T1,2} = Array{Any,2}(undef,0,0);
        labels::AbstractArray{T2,2} = Array{Any,2}(undef,0,0),
        kwargs...,
    ) where {T1,T2}

Perform the back propagation of `FCLayer` type on the activations and trainable parameters `W` and `B`

# Argument

- `cLayer` := the layer to perform the backprop on

- `model` := the `Model`

- `FCache` := the cache values of the forprop

- `BCache` := the cache values of the backprop from the front `Layer`(s)

- `dAo` := (for test purpose) the derivative of the front `Layer`

- `labels` := in case this is the output layer

# Return

- A `Dict{Symbol, AbstractArray}(:dA => dAi)`
"""
function layerBackProp(
    cLayer::FCLayer,
    model::Model,
    FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    dAo::AbstractArray{T1,2} = Array{Any,2}(undef,0,0);
    labels::AbstractArray{T2,2} = Array{Any,2}(undef,0,0),
    kwargs...,
) where {T1,T2}

    prevLayer = cLayer.prevLayer
    lossFun = model.lossFun

    Ao = FCache[cLayer][:A]

    regulization, λ = model.regulization, model.λ

    actFun = cLayer.actFun

    dZ = []
    if cLayer.actFun == model.outLayer.actFun

        dZ = layerBackProp(cLayer, model, eval(:($actFun)), Ao, labels)

    elseif length(dAo) != 0
        keepProb = cLayer.keepProb
        if keepProb < 1.0 #to save time of multiplication in case keepProb was one
            D = rand(cLayer.channels, 1) .< keepProb
            dAo .*= D
            dAo ./= keepProb
        end

        dActFun = Symbol("d", cLayer.actFun)

        Z = FCache[cLayer][:Z]

        dZ = dAo .* eval(:($dActFun($Z)))

    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dAo = []
        for nextLayer in cLayer.nextLayers
            try
                dAo .+= getLayerSlice(cLayer, nextLayer, BCache)
            catch e #in case first time DimensionMismatch
                dAo = getLayerSlice(cLayer, nextLayer, BCache)
            end
        end #for

        keepProb = cLayer.keepProb
        if keepProb < 1.0 #to save time of multiplication in case keepProb was one
            D = rand(cLayer.channels, 1) .< keepProb
            dAo .*= D
            dAo ./= keepProb
        end

        dActFun = Symbol("d", cLayer.actFun)

        Z = FCache[cLayer][:Z]

        dZ = dAo .* eval(:($dActFun($Z)))

    else #in case not every next layer done backprop
        return Dict(:dA => Array{Any,1}(undef,0))

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    m = size(dZ)[end]

    cLayer.dW = dZ * FCache[cLayer.prevLayer][:A]' ./ m

    if regulization > 0
        if regulization == 1
            cLayer.dW .+= (λ / 2m)
        else
            cLayer.dW .+= (λ / m) .* cLayer.W
        end
    end

    cLayer.dB = 1 / m .* sum(dZ, dims = 2)

    dAi = cLayer.W'dZ

    cLayer.backCount += 1

    return Dict(:dA => dAi)

end #function layerBackProp(cLayer::FCLayer)

@doc raw"""
    layerBackProp(
        cLayer::Activation,
        model::Model,
        FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
        labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
        kwargs...,
    ) where {T1,T2,N}

Perform the back propagation of `Activation` type

# Argument

- `cLayer` := the layer to perform the backprop on

- `model` := the `Model`

- `FCache` := the cache values of the forprop

- `BCache` := the cache values of the backprop from the front `Layer`(s)

- `dAo` := (for test purpose) the derivative of the front `Layer`

- `labels` := in case this is the output layer

# Return

- A `Dict{Symbol, AbstractArray}(:dA => dAi)`
"""
function layerBackProp(
    cLayer::Activation,
    model::Model,
    FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
    labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
    kwargs...,
) where {T1,T2,N}

    prevLayer = cLayer.prevLayer
    lossFun = model.lossFun

    Ao = FCache[cLayer][:A]

    regulization, λ = model.regulization, model.λ

    actFun = cLayer.actFun

    dZ = []
    if cLayer.actFun == model.outLayer.actFun

        dZ = layerBackProp(cLayer, model, eval(:($actFun)), Ao, labels)

    elseif length(dAo) != 0
        keepProb = cLayer.keepProb
        if keepProb < 1.0 #to save time of multiplication in case keepProb was one
            D = rand(cLayer.channels, 1) .< keepProb
            dAo .*= D
            dAo ./= keepProb
        end

        dActFun = Symbol("d", cLayer.actFun)

        Z = FCache[cLayer][:Z]

        dZ = dAo .* eval(:($dActFun($Z)))

    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dAo = []
        for nextLayer in cLayer.nextLayers
            try
                dAo .+= getLayerSlice(cLayer, nextLayer, BCache)
            catch e #in case first time DimensionMismatch
                dAo = getLayerSlice(cLayer, nextLayer, BCache)
            end
        end #for

        dActFun = Symbol("d", cLayer.actFun)

        ## get the input as the output of the previous layer
        Z = FCache[cLayer.prevLayer][:A]

        dZ = dAo .* eval(:($dActFun($Z)))

    else #in case not every next layer done backprop
        return Dict(:dA => Array{Any,1}(undef,0))

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    # cLayer.dW = dZ * FCache[cLayer.prevLayer][:A]' ./ m
    #
    # if regulization > 0
    #     if regulization == 1
    #         cLayer.dW .+= (λ / 2m)
    #     else
    #         cLayer.dW .+= (λ / m) .* cLayer.W
    #     end
    # end
    #
    # cLayer.dB = 1 / m .* sum(dZ, dims = 2)

    dAi = dZ

    cLayer.backCount += 1

    return Dict(:dA => dAi)

end #function layerBackProp(cLayer::Activation

### Flatten

@doc raw"""
    layerBackProp(
        cLayer::Flatten,
        model::Model,
        FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
        labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
        kwargs...,
    ) where {T1,T2,N}

Perform the back propagation of `Flatten` type

# Argument

- `cLayer` := the layer to perform the backprop on

- `model` := the `Model`

- `FCache` := the cache values of the forprop

- `BCache` := the cache values of the backprop from the front `Layer`(s)

- `dAo` := (for test purpose) the derivative of the front `Layer`

- `labels` := in case this is the output layer

# Return

- A `Dict{Symbol, AbstractArray}(:dA => dAi)`
"""
function layerBackProp(
    cLayer::Flatten,
    model::Model,
    FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
    labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
    kwargs...,
) where {T1,T2,N}

    prevLayer = cLayer.prevLayer

    Ao = FCache[cLayer][:A]

    regulization, λ = model.regulization, model.λ

    if length(dAo) == 0 && all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dAo = []
        for nextLayer in cLayer.nextLayers
            try
                dAo .+= getLayerSlice(cLayer, nextLayer, BCache)
            catch e #in case first time DimensionMismatch
                dAo = getLayerSlice(cLayer, nextLayer, BCache)
            end
        end #for

    else #in case not every next layer done backprop
        return Dict(:dA => Array{Any,1}(undef,0))

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    dAi = reshape(dAo, cLayer.inputS)

    cLayer.backCount += 1

    return Dict(:dA => dAi)

end #function layerBackProp(cLayer::Activation

### AddLayer

@doc raw"""
    layerBackProp(
        cLayer::AddLayer,
        model::Model,
        FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
        labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
        kwargs...,
    ) where {T1,T2,N}

Perform the back propagation of `AddLayer` type

# Argument

- `cLayer` := the layer to perform the backprop on

- `model` := the `Model`

- `FCache` := the cache values of the forprop

- `BCache` := the cache values of the backprop from the front `Layer`(s)

- `dAo` := (for test purpose) the derivative of the front `Layer`

- `labels` := in case this is the output layer

# Return

- A `Dict{Symbol, AbstractArray}(:dA => dAi)`
"""
function layerBackProp(
    cLayer::AddLayer,
    model::Model,
    FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
    labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
    kwargs...,
) where {T1,T2,N}

    if length(dAo) != 0
        return Dict(:dA => dAo)

    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dAo = []
        for nextLayer in cLayer.nextLayers
            try
                dAo .+= getLayerSlice(cLayer, nextLayer, BCache)
            catch e
                dAo = getLayerSlice(cLayer, nextLayer, BCache) #to initialize dA
            end #try/catch
        end #for

    else #in case not every next layer done backprop
        return Dict(:dA => Array{Any,1}(undef,0))

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    cLayer.backCount += 1

    return Dict(:dA => dAo)
end #function layerBackProp(cLayer::AddLayer

### ConcatLayer

function layerBackProp(
    cLayer::ConcatLayer,
    model::Model,
    FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
    labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
    kwargs...,
) where {T1,T2,N}

    if length(dAo) != 0
        return Dict(:dA => dAo)

    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dAo = []
        for nextLayer in cLayer.nextLayers
            try
                dAo .+= getLayerSlice(cLayer, nextLayer, BCache)
            catch e
                dAo = getLayerSlice(cLayer, nextLayer, BCache) #to initialize dA
            end #try/catch
        end #for

    else #in case not every next layer done backprop
        return Dict(:dA => Array{Any,1}(undef,0))

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    cLayer.backCount += 1

    return Dict(:dA => dAo)
end #function layerBackProp(cLayer::AddLayer

@doc raw"""
    layerBackProp(
        cLayer::Input,
        model::Model,
        FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
        labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
        kwargs...,
    ) where {T1,T2,N}

Perform the back propagation of `Input` type

# Argument

- `cLayer` := the layer to perform the backprop on

- `model` := the `Model`

- `FCache` := the cache values of the forprop

- `BCache` := the cache values of the backprop from the front `Layer`(s)

- `dAo` := (for test purpose) the derivative of the front `Layer`

- `labels` := in case this is the output layer

# Return

- A `Dict{Symbol, AbstractArray}(:dA => dAi)`
"""
function layerBackProp(
    cLayer::Input,
    model::Model,
    FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    dAo::AbstractArray{T1,N} = Array{Any,1}(undef,0);
    labels::AbstractArray{T2,N} = Array{Any,1}(undef,0),
    kwargs...,
) where {T1,T2,N}

    if length(dAo) != 0
        return Dict(:dA => dAo)
    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dAo = []
        for nextLayer in cLayer.nextLayers
            try
                dAo .+= getLayerSlice(cLayer, nextLayer, BCache)
            catch e
                dAo = getLayerSlice(cLayer, nextLayer, BCache) #need to initialize dA
            end #try/catch
        end #for

    else #in case not every next layer done backprop
        return Dict(:dA => Array{Any,1}(undef,0))

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    cLayer.backCount += 1

    return Dict(:dA => dAo)
end #function layerBackProp(cLayer::Input

@doc raw"""
    layerBackProp(
        cLayer::BatchNorm,
        model::Model,
        FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
        dAo::AbstractArray = Array{Any,1}(undef,0);
        labels::AbstractArray = Array{Any,1}(undef,0),
        kwargs...,
    )

Perform the back propagation of `BatchNorm` type on the activations and trainable parameters `W` and `B`

# Argument

- `cLayer` := the layer to perform the backprop on

- `model` := the `Model`

- `FCache` := the cache values of the forprop

- `BCache` := the cache values of the backprop from the front `Layer`(s)

- `dAo` := (for test purpose) the derivative of the front `Layer`

- `labels` := in case this is the output layer

# Return

- A `Dict{Symbol, AbstractArray}(:dA => dAi)`
"""
function layerBackProp(
    cLayer::BatchNorm,
    model::Model,
    FCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    BCache::Dict{Layer, Dict{Symbol, AbstractArray}},
    dAo::AbstractArray = Array{Any,1}(undef,0);
    labels::AbstractArray = Array{Any,1}(undef,0),
    kwargs...,
)

    prevLayer = cLayer.prevLayer
    lossFun = model.lossFun

    Ao = FCache[cLayer][:A]
    Ai = FCache[cLayer.prevLayer][:A]

    m = size(Ao)[end]

    normDim = cLayer.dim+1

    regulization, λ = model.regulization, model.λ

    N = length(cLayer.outputS)

    dZ = []
    if length(dAo) != 0
        dAo = permutedims(dAo, [N, (1:N-1)...])
        dZ = cLayer.W .* dAo
    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dAo = []
        for nextLayer in cLayer.nextLayers
            try
                dAo .+= getLayerSlice(cLayer, nextLayer, BCache)
            catch e
                dAo = getLayerSlice(cLayer, nextLayer, BCache) #need to initialize dA
            end #try/catch
        end #for
        dAo = permutedims(dAo, [N, (1:N-1)...])
        dZ = cLayer.W .* dAo

    else #in case not every next layer done backprop
        return Dict(:dA => Array{Any,1}(undef,0))

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    #Z is already flipped
    cLayer.dW = sum(dAo .* FCache[cLayer][:Z], dims = 1:normDim)

    if regulization > 0
        if regulization == 1
            cLayer.dW .+= (λ / 2m)
        else
            cLayer.dW .+= (λ / m) .* cLayer.W
        end
    end

    tA = eltype(dAo)

    cLayer.dB = sum(dAo, dims = 1:normDim)

    Num = prod(size(dAo)[1:normDim])

    varep = FCache[cLayer][:var] .+ cLayer.ϵ
    Ai_μ = FCache[cLayer][:Ai_μ]
    Ai_μ_s = FCache[cLayer][:Ai_μ_s]

    svarep = sqrt.(varep)

    dZ1 = dZ ./ svarep

    dZ2 = sum(dZ .* Ai_μ; dims = 1:normDim)
    dZ2 .*= (-1 ./ (svarep .^ 2))
    dZ2 .*= (0.5 ./ svarep)
    dZ2 = dZ2 .* (1.0 / Num) .* ones(promote_type(eltype(Ai_μ), eltype(dZ2)), size(Ai_μ))
    dZ2 .*= (2 .* Ai_μ)

    dẐ = dZ1 .+ dZ2

    dZ3 = dẐ

    dZ4 = .- sum(dẐ; dims = 1:normDim)

    dZ4 = dZ4 .* (1.0 / Num) .* ones(promote_type(eltype(Ai), eltype(dZ4)), size(dZ3))

    dAip = dZ3 .+ dZ4

    dAi = tA.(permutedims(dAip, [(2:N)..., 1]))

    # dẐ =
    #     dZ ./ sqrt.(varep) .*
    #     (.-(Ai_μ_s) .* ((N - 1) / N^2) ./ (varep) .^ 2 .+ 1) .*
    #     (1 - 1 / N)
    #
    # dAi = dẐ

    cLayer.backCount += 1

    return Dict(:dA => dAi)

end #function layerBackProp(cLayer::BatchNorm)

export layerBackProp


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((1, 2),),
    ((2, 1),),
    ]
insts["fast"] = [
    ((2, 3),),
    ((3, 3),),
    ]
insts["various"] = [
    ((2, 2),),
    ((4, 2),),
    ((6, 2),),
    ((2, 4),),
    ((3, 4),),
    ((4, 4),),
    ]
return (RelEntrEntanglementJuMP, insts)


text:
using PERMANOVA, DataFrames,Distances, RCall, BenchmarkTools
using Test
R"library(vegan)"

x = rand(1:4,100)
z = rand(1:4,100)
y = rand(100,5) .+ (x .* (x .==3) .*(z .>2))

preds = [["a","b","c","d"][i] for i in x]
preds2 = [["a","b","c","d"][i] for i in z]
preds3 = [rand(["a","b","c","d"]) for i in x]
df = DataFrame([preds,preds2,preds3],[:X,:Z,:W])
Hyde = hydra(df,y,BrayCurtis,@formula(1~X+Z+W),1000)
R"ad =adonis2($y ~X+Z+W,$df,1000)"
n_terms = 3

@test sum(R"ad$Df" .== Hyde.results[:,"Df"]) == n_terms +2
@test sum(R"ad$SumOfSqs".≈ Hyde.results[:,"SumOfSqs"]) == n_terms +2
@test sum(R"ad$R2".≈ Hyde.results[:,"R2"]) == n_terms +2
@test sum(R"ad$F".≈ Hyde.results[:,"F"]) == n_terms

x = rand(1:4,1000)
z = rand(1:4,1000)
y = rand(1000,5) .+ (x .* (x .==3) .*(z .>2))

preds = [["a","b","c","d"][i] for i in x]
preds2 = [["a","b","c","d"][i] for i in z]
preds3 = [rand(["a","b","c","d"]) for i in x]
df = DataFrame([preds,preds2,preds3],[:X,:Z,:W])
Hyde = hydra(df,y,BrayCurtis,@formula(1~X+Z+W),10000)
R"ad =adonis2($y ~X+Z+W,$df,10000)"
n_terms = 3

@test sum(R"ad$Df" .== Hyde.results[:,"Df"]) == n_terms +2
@test sum(R"ad$SumOfSqs".≈ Hyde.results[:,"SumOfSqs"]) == n_terms +2
@test sum(R"ad$R2".≈ Hyde.results[:,"R2"]) == n_terms +2
@test sum(R"ad$F".≈ Hyde.results[:,"F"]) == n_terms


text:
using Polynomials
using ARMA

"""`noisePSD(tes, freq, [SI_amp])`

Returns (noise,A,B,C,D) where `noise` is a noise power spectral
density in A^2 per Hz, at each input frequency (Hz) given by `freq`.
A,B,C,D are the components of the noise, such that A+B+C+D=noise.
In order, they are TES Johnson noise; load resistor Johnson noise;
thermal fluctuation noise; and amplifier noise.

You can input the amplifier current noise `SI_amp`, in A^2/Hz, or you
can take the default value."""
function noisePSD(tes::IrwinHiltonTES, freq::AbstractVector, SI_amp=5e-22)
    # This term F goes from 0 to one and depends on whether the thermal conductivity is ballistic
    # or diffusive, hardcoded as 1 for now
    F  = 1
    SP_TFN = 4kb*tes.T0^2*tes.G0*F
    SV_TES = 4kb*tes.T0*tes.R0*(1+2*tes.beta) # TES voltage noise
    SV_L   = 4kb*tes.T0*tes.Rl  # Load voltage noise

    omega = 2pi*freq  # Radians / sec
    sIomeg = (1-tes.tauplus/tes.taucc)*(1-tes.tauminus/tes.taucc)./((1 .+im*omega*tes.tauplus).*(1 .+im*omega*tes.tauminus)) /(tes.I0*tes.R0*(2+tes.beta))
    sIomeg2 = abs2.(sIomeg)

    Inoise_TES = SV_TES*tes.I0^2/tes.loopgain^2 * (1 .+(tes.tauthermal*omega).^2) .* sIomeg2
    Inoise_load = SV_L*tes.I0^2*(tes.loopgain-1)^2/tes.loopgain^2 * (1 .+(tes.taucc*omega).^2) .* sIomeg2
    Inoise_TFN = SP_TFN*sIomeg2
    Inoise_amp = SI_amp

    Inoise = Inoise_TFN .+ Inoise_amp .+ Inoise_TES .+ Inoise_load
    Inoise, Inoise_TES, Inoise_load, Inoise_TFN, Inoise_amp .+ zeros(Float64, length(Inoise))
end

noisePSD(tes::BiasedTES, freq::AbstractVector, SI_amp::Float64=5e-22) =
    noisePSD(IrwinHiltonTES(tes), freq, SI_amp)

"""`NoiseModel(tes, sampleTime, [SI_amp])`

Returns noise of an IrwinHiltonTES modeled as an ARMA(2,2) process, specifically as an
ARMAModel object from package ARMA.jl"""
function NoiseModel(tes::IrwinHiltonTES, sampleTime::Float64, SI_amp=5e-22)
    # This term F goes from 0 to one and depends on whether the thermal conductivity is ballistic
    # or diffusive, hardcoded as 1 for now
    F  = 1
    SP_TFN = 4kb*tes.T0^2*tes.G0*F
    SV_TES = 4kb*tes.T0*tes.R0*(1+2*tes.beta) # TES voltage noise
    SV_L   = 4kb*tes.T0*tes.Rl  # Load voltage noise

    A = SI_amp
    B = abs2((1-tes.tauplus/tes.taucc)*(1-tes.tauminus/tes.taucc) /(tes.I0*tes.R0*(2+tes.beta)))
    C = SP_TFN
    D = SV_TES*tes.I0^2/tes.loopgain^2
    E = SV_L*tes.I0^2*(tes.loopgain-1)^2/tes.loopgain^2

    # (u,v) are the solutions to the quadratic equation ax^2+bx+c=0
    a = A * (tes.tauplus*tes.tauminus)^2
    b = -(A*(tes.tauplus^2+tes.tauminus^2)+B*(D*tes.tauthermal^2+E*tes.taucc^2))
    c = A+B*(C+D+E)
    discr = sqrt(b^2-4a*c)
    u2 = (-b-discr)/(2a)
    v2 = (-b+discr)/(2a)

    u,v = sqrt(u2), sqrt(v2)
    K = sqrt(a)
    theta = [(sampleTime*u*.5+1)*(sampleTime*v*0.5+1),
            sampleTime^2*u*v*.5-2,
            (sampleTime*u*.5-1)*(sampleTime*v*.5-1)]
    phi = [(2tes.tauplus/sampleTime+1)*(2tes.tauminus/sampleTime+1),
            2-8tes.tauplus*tes.tauminus/sampleTime^2,
            (1-2tes.tauplus/sampleTime)*(1-2tes.tauminus/sampleTime)]

    # Rescale so that phi has unit coefficients for the z^0 term; theta[1] is sigma.
    sigma = theta[1]/phi[1]*4K/(sampleTime^2)
    sigma /= sqrt(2sampleTime)
    theta .*= sigma/theta[1]
    phi ./= phi[1]

    ARMAModel(real(theta), real(phi))
end

NoiseModel(tes::BiasedTES, sampleTime::Float64, SI_amp::Float64=5e-22) =
    NoiseModel(IrwinHiltonTES(tes), sampleTime, SI_amp)


text:
@testset "$C" for C in (HSV, HSL, HSI)
    @testset "$C over $C: $T" for T in (Float64, Float32)
        @test BlendNormal(    C{T}(100,  0.75,  0), C{T}(0, 0.5, 1)) ≈ C{T}(0, 0.5, 1)
        @test BlendHue(       C{T}(100,  0.75,  0), C{T}(0, 0.5, 1)) ≈ C{T}(0, 0.75, 0)
        @test BlendSaturation(C{T}(100,  0.75,  0), C{T}(0, 0.5, 1)) ≈ C{T}(100, 0.5, 0)
        @test BlendColor(     C{T}(100,  0.75,  0), C{T}(0, 0.5, 1)) ≈ C{T}(0, 0.5, 0)
        @test BlendLuminosity(C{T}(100,  0.75,  0), C{T}(0, 0.5, 1)) ≈ C{T}(100, 0.75, 1)
    end

    A = coloralpha(C)
    @testset "$A over $C: $T" for T in (Float64, Float32)
        @test BlendNormal(C{T}(100,  0.75,  0), A{T}(0, 0.5, 1, 0.6)) ≈ C{T}(40, 0.6, 0.6)
    end

    @testset "$C over $A: $T" for T in (Float64, Float32)
        @test BlendNormal(A{T}(100,  0.75, 0, 0.6), C{T}(0, 0.5, 1)) ≈ A{T}(0, 0.5, 1, 1)
    end

    @testset "$A over $A: $T" for T in (Float64, Float32)
        @test BlendNormal(A{T}(100,  0.75, 0, 0.6), A{T}(0, 0.5, 1, 0.6)) ≈ A{T}(200/7, 4/7, 5/7, 0.84)
    end

    @testset "$C over gray $C: $T" for T in (Float64, Float32)
        @test blend(C{T}(100, 0, 1), C{T}(200, 1, 0.5), opacity=0.5) ≈ C{T}(200, 0.5, 0.75)
    end
end


text:
function build_new(obj::Region, arg0::BSPTree)
    return jcall(obj, "buildNew", Region, (BSPTree,), arg0)
end

function check_point(obj::Region, arg0::Point)
    return jcall(obj, "checkPoint", Region_Location, (Point,), arg0)
end

function contains(obj::Region, arg0::Region)
    return jcall(obj, "contains", jboolean, (Region,), arg0)
end

function copy_self(obj::Region)
    return jcall(obj, "copySelf", Region, ())
end

function get_barycenter(obj::Region)
    return jcall(obj, "getBarycenter", Point, ())
end

function get_boundary_size(obj::Region)
    return jcall(obj, "getBoundarySize", jdouble, ())
end

function get_size(obj::Region)
    return jcall(obj, "getSize", jdouble, ())
end

function get_tree(obj::Region, arg0::jboolean)
    return jcall(obj, "getTree", BSPTree, (jboolean,), arg0)
end

function intersection(obj::Region, arg0::SubHyperplane)
    return jcall(obj, "intersection", SubHyperplane, (SubHyperplane,), arg0)
end

function is_empty(obj::Region)
    return jcall(obj, "isEmpty", jboolean, ())
end

function is_empty(obj::Region, arg0::BSPTree)
    return jcall(obj, "isEmpty", jboolean, (BSPTree,), arg0)
end

function is_full(obj::Region)
    return jcall(obj, "isFull", jboolean, ())
end

function is_full(obj::Region, arg0::BSPTree)
    return jcall(obj, "isFull", jboolean, (BSPTree,), arg0)
end

function project_to_boundary(obj::Region, arg0::Point)
    return jcall(obj, "projectToBoundary", BoundaryProjection, (Point,), arg0)
end