text:
# converts exhaustive xrd data to xrd gradients
using LinearAlgebra
using GaussianDistributions
using SARA
using SARA: TemperatureProfile, stripe_to_global, get_relevant_T, iSARA
using Base.Threads: @threads
using SARA: iSARA, oSARA
using JSON
using HDF5

path = "../inner/"
# path = "SARA/NatCom2020/inner/"
include(path * "inner_estimate_noise.jl")
include(path * "inner_load_data.jl")

datapath = "SARA/NatCom2020/outer/data/"
# datapath = ".."
file = "Bi2O3_19F44_01_outer_xrd_data.h5"
f = h5open(datapath * file)
positions = read(f, "positions")
temperatures = read(f, "temperatures")
dwelltimes = read(f, "dwelltimes")
spectrograms = read(f, "spectrograms")
close(f)
conditions = tuple.(temperatures, dwelltimes)

########### match rescaling parameters from optical coefficients to xrd data
# returns true if dwell time and peak temperature are the same
same_conditions(f) = g->same_conditions(f, g)
same_conditions(f, g) = all(f .== g)

innerpath = "SARA/NatCom2020/inner/"
# path = ".."
file = "Bi2O3_19F44_01_inner_loop_data.json"
f = JSON.parsefile(innerpath * "data/" * file)
_, _, optical_rescaling, optical_temperatures, optical_dwelltimes = load_data(f)
optical_conditions = tuple.(optical_temperatures, optical_dwelltimes)

rescaling_parameters = []
for i in eachindex(conditions)
    j = findfirst(same_conditions(conditions[i]), optical_conditions)
    if !isnothing(j)
        push!(rescaling_parameters, optical_rescaling[j])
    end
end

# 2. dimensionality reduction
using SARA: legendre_coefficients, svd_coefficients
ncoeff = 16
coefficients = svd_coefficients(spectrograms, ncoeff)
# coefficients = legendre_coefficients(spectrograms, ncoeff)
coefficients ./= maximum(coefficients)

########### estimating noise variance of normalized coefficient
i = 3
x = positions[:, i]
y = coefficients[1, :, i]
using CovarianceFunctions
const Kernel = CovarianceFunctions

l = .1 # from inner loop benchmarks
θ = rescaling_parameters[i]
k = iSARA(l, θ)
σ = estimate_noise(x, y, k)
println(σ)

# using Plots
# plotly()
# C = Gaussian(k) | (x, y, σ^2)
# plot(x, mean(C, x), ribbon = 2std(C).(x))
# scatter!(x, y)

############## convert to global data
# dT = 20 # uncertainty in peak temperature in C
# dx = 10/1000 # uncertainty in position in mm
dT = 25 # uncertainty in peak temperature in C
dx = 50/1000 # uncertainty in position in mm
P = TemperatureProfile(dT, dx)
nout = 32 # we can still change this after the creation of the data
# constant_offset = Val(false)
# T_proportions = (.75, .99) # as a proportion of T_max, what data to convert to gradients
constant_offset_bool = true
constant_offset = Val(constant_offset_bool)
T_offset = [200, 10] # number of degrees from T_max we are generating data for
relevant_T = get_relevant_T(constant_offset, T_offset..., nout)

input_noise_bool = true
input_noise = Val(input_noise_bool)
nstripes = size(coefficients, 3)
gradients = zeros(nout, nstripes)
var_gradients = zeros(nout, nstripes)
@threads for i in 1:nstripes
    println(i)
    pos = @view positions[:, i]
    coeff = [@view coefficients[j, :, i] for j in 1:size(coefficients, 1)]
    kernel = iSARA(l, rescaling_parameters[i])
    gradients[:, i], var_gradients[:, i] = stripe_to_global(pos, coeff,
                                            σ, kernel,
                                            P, temperatures[i], dwelltimes[i],
                                            relevant_T, input_noise)
end

outer_temperatures = zeros(nout, nstripes)
outer_dwelltimes = zeros(nout, nstripes)
for i in 1:nstripes
    outer_dwelltimes[:, i] .= dwelltimes[i]
    outer_temperatures[:, i] = relevant_T(temperatures[i], dwelltimes[i])
end

savepath = "SARA/NatCom2020/outer/data/"
# savefile = "Bi2O3_19F44_01_outer_xrd_gradients_large_input_noise.h5"
savefile = "Bi2O3_19F44_01_outer_xrd_gradients_l_1_input_noise_50_new_profile.h5"
# savefile = "Bi2O3_19F44_01_outer_xrd_gradients_no_input_noise.h5"
f = h5open(savepath * savefile, "w")
f["dwelltimes"] = outer_dwelltimes
f["temperatures"] = outer_temperatures
f["gradients"] = gradients
f["var(gradients)"] = var_gradients
f["input noise"] = input_noise_bool
f["nstripes"] = nstripes
f["constant offset"] = constant_offset_bool
f["T_offset"] = T_offset
f["coefficients"] = coefficients
f["coefficient type"] = "svd" # "legendre"
f["nout"] = nout
f["l"] = l
f["coefficient noise variance"] = σ^2
f["temperature uncertainty parameters"] = [dT, dx]
close(f)

# using GaussianDistributions: optimize
# k = (x, y) -> Kernel.NN(100.)(x./l, y./l)
# kernel = Kernel.Lengthscale(Kernel.MaternP(2), l)
# k = (θ) -> ((x, y) -> exp(θ[1]) * Kernel.NN(exp(θ[2]))(x./exp(θ[3]), y./exp(θ[3])))
# θ = zeros(3)
# println(θ)
# θ = optimize(k, θ, x, y, σ^2)
# println(θ)
# k = k(θ)


text:
# This file was generated, do not modify it. # hide
r_maxD = range(rfr, :n_trees, lower=9, upper=15)
r_samF = range(rfr, :sampling_fraction, lower=0.6, upper=0.8)
r = [r_maxD, r_samF];

ext: jl


text:
using BinaryBuilder

cuda_version = v"11.0.3"   # NOTE: could be less specific

sources_linux_x64 = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.1.0.77/11.2_20210127/cudnn-11.2-linux-x64-v8.1.0.77.tgz",
                  "dbe82faf071d91ba9bcf00480146ad33f462482dfee56caf4479c1b8dabe3ecb")
]
sources_linux_ppc64le = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.1.0.77/11.2_20210127/cudnn-11.2-linux-ppc64le-v8.1.0.77.tgz",
                  "0d3f8fa21959e9f94889841cc8445aecf41d2f3c557091b447313afb43034037")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.1.0.77/11.2_20210127/cudnn-11.2-windows-x64-v8.1.0.77.zip",
                  "f5a23402ebed5b2add2ae0dc79df42f055f5851bbb3b1d0f6aa9dcdeeb9346ce")
]

include("../common.jl")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
@defcomp TotalForcing begin
    f_CO2forcing = Parameter(index=[time], unit="W/m2")
    f_CH4forcing = Parameter(index=[time], unit="W/m2")
    f_N2Oforcing = Parameter(index=[time], unit="W/m2")
    f_lineargasforcing = Parameter(index=[time], unit="W/m2")
    exf_excessforcing = Parameter(index=[time], unit="W/m2")

    ft_totalforcing = Variable(index=[time], unit="W/m2")

    function run_timestep(p, v, d, tt)

        # From equation 16 of Hope (2006)
        v.ft_totalforcing[tt] = p.f_CO2forcing[tt] + p.f_CH4forcing[tt] + p.f_N2Oforcing[tt] + p.f_lineargasforcing[tt] + p.exf_excessforcing[tt]
    end
end


text:
module WebSocketUpgrade

using DandelionWebSockets.UniformResourceIdentifiers
using Sockets
using MbedTLS

export ResponseParser, dataread, hascompleteresponse, parseresponse, findheaderfield
export websocketupgrade

const HeaderList = Vector{Pair{String, String}}

struct Request
    host::String
    abs_path::String
    headers::HeaderList
end

function Base.write(io::IO, req::Request)
    write(io, "GET $(req.abs_path) HTTP/1.1\r\n")
    write(io, "Host: $(req.host)\r\n")
    for (k, v) in req.headers
        write(io, "$(k): $(v)\r\n")
    end
    write(io, "\r\n")
end

function findfirstsubstring(needle::AbstractVector{UInt8}, haystack::AbstractVector{UInt8}) :: Union{Int, Nothing}
    lastpossibleindex = length(haystack) - length(needle) + 1
    for i = 1:lastpossibleindex
        substring = haystack[i:i + length(needle) - 1]
        if substring == needle
            return i
        end
    end
    nothing
end

struct ResponseParser
    data::Vector{UInt8}

    ResponseParser() = new(Vector{UInt8}())
end

struct HTTPVersion
    major::Int
    minor::Int
end
struct HTTPResponse
    status::Int
    reasonphrase::String
    headers::HeaderList
    httpversion::HTTPVersion
end

function findheaderfield(response::HTTPResponse, name::String) :: Union{String, Nothing}
    for h in response.headers
        if h[1] == name
            return h[2]
        end
    end
end

struct BadHTTPResponse <: Exception end
struct InvalidHTTPResponse <: Exception end

hastoken(response::HTTPResponse, fieldname::String, token::String) = occursin(token, findheaderfield(response, fieldname))

function validateresponse(response::HTTPResponse) :: HTTPResponse
    hasupgrade = findheaderfield(response, "Upgrade") != nothing
    hasconnection = findheaderfield(response, "Connection") != nothing
    if hasupgrade && !hasconnection
        @error "Found Upgrade header, but not Connection header."
        throw(InvalidHTTPResponse())
    end

    if hasupgrade && hasconnection
        hasconnectionupgrade = hastoken(response, "Connection", "Upgrade")
        if !hasconnectionupgrade
            @error "Connection header does not have the token 'upgrade'"
            throw(InvalidHTTPResponse())
        end
    end

    response
end

dataread(parser::ResponseParser, data::AbstractVector{UInt8}) = append!(parser.data, data)
hascompleteresponse(parser::ResponseParser) = findfirstsubstring(b"\r\n\r\n", parser.data) != nothing
function parseresponse(parser::ResponseParser)
    endofheader = findfirstsubstring(b"\r\n\r\n", parser.data)
    headerbytes = parser.data[1:endofheader-1]
    header = String(headerbytes)
    headerlines = split(header, "\r\n")
    statusline = headerlines[1]

    statuslinematch = match(r"^HTTP/([0-9]+)\.([0-9]+) +([0-9]+) +([^\r\n]*)", statusline)
    if statuslinematch != nothing
        status = parse(Int, statuslinematch.captures[3])
        reasonphrase = statuslinematch.captures[4]
        major = parse(Int, statuslinematch.captures[1])
        minor = parse(Int, statuslinematch.captures[2])

        headers = []
        headerregex = r"([^:]+): (.*)"
        for line in headerlines[2:end]
            m = match(headerregex, line)
            push!(headers, m.captures[1] => m.captures[2])
        end

        #return validateresponse(HTTPResponse(status, reasonphrase, headers, HTTPVersion(major, minor)))
        # Skip the validation until it conforms to requirements.
        return HTTPResponse(status, reasonphrase, headers, HTTPVersion(major, minor))
    end

    throw(BadHTTPResponse())
end

function tlshandshake(socket::TCPSocket)
    entropy = MbedTLS.Entropy()
    rng = MbedTLS.CtrDrbg()
    MbedTLS.seed!(rng, entropy)

    ctx = MbedTLS.SSLContext()
    conf = MbedTLS.SSLConfig()

    MbedTLS.config_defaults!(conf)
    MbedTLS.authmode!(conf, MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED)
    MbedTLS.rng!(conf, rng)

    function show_debug(level, filename, number, msg)
        @show level, filename, number, msg
    end

    MbedTLS.dbg!(conf, show_debug)

    MbedTLS.ca_chain!(conf)

    MbedTLS.setup!(ctx, conf)
    MbedTLS.set_bio!(ctx, socket)

    MbedTLS.handshake(ctx)

    ctx
end

function websocketupgrade(suri::String, headers::HeaderList)
    uri = URI(suri)
    socket = connect(uri.host, uri.port)
    if uri.issecure
        socket = tlshandshake(socket)
    end
    upgraderequest = Request(uri.host, uri.abs_path, headers)
    write(socket, upgraderequest)
    isheadercomplete = false
    parser = ResponseParser()
    while !isheadercomplete
        if eof(socket)
            throw(EOFError())
        end
        data = readavailable(socket)
        dataread(parser, data)
        isheadercomplete = hascompleteresponse(parser)
    end
    # TODO: No excess returned
    response = parseresponse(parser)
    (response, socket)
end

end


text:
using ElTopo
using GeometryTypes

t = ( 1 + sqrt( 5 ) ) / 2;

vertices = Point{3,Float64}[
    [ -1,  t,  0 ], [  1, t, 0 ], [ -1, -t,  0 ], [  1, -t,  0 ],
    [  0, -1,  t ], [  0, 1, t ], [  0, -1, -t ], [  0,  1, -t ],
    [  t,  0, -1 ], [  t, 0, 1 ], [ -t,  0, -1 ], [ -t,  0,  1 ]
]

faces = Face{3,Int64}[
    [1, 12, 6], [1, 6, 2], [1, 2, 8], [1, 8, 11], [1, 11, 12], [2, 6, 10], [6, 12, 5],
    [12, 11, 3], [11, 8, 7], [8, 2, 9], [4, 10, 5], [4, 5, 3], [4, 3, 7], [4, 7, 9],
    [4, 9, 10], [5, 10, 6], [3, 5, 12], [7, 3, 11], [9, 7, 8], [10, 9, 2]
]

msh = HomogenousMesh(vertices ./ sqrt(1+t^2),faces)

par = SurfTrack(allow_vertex_movement=true)
mshnew = stabilize(msh,par)


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# ╔═╡ 10ebd434-adba-11eb-048f-2d084049d48f
html"""

Section 3.5

 How to collaborate on software

Lecture Video

"""

# ╔═╡ ef4eea24-bc1c-43be-b9b6-d073ac2433cf
md"""

# How to collaborate on software

What is the simplest way to collaborate on project online? You can _email files to each other!_ This works for small projects, but at some point, you will need something more suited to the task. For example:

- Encyclopedia -- a wiki system
- School essay -- Google Drive
- **Software project -- GitHub**

For software projects, people of use GitHub, or similar (git-based) tools to collaborate on their work. Why is that?

"""

# ╔═╡ 69b8490f-cf04-4e73-bc7b-639b1fc0e2d6
TableOfContents(; aside=false, depth=1)

# ╔═╡ cbe5fcba-3ed7-41a6-8932-2693e370c16c
md"""

# Why not Google Drive?
"""

# ╔═╡ 812002d3-8603-4ffa-8695-2b2da7f0766a
html"""

If you have not used Google Drive before, here is a small demonstration:

<video src="https://user-images.githubusercontent.com/6933510/117038375-d8497080-ad07-11eb-8260-34e96414131a.mov" data-canonical-src="https://user-images.githubusercontent.com/6933510/117038375-d8497080-ad07-11eb-8260-34e96414131a.mov" controls="controls" muted="muted" class="d-block rounded-bottom-2 width-fit" style="max-height:640px;"></video>
""" |> aside

# ╔═╡ 56db20a5-0e8d-4d34-b3ba-c3ab1b8de58e
md"""
Google Drive (or Dropbox, Nextcloud, etc) is a platform for _realtime collaboration_ on files. Besides synchronizing single files (like the video above), Google Drive also allows you to synchronize entire folders between multiple computers and the cloud. That's awesome! So do people use it for collaborative software projects?

Well, generally not. Software projects have some unique properties that require different tools. **As today's introduction into _git_, let's look at the differences between a software project and a non-software project**, and why platforms like Google Drive don't work well for collaborating on code.

Note that we will take a different perspective that other introductions to _git_. Let's start out with a realtime-synchronization system like Google Drive as our default, and think about why we might want to change its functionality to make it more suitable for software.

"""

# ╔═╡ 83d0162f-5960-4938-8353-91c4cd220459
md"""

## Reason 1: sensitive to change

> **Reason 1:** In software projects, small changes often make a big difference. **Changes to one file affect the entire program.**

For example, here is what Pluto.jl looks like after adding **a single character** to the source code:

"""

# ╔═╡ 51448106-3e0e-4abf-84fc-7e6e81425d12
flex(
	md"""
	#### Before:

	![Schermafbeelding 2021-05-04 om 18 51 27](https://user-images.githubusercontent.com/6933510/117040056-c10b8280-ad09-11eb-9384-d211156440b1.png)
	""",
	md"""
	#### After a single change:

	![Schermafbeelding 2021-05-04 om 18 50 00](https://user-images.githubusercontent.com/6933510/117039958-a46f4a80-ad09-11eb-90fa-c1264d896648.png)
	"""
	) |> outline

# ╔═╡ 3c0059b7-7a99-4c1c-a2bf-e47bdd06a252
outline(x) = Div(x, Dict(
		"border" => "3px solid rgba(0,0,0,.3)",
		"border-radius" => "3px",
		"padding" => "5px",
		))

# ╔═╡ 604837c5-b017-4d6c-a5c5-dab50d5f3f61
md"""
While a single character can mean disaster, this _sensitivity_ also means that programming languages are quite powerful! A small change will often do the trick.
"""

# ╔═╡ 8678b5e1-0b67-4097-82ba-0daa4e878032
md"""
### Synchronizing every change

This sensitivity to changes makes realtime synchronization unsuitable for software projects. To see why, let's say that I want to change `sqrt` to `log` in the following code:
"""

# ╔═╡ 0bf51d8c-7adf-4231-8b2d-c976484a3e7c
Div(
	[md"""
	```julia
	function height(p)
		c1 * sqrt(p * c2)
	end
	```
	""",
	md"to",
	md"""
	```julia
	function height(p)
		c1 * log(p * c2)
	end
	```
	"""],
	Dict(
		"display" => "flex",
		"justify-content" => "space-evenly",
		"align-items" => "center",
	))

# ╔═╡ 02905480-8864-4e56-af3a-6c7c0789ce6f

# ╔═╡ 2ef5db21-0092-4523-b930-0ec99c459ffa
md"""
If you would synchronize every change while making that edit (removing the letters `s`, `q`, `r`, `t`, typing `l`, `o`, `g`), then you are **publishing broken code** to your online code repository. If you were sharing changes in realtime (like Google Drive), then other people would get those broken states while working on their features.
"""

# ╔═╡ 51114dc9-cb32-4d31-b780-6f5e372f8763
md"""
##### Synchronizing every _keystroke_
"""

# ╔═╡ 8b24ce23-4ead-4fbc-875e-a8261f671abe
grid([
	nothing                  md"`sqrt`" md"`sq`" md"` `" md"`lo`" md"`log`"
	"Your computer (local):" pass fail fail fail pass
	"Online (remote): " pass fail fail fail pass
		])

# ╔═╡ 3c8abdfc-cf68-4f10-a3e9-08d24803535b

# ╔═╡ dfb316ca-0502-4419-93c2-6d455b7b2f98
md"""
However, if we only publish code _manually_, we can make sure that we never publish broken code. This means that my collaborators, who might be working on the project at the same time, will not have to deal with errors that are coming from my work-in-progress changes.
"""

# ╔═╡ ae28cad6-2ffe-48b4-895b-fab1bd2f2443
md"""
##### Synchronizing manually
"""

# ╔═╡ 86e4b0dd-9bf5-4711-bf3f-b55ed2627a03
grid([
	nothing                  md"`sqrt`" md"`sq`" md"` `" md"`lo`" md"`log`"
	md"Your computer _(local)_:" pass fail fail fail pass
	md"Online _(remote)_: " pass nothing nothing nothing pass
		])

# ╔═╡ 8ea1ca6b-4bb7-4f53-907d-0e5ca83e5761
md"""
[^sidenote]:

    _Language design side note:_ There _are_ languages that are designed to be robust against small changes, and you can actually use realtime collaboration with those! Examples are [glitch.com](glitch.com) for collaborative HTML and CSS editing _(these languages can ignore syntax errors and continue)_, and the more experimental language [_Dark_](https://darklang.com/) _(which uses a special editor that does not allow you to type errors)_.
""" |> aside

# ╔═╡ ee5a3219-4547-4a9d-b527-3489a2925f68
md"""
### Even larger chunks

In the previous example, we found that we can avoid publishing _invalid syntax_ by not synchronizing on every keystroke like Google Docs, but by synchronizing manually.

We also don't want to synchronize the file every time it is _saved_, because you often try multiple things before finding something that works.
"""

# ╔═╡ e5f2b4e1-f05a-4457-b2e9-091a86ef86bb
md"""
##### Synchronizing every _save_
"""

# ╔═╡ efbd6d1c-d2c9-48a9-9d6c-dc9ce6af0b5b
grid([
	nothing                  md"`sqrt`" md"Let's try `sin`" md"Let's try `cos`" md"Try `log`"
	"Your computer (local):" pass fail fail pass
	"Online (remote): " pass fail fail pass
		])

# ╔═╡ 1421f832-9f97-4ec3-b967-64618983349b
md"""
##### Synchronizing manually
"""

# ╔═╡ b139e8ea-88ad-4df4-9f13-c867edfc2db0
grid([
	nothing                  md"`sqrt`" md"Let's try `sin`" md"Let's try `cos`" md"Try `log`"
	"Your computer (local):" pass fail fail pass
	"Online (remote): " pass nothing nothing pass
		])

# ╔═╡ e5bab2e1-be9e-4654-844a-d50285e330c8
pass2 = html"✅"

# ╔═╡ 47c67fa6-1490-4134-9e0c-5754cb273d1e
pass = "✅" |> HTML

# ╔═╡ 79ab206d-2140-4c78-8fd4-a874fe2551e1
fail = "❌" |> HTML

# ╔═╡ 4ce38ec8-4084-4fe8-a248-87408cd0d39b
bigbreak

# ╔═╡ 38261f04-f410-440e-a1cf-218fa240a0ae
md"""

### Reason 2: _branches_ and _forks_

> **Reason 2:** It should be easy to split a project into multiple _branches_: **copies of the codebase** that can be worked on independently. And after some time, we want to **compare branches** and merge the new additions into one.

Something that you may have done: create a copy of your `presentation.ppt` called `presentation-with-pictures.ppt` before adding your pictures, and if it didn't work, you can go back to yesterday's backup.

Now that's great, but what if your teammate also made a copy `presentation-with-titlepage.ppt` where they also made a change? In this situation, you have created two _branches_ of the same file in parallel, and it will be **difficult to combine those changes** back into a single `presentation-this-one-really.ppt` file.

"""

# ╔═╡ 93ce9618-5484-4572-97f1-1cb9c9367989
bigbreak

# ╔═╡ a0a97cd2-838c-4a2c-9233-969b3274764c
md"""
### Reason 3: automation

> **Reason 3:** Git is a platform for **automation**, and can be a mechanism to automatically test and review changes, release weekly versions, etc etc.

Programmers _love_ to automate things, and git allows all sorts of automation. At first, this can be **intimidating**, and it will take years of experience to learn all the possible tricks. Someone who does 'git magic' as a profession is called a _[DevOps](https://en.wikipedia.org/wiki/DevOps) engineer_. Don't worry about this too much when getting started with GitHub, and ask others for help if necessary.
"""

# ╔═╡ 0faa4042-42f5-4c74-9270-fbf2205920ca
md"""
##### Automatic testing

While we don't recommend learning about _git automation_ right away, one thing you will probably encounter when contributing to an open source project is **automatic testing**. More on this later!
"""

# ╔═╡ 167cbd94-851f-4147-9b9d-9b46c7959785
bigbreak = html"""

"""

# ╔═╡ e3fe649e-2733-455a-897d-d4d2d70b9dc6
md"""
# Our first contribution: _documentation_

Today will be a kick-start into contributing to open source projects! The easiest way to make a great contribution to packages is to contribute a _documentation improvement_. For our example today, we will look at the package [ClimateMARGO.jl](https://github.com/ClimateMARGO/ClimateMARGO.jl), a climate-economic model.

This package is hosted on GitHub, which means that all code and previous versions of code is available at [github.com/ClimateMARGO/ClimateMARGO.jl](https://github.com/ClimateMARGO/ClimateMARGO.jl). The authors also use the same GitHub site to manage their TO-DOs — called [Issues](https://github.com/ClimateMARGO/ClimateMARGO.jl/issues) — and to accept contributions from others — called [Pull Requests](https://github.com/ClimateMARGO/ClimateMARGO.jl/pulls). Today, we are going to make a Pull Request!

"""

# ╔═╡ 2dd5d359-eb44-4696-bc83-f257eff95618
md"""
### What is a Pull Request?

Julia packages are _open source_, which means that you are allowed (and encouraged) to make your own copy of the source code (we call this a _**fork**_) and modify it.

If you made some good changes on your fork then you can suggest to changes to the original project by submitting your code as a _**Pull Request (PR)**_. This is a feature built-in to GitHub, and this is how most code collaboration happens online.
"""

# ╔═╡ 3229f0aa-54d1-43f4-a6aa-cf212620ae13
md"""
> Video tutorial will be available here after the live lecture.
"""

# ╔═╡ bbb59b45-de20-41c1-b5dd-079fc51655f2
md"""
## _After_ the Pull Request

Submitting your Pull Request is a great feeling, but it's not the last step! After your submission, the PR will be reviewed by the project authors, and they will leave their feedback. It is common that a PR is not _merged_ (accepted) directly, but that you are asked to make some small changes.

This is possible because a Pull Request is not static — you can make changes to the PR after submitting it! Remember that a PR is _from one branch to another_. You can make changes to your PR simply by adding additional commits to the _from_ branch.

"""

# ╔═╡ 7af9e69c-2b81-4a90-861c-ed737a4a9ec4
md"""
[^1]:
    If you are working on a _fork_, then by creating a PR, you also give the original project authors access to make changes to the _from_ branch of the PR. This is a useful feature, allowing you to work together on the same codebase.
""" |> aside

# ╔═╡ b83ade3d-6f8d-4ac8-9255-956d0a348416
md"""
As an example of a Pull Request in the wild, let's look at this PR to `JuliaLang/julia`: [github.com/JuliaLang/julia/pull/40596](https://github.com/JuliaLang/julia/pull/40596). The code change is not so important to us right now, but pay attention to the communication around the code change.

"""

# ╔═╡ 4d303641-0299-44d7-ba74-6daea0026b09
md"""
# How to use GitHub

Enough talk, let's create our first repository!
"""

# ╔═╡ 2c456900-047d-432c-a62a-87f8eeaba8d5
md"""
##### Requirements

To follow this introduction, you will need:
- Create an account at [github.com](github.com). Add a friendly profile picture!
- For Windows & MacOS users, download [GitHub Desktop](https://desktop.github.com/), a graphical program to easily manage git repositories on your disk. Linux users can use the slightly more advanced program [GitKraken](https://www.gitkraken.com/).
- A code editor. If you don't have a code editor yet, or if you are using Notepad, gedit or Notepad++, we highly recommend [VS Code](https://code.visualstudio.com/), an open source, beginner-friendly programming environment. _It may seem like something new and difficult to learn, but it is easy to use, and because it **understands** the code that you are writing, it will help you as a beginner programmer._
"""

# ╔═╡ dc112303-061e-4e53-8f58-cf9ea0f556f1
md"""
> Video tutorial will be available here after the live lecture.
"""

# ╔═╡ b2e49cd5-49d5-4ac7-a3ae-9820a97720fb
[
	@htl("remote") => md"The version that is on the internet, you browse the remote on github.com.",
	@htl("local") => md"What's on your computer. Use `pull`, `commit` and `push` to synchronize this with the remote. Google Drive does this all the time."
	] |> vocabulary

# ╔═╡ d400d538-4f73-4490-ad68-aedbb57cd70f

# ╔═╡ a98993b9-a5c0-4260-b96e-2655c472ccba
[
	"fetch" => md"Make your local git aware of any changes online. Do this often!",
	"pull" => md"Apply any changes on the remote version to your local copy. This will get the two _in sync_. Do this often!",
	"commit" => md"Create a collection of changes to files, ready to be `push`ed.",
	"push" => md"Publish any local `commit`s to the remote version.",
	] |> vocabulary

# ╔═╡ 67cf205a-3d89-4cd9-ab5e-febc85ea8af1
md"""
## ⛑ Git trouble!

Sooner or later, you will find yourself in _git trouble_. Unlike Google Drive, it can happen that synchronization requires manual intervention, for example, git might say:

#### `🤐 merge conflict`

Which means that someone changed a file right before you were going to commit.

"""

# ╔═╡ c4442667-072a-4e17-94a4-104a8ec33bd0

# ╔═╡ d8c5273f-7ebc-4399-b893-36f742162938
md"""
### Fonsi's Git Trouble Fix™

Many git problems have an 'official' solution, and it will take time to learn all of these techniques. In the meantime, here is my fool-proof way of solving many tricky git situations:
"""

# ╔═╡ cc1149de-4895-48aa-8335-6dcc78d882c9
md"""
##### Step 1

Use GitHub Desktop to find out which files you have changed. Are the changes important?

##### Step 2

Take any files that you want to keep, and copy them to a different folder, e.g. your desktop.

##### Step 3

**Remove your local clone of the repository, and move any leftover files to the trash.**

##### Step 4

Clone the repository again, and copy your files back to their original place, overwriting the cloned files.
"""

# ╔═╡ a7df39f0-0e97-4fff-9202-cfc629b68f46
bigbreak

# ╔═╡ d1e48204-79b8-4b0f-8dc7-eb69244068de
md"""
# How to fork and contribute code

"""

# ╔═╡ ea12c669-5429-4bf3-af03-378843ca8838
md"""
> Video tutorial will be available here after the live lecture.
"""

# ╔═╡ aeafbdbb-4aeb-452e-a21b-d2c8ec48c64e
md"""
### Tests

People often write _tests_ for their code! As part of the codebase, there will be a `test` folder containing scripts that will import the project and run some basic checks.
"""

# ╔═╡ 40f9fe4d-ddae-4bdb-aee2-7999e288931a
function double(x)
	x + x
end

# ╔═╡ f7425775-55aa-4e46-a11f-7d981a4cfacc
@test double(3) == 6

# ╔═╡ 91b6151d-284e-4f06-954c-ce648fec3327
@test double(0) == 0

# ╔═╡ ddab6132-7d4d-41e8-81f3-ab0fc24cbeeb
md"""
One reason to write tests is to _pin_ specific behaviour, protecting yourself from accidentally changing it later. For example, if you fix a bug that `double(2)` returns `40` instead of `4`, you would also add a test for `double(2) == 4`. Months later, when you are changing `double` for another use case, the old test insures that you are not accidentally breaking something that was once fixed.
"""

# ╔═╡ 7174076d-5eba-4380-8d76-292935014d90
md"""
> ##### Test-driven design
>
> Some people like to _first_ write their tests, which will fail initially, and _then_ write the code to solve the problem. While working on the code, they keep re-running the tests, until all checks are green! This can be an effective and rewarding way to work on a software problem!
>
> If you are following this course, then you will already be familiar with this concept! The homework exercises were all designed with the test-driven principle.
""" |> aside

# ╔═╡ 0d76ea2f-18a9-46d1-8328-f077482d5d1f
md"""
#### Running tests

You can run the tests of any Julia package by opening the Julia REPL, and typing:
```julia
julia> ]

(v1.6) pkg> test Example
```

As mentioned before, many project use _github magic_ to automatically run their tests on a server for every change, and you can view the test results online. An important application is running tests for every Pull Request. This means that package authors can quickly review your changes, knowing that it does not break anything.
"""

# ╔═╡ de95e033-932a-4de9-8e1b-36fcf22c7e20
bigbreak

# ╔═╡ 6efa4ee8-b477-4d47-8d42-5b87f3aa02d2
md"""
# Extra tips
"""

# ╔═╡ f26e39cc-c175-4439-868a-0686250e8e29
md"""

Contributing does not only mean contributing code! [https://opensource.guide/how-to-contribute/](https://opensource.guide/how-to-contribute/)

"""

# ╔═╡ b16a228e-5056-44a0-ab57-0dea5082669d
md"""

Create a "test" repository! Create multiple!

"""

# ╔═╡ 436c9fa2-b770-4dee-82a6-23e9baa551e4
md"""
# Appendix
"""

# ╔═╡ 75c8f825-d988-4f9e-8038-6b4dd2e24181
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="HypertextLiteral", version="0.6"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(url="https://github.com/JuliaPluto/PlutoTest.jl", rev="0c8b9b0"),
    ])
    using HypertextLiteral
	using PlutoTest
	using PlutoUI
end

# ╔═╡ d92d55a3-8fbc-4178-81b4-7ddc379ef7c7
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ 089c9a6c-f4ef-4815-9a99-cb4023b42225
# Layout = ingredients(download("https://raw.githubusercontent.com/fonsp/disorganised-mess/5a59ad7cff1e760b997a54ffa0f8fa202ac16db3/Layout.jl"))

# ╔═╡ 4ea69625-0064-42da-a75a-a54fbd106f78
stackrows(x) = permutedims(hcat(x...),(2,1))

# ╔═╡ 13c0fbf3-08c6-4515-b710-f16b55165a2d
vocabulary(x) = grid(stackrows((
		[@htl("$(k)"), v]
		for (k,v) in x
		)); fill_width=false)

# ╔═╡ 574448b8-2ff1-4bff-8580-33bfcba860e8
function flex(args...; kwargs...)
	Div(;
		contents=collect(args),
		style=Dict("display" => "flex", ("flex-" * String(k) => string(v) for (k,v) in kwargs)...)
		)
end

# ╔═╡ dd4855a0-0b7c-40a5-8565-94b40948f86d
flex(x::Union{AbstractVector,Base.Generator}; kwargs...) = flex(x...; kwargs...)

# ╔═╡ d31f0e84-dce9-4f81-8643-ef08684530d2
begin
	Base.@kwdef struct Div
		contents
		style=Dict()
	end

	Div(x) = Div(contents=x)

	function Base.show(io::IO, m::MIME"text/html", d::Div)
		h = @htl("""

			$(d.contents)

			""")
		show(io, m, h)
	end
end

# ╔═╡ 676ac6ff-1b7e-4c88-b850-45f4375a8d58
function grid(items::AbstractMatrix; fill_width::Bool=true)
	Div(
		contents=Div.(vec(permutedims(items, [2,1]))),
		style=Dict(
			"display" => fill_width ? "grid" : "inline-grid",
			"grid-template-columns" => "repeat($(size(items,2)), auto)",
			"column-gap" => "1em",
		),
	)
end

# ╔═╡ d43abe78-5a9d-4a22-999d-0ee85eb5ab7f
function aside(x)
	@htl("""

		$(x)

		""")
end

# ╔═╡ Cell order:
# ╟─10ebd434-adba-11eb-048f-2d084049d48f
# ╟─ef4eea24-bc1c-43be-b9b6-d073ac2433cf
# ╟─69b8490f-cf04-4e73-bc7b-639b1fc0e2d6
# ╟─cbe5fcba-3ed7-41a6-8932-2693e370c16c
# ╟─812002d3-8603-4ffa-8695-2b2da7f0766a
# ╟─56db20a5-0e8d-4d34-b3ba-c3ab1b8de58e
# ╟─83d0162f-5960-4938-8353-91c4cd220459
# ╟─51448106-3e0e-4abf-84fc-7e6e81425d12
# ╟─3c0059b7-7a99-4c1c-a2bf-e47bdd06a252
# ╟─604837c5-b017-4d6c-a5c5-dab50d5f3f61
# ╟─8678b5e1-0b67-4097-82ba-0daa4e878032
# ╟─0bf51d8c-7adf-4231-8b2d-c976484a3e7c
# ╟─02905480-8864-4e56-af3a-6c7c0789ce6f
# ╟─2ef5db21-0092-4523-b930-0ec99c459ffa
# ╟─51114dc9-cb32-4d31-b780-6f5e372f8763
# ╟─8b24ce23-4ead-4fbc-875e-a8261f671abe
# ╟─3c8abdfc-cf68-4f10-a3e9-08d24803535b
# ╟─dfb316ca-0502-4419-93c2-6d455b7b2f98
# ╟─ae28cad6-2ffe-48b4-895b-fab1bd2f2443
# ╟─86e4b0dd-9bf5-4711-bf3f-b55ed2627a03
# ╟─8ea1ca6b-4bb7-4f53-907d-0e5ca83e5761
# ╟─ee5a3219-4547-4a9d-b527-3489a2925f68
# ╟─e5f2b4e1-f05a-4457-b2e9-091a86ef86bb
# ╟─efbd6d1c-d2c9-48a9-9d6c-dc9ce6af0b5b
# ╟─1421f832-9f97-4ec3-b967-64618983349b
# ╟─b139e8ea-88ad-4df4-9f13-c867edfc2db0
# ╟─e5bab2e1-be9e-4654-844a-d50285e330c8
# ╟─47c67fa6-1490-4134-9e0c-5754cb273d1e
# ╟─79ab206d-2140-4c78-8fd4-a874fe2551e1
# ╟─4ce38ec8-4084-4fe8-a248-87408cd0d39b
# ╟─38261f04-f410-440e-a1cf-218fa240a0ae
# ╟─93ce9618-5484-4572-97f1-1cb9c9367989
# ╟─a0a97cd2-838c-4a2c-9233-969b3274764c
# ╟─0faa4042-42f5-4c74-9270-fbf2205920ca
# ╟─167cbd94-851f-4147-9b9d-9b46c7959785
# ╟─e3fe649e-2733-455a-897d-d4d2d70b9dc6
# ╟─2dd5d359-eb44-4696-bc83-f257eff95618
# ╟─3229f0aa-54d1-43f4-a6aa-cf212620ae13
# ╟─bbb59b45-de20-41c1-b5dd-079fc51655f2
# ╟─7af9e69c-2b81-4a90-861c-ed737a4a9ec4
# ╟─b83ade3d-6f8d-4ac8-9255-956d0a348416
# ╟─4d303641-0299-44d7-ba74-6daea0026b09
# ╟─2c456900-047d-432c-a62a-87f8eeaba8d5
# ╟─dc112303-061e-4e53-8f58-cf9ea0f556f1
# ╟─b2e49cd5-49d5-4ac7-a3ae-9820a97720fb
# ╟─d400d538-4f73-4490-ad68-aedbb57cd70f
# ╟─a98993b9-a5c0-4260-b96e-2655c472ccba
# ╟─67cf205a-3d89-4cd9-ab5e-febc85ea8af1
# ╟─c4442667-072a-4e17-94a4-104a8ec33bd0
# ╟─d8c5273f-7ebc-4399-b893-36f742162938
# ╟─cc1149de-4895-48aa-8335-6dcc78d882c9
# ╟─a7df39f0-0e97-4fff-9202-cfc629b68f46
# ╟─d1e48204-79b8-4b0f-8dc7-eb69244068de
# ╟─ea12c669-5429-4bf3-af03-378843ca8838
# ╟─aeafbdbb-4aeb-452e-a21b-d2c8ec48c64e
# ╠═40f9fe4d-ddae-4bdb-aee2-7999e288931a
# ╠═f7425775-55aa-4e46-a11f-7d981a4cfacc
# ╠═91b6151d-284e-4f06-954c-ce648fec3327
# ╟─ddab6132-7d4d-41e8-81f3-ab0fc24cbeeb
# ╟─7174076d-5eba-4380-8d76-292935014d90
# ╟─0d76ea2f-18a9-46d1-8328-f077482d5d1f
# ╟─de95e033-932a-4de9-8e1b-36fcf22c7e20
# ╟─6efa4ee8-b477-4d47-8d42-5b87f3aa02d2
# ╟─f26e39cc-c175-4439-868a-0686250e8e29
# ╟─b16a228e-5056-44a0-ab57-0dea5082669d
# ╟─436c9fa2-b770-4dee-82a6-23e9baa551e4
# ╠═75c8f825-d988-4f9e-8038-6b4dd2e24181
# ╟─d92d55a3-8fbc-4178-81b4-7ddc379ef7c7
# ╟─089c9a6c-f4ef-4815-9a99-cb4023b42225
# ╟─4ea69625-0064-42da-a75a-a54fbd106f78
# ╟─13c0fbf3-08c6-4515-b710-f16b55165a2d
# ╟─574448b8-2ff1-4bff-8580-33bfcba860e8
# ╟─dd4855a0-0b7c-40a5-8565-94b40948f86d
# ╟─d31f0e84-dce9-4f81-8643-ef08684530d2
# ╟─676ac6ff-1b7e-4c88-b850-45f4375a8d58
# ╟─d43abe78-5a9d-4a22-999d-0ee85eb5ab7f


text:
##
##!------------------------------------------------------------------------------------------------!
##  common.jl
##
##  Common types and functions used by code elsewhere.
##
##  Creation date: 18/08/15
##  Author: Arindam Biswas <arindam dot b at eml dot cc>
##!------------------------------------------------------------------------------------------------!
##

bitstype 128 Reference


text:
@testset "Tensor and tensor_size gives correct size of tensor for the simplest possible system of four spins" begin

    # Model's parameters
    J12 = -1.0
    J13 = -0.5
    J14 = -0.1
    h1 = 0.5
    h2 = 0.75
    h3 = 0.1
    h4 = 0.1

    # dict to be read
    D = Dict(
        (1, 2) => J12,
        (1, 3) => J13,
        (1, 4) => J14,
        (1, 1) => h1,
        (2, 2) => h2,
        (3, 3) => h3,
        (4, 4) => h4,
    )

    # control parameters
    m, n = 2, 2
    t = 1
    β = 1.
    num_states = 4

    # read in pure Ising
    ig = ising_graph(D)

    # construct factor graph with no approx
    fg = factor_graph(
        ig,
        spectrum=full_spectrum,
        cluster_assignment_rule=super_square_lattice((m, n, t)),
    )

    peps = PEPSNetwork{true}(m, n, fg, rotation(0), β=β)
    update_gauges!(peps, :rand)

    @testset "for site tensors" begin
        expected_traced = [(1, 1, 2, 2) (2, 1, 1, 1); (1, 2, 1, 1) (2, 1, 1, 1)]
        TA = tensor_assignment(peps, :site)
        for i ∈ 1:peps.nrows, k ∈ 1:peps.ncols
            A = tensor(peps, (i, k))
            @test size(A) == tensor_size(peps, (i, k)) == expected_traced[i,k]
            @test TA[i, k] == :site
        end
    end

    @testset "for connecting tensors" begin
        expected_connecting = [(2, 1, 2, 2) (1, 1, 1, 1); (1, 2, 2, 1) (1, 1, 1, 1)]
        TA = tensor_assignment(peps, :central_h)
        for i ∈ 1:peps.nrows, k ∈ 1:peps.ncols
            A = tensor(peps, (i, k+1//2))
            @test size(A) == tensor_size(peps, (i, k+1//2)) == expected_connecting[i,k]
            @test TA[i, k+1//2] == :central_h
        end
    end

    @testset "for horizontal and vertical tensors" begin
        expected_vhcentral = [(1, 2, 1, 2) (1, 1, 1, 1)]
        TA = tensor_assignment(peps, :central_v)
        for i ∈ 1:peps.nrows-1, k ∈ 1:peps.ncols
            A = tensor(peps, (i+1//2, k))
            @test size(A) == tensor_size(peps, (i+1//2, k)) == expected_vhcentral[i,k]
            @test TA[i+1//2, k] == :central_v
        end
    end

    @testset "for diagonal tensors" begin
        expected_dcentral = [(1, 2, 1, 2) (1, 1, 1, 1)]
        TA = tensor_assignment(peps, :central_d)
        for i ∈ 1:peps.nrows-1, k ∈ 1:peps.ncols-1
            A = tensor(peps, (i+1//2, k+1//2))
            @test size(A) == tensor_size(peps, (i+1//2, k+1//2)) == expected_dcentral[i,k]
            @test TA[i+1//2, k+1//2] == :central_d
        end
    end

    @testset "for gauges" begin
        expected_gauges = [(1, 2, 1, 2) (1, 2, 1, 2) (1, 1, 1, 1)]
        TA = tensor_assignment(peps, :gauge_h)
        for i ∈ 1:peps.nrows-1, k ∈ 1:1//2:peps.ncols - 1, g ∈ [1//6, 2//6]
            A = tensor(peps, (i+g, k))
            @test size(A) == tensor_size(peps, (i+g, k)) == expected_gauges[Int(i), Int(k)]
            @test TA[i+g, k] == :gauge_h
        end
    end
end

@testset "Tensors in pathological instance have correct sizes" begin
    m = 3
    n = 4
    t = 3

    β = 1.

    schedule = 1.

    L = n * m * t
    states_to_keep = 20

    instance = "$(@__DIR__)/instances/pathological/cross_$(m)_$(n)_dd.txt"

    ig = ising_graph(instance)

    fg = factor_graph(
        ig,
        spectrum=full_spectrum,
        cluster_assignment_rule=super_square_lattice((m, n, t))
    )

    peps = PEPSNetwork{true}(m, n, fg, rotation(0), β=β)
    update_gauges!(peps, :rand)

    expected_site = [
        (1,1,4,2) (4,1,4,2) (8,1,8,2) (4,1,1,2);
        (1,2,8,2) (8,2,8,2) (8,2,4,2) (8,2,1,1);
        (1,2,2,1) (4,2,4,1) (8,2,8,1) (1,1,1,1)
    ]

    expected_connecting = [
        (4,1,4,8) (4,1,8,2) (8,1,4,4) (1,1,1,1);
        (8,8,8,8) (8,2,8,2) (4,4,8,8) (1,1,1,1);
        (2,8,4,1) (4,2,8,1) (8,8,1,1) (1,1,1,1)
    ]

    for i ∈ 1:peps.nrows, k ∈ 1:peps.ncols
        j = denominator(k) == 1 ? numerator(k) : k
        A = tensor(peps, (i, j))
        B = tensor(peps, (i, j+1//2))
        @test size(A) == tensor_size(peps, (i, j)) == expected_site[i,j]
        @test size(B) == tensor_size(peps, (i, j+1//2)) == expected_connecting[i,j]
    end
end


text:
# Simple test of function run_mut_evolution()
# include("../src/CGP.jl") # Assumes that CGP.jl is loaded
# The following lines allow this file to included from CGP.jl, CGP.jl/test, or CGP.jl/src.
#include("../src/Env_evolution.jl")
#include("../src/Assignment.jl")
cwd = pwd()
date = "7_25"
if cwd[end-5:end] == "CGP.jl"
  csvfile = "$cwd/data/$date/test_env_evolveA.csv"
else
  csvfile = "../data/$date/test_env_evolveA.csv"
end
iterations = 1
numinputs = 3
numoutputs = 6
nodearity = 2
numinteriors = 15
ngoals = 1
levelsback=15
#hamming_rng = true:true
maxsteps = 100000:100000
gl_repetitions = 3
num_flip_bits = 2:2
perm_heuristic=true
perturb_goal_range = false:true
avgfitness = true
fit_limit_list = [numoutputs-2.0, numoutputs-1.0, Float64(numoutputs)]
#run_env_evolution( iterations, numinputs, numoutputs, numinteriors, goallistlength, maxsteps, levelsback, csvfile )
df = run_env_evolution( iterations, numinputs, numoutputs, numinteriors, ngoals, maxsteps, levelsback,
    gl_repetitions, num_flip_bits, avgfitness, perm_heuristic, perturb_goal_range, fit_limit_list, csvfile )


text:
"""
    MatrixHamiltonian(
        mat::AbstractMatrix{T};
        starting_address::Int = starting_address(mat)
    ) <: AbstractHamiltonian{T}
Wrap an abstract matrix `mat` as an [`AbstractHamiltonian`](@ref) object for use with
regular `Vector`s indexed by integers. Works with stochastic methods of
[`lomc!()`](@ref). Optionally, a [`starting_address`](@ref) can be provided.

Specialised methods are implemented for sparse matrices of type `AbstractSparseMatrixCSC`.
"""
struct MatrixHamiltonian{T,AM,H} <: AbstractHamiltonian{T}
    m::AM
    starting_index::Int
end

function MatrixHamiltonian(
    m::AM;
    starting_address=starting_address(m)
) where AM <:AbstractMatrix
    s = length.(axes(m))
    @assert s[1]==s[2] "Matrix needs to be square, got $s."
    i = Int(starting_address)
    @assert minimum(axes(m, 2)) ≤ i ≤ maximum(axes(m, 2)) "Invalid index $starting_address."
    MatrixHamiltonian{eltype(m), AM, ishermitian(m)}(m, i)
end

LOStructure(::Type{<:MatrixHamiltonian{<:Any,<:Any,true}}) = IsHermitian()
LOStructure(::Type{<:MatrixHamiltonian}) = AdjointKnown()
LinearAlgebra.adjoint(mh::MatrixHamiltonian) = MatrixHamiltonian(mh.m')
LinearAlgebra.adjoint(mh::MatrixHamiltonian{<:Any,<:Any,true}) = mh

starting_address(mh::MatrixHamiltonian) = mh.starting_index

dimension(::Type{T}, mh::MatrixHamiltonian) where T = T(size(mh.m,2))

num_offdiagonals(mh::MatrixHamiltonian, _) = dimension(mh) - 1

function get_offdiagonal(mh::MatrixHamiltonian, add, chosen)
    newadd = ifelse(chosen < add, chosen, chosen+1)
    return newadd, mh.m[newadd, add]
end
diagonal_element(mh::MatrixHamiltonian, i) = mh.m[i,i]

# specialised methods for sparse matrices - avoid spawning via zero value matrix elements
function num_offdiagonals(
    mh::MatrixHamiltonian{<:Any, <:SparseArrays.AbstractSparseMatrixCSC},
    j
)
    nnz(mh.m[:, j]) - 1
end

function offdiagonals(
    mh::MatrixHamiltonian{T, <:SparseArrays.AbstractSparseMatrixCSC},
    col::Integer
) where T
    rows = rowvals(mh.m)[nzrange(mh.m, col)]
    vals = nonzeros(mh.m)[nzrange(mh.m, col)]
    drow = findprev(x->x==col, rows, min(col,length(rows))) # rows[drow]==col should be true
    return SparseMatrixOffdiagonals{Int,T,typeof(rows),typeof(vals)}(rows, vals, drow, col)
end

struct SparseMatrixOffdiagonals{A,T,R,V} <: AbstractOffdiagonals{A,T}
    rows::R # indices of rows with nonzero values
    vals::V # nonzero values
    drow::Int # index of col in rows; rows[drow] == col should be true
    col::Int # colum of the matrix
end
function Base.getindex(smo::SparseMatrixOffdiagonals, chosen)
    ind = ifelse(chosen < smo.drow, chosen, chosen + 1)
    return smo.rows[ind], smo.vals[ind]
end
Base.size(smo::SparseMatrixOffdiagonals) = (length(smo.rows) - 1,)

function get_offdiagonal(
    mh::MatrixHamiltonian{<:Any, <:SparseArrays.AbstractSparseMatrixCSC},
    add,
    chosen
)
    return offdiagonals(mh, add)[chosen]
end


text:
"""
    'ODC!(pokes::AbstractDataFrame)'

Given the pokes Dataframe adds 2 columns
ODC: the ratio between each nose-poke’s duration and the sum of the duration and
the preceding inter-poke interval
PokeFromLeaving: An index about how many pokes are left in a given trial before
the animal left
"""

function ODC(pokes::AbstractDataFrame)
    df = copy(pokes)
    pn = propertynames(df)
    if !(:Pre_Interpoke in pn) || !(:PokeDuration in pn)
        error("Error in odc calculation, missing PreInterpoke or PokeDuration info")
    end
    df[!,:ODC] = [d/(d+i) for (d,i) in zip(df[:,:PokeDuration],df[:,:Pre_Interpoke])]
    gd = groupby(df,[:Day,:MouseID,:Trial])
    transform!(gd,:Poke => p -> collect(1:length(p)))
    rename!(df, :Poke_function => :PokeFromTrial)
    transform!(gd,:Poke => p -> collect(length(p)-1:-1:0))
    rename!(df, :Poke_function => :PokeFromLeaving)
    transform!(gd,:Reward => r -> count_from_last(r))
    rename!(df, :Reward_function => :PokeFromLastRew)
    return df
end

function count_from_last(v::AbstractVector)
    l = length(v)
    if l == 1
        return [missing]
    end
    last_r  = findlast(v)
    if isnothing(last_r)
        return missings(l)
    else
        return collect(1:l) .- last_r
    end
end

function calculate_odc(data::AbstractDataFrame,Phase::AbstractString,Allignment::Symbol)
    df0 = dropmissing(data,Allignment)
    # filter!(Allignment => a -> -limit < a < limit, df0)

    df1 =  filter(r -> !ismissing(r.ODC) &&
        !r.Reward &&
        r.Phase in [Phase],
        df0)

    df2 = combine(groupby(df1,:Treatment)) do dd
        summarize(dd,Allignment,:ODC)
    end
end

function calculate_odc(data::AbstractDataFrame,Allignment::Symbol; normalize = nothing)
    if !in(:ODC,propertynames(data))
        error("ODC column missing from data")
    end
    df0 = dropmissing(data,Allignment)
    df1 =  filter(r -> !ismissing(r.ODC) &&
        !r.Reward,
        df0)
    gd = groupby(df1,[:Phase,:Treatment,:MouseID, Allignment])
    df2 = combine(:ODC => mean,gd)
end


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

const CPU_PHYSICAL_CORES = num_physical_cores()

function get_imfilter_impl(GPU)
  if GPU ≠ nothing
    imfilter_gpu
  else
    imfilter_cpu
  end
end

function convdist(Xs::AbstractArray, masks::AbstractArray; weights=nothing)
  # choose among imfilter implementations
  imfilter_impl = get_imfilter_impl(GPU)

  # default to uniform weights
  weights == nothing && (weights = ones(masks[1]))

  result = []
  for (X, mask) in zip(Xs, masks)
    wmask = weights.*mask

    A² = imfilter_impl(X.^2, weights)
    AB = imfilter_impl(X, wmask)
    B² = sum(abs2, wmask)

    push!(result, abs.(A² - 2AB + B²))
  end

  D = sum(result)

  # always return a plain simple array
  parent(D)
end

function genpath(extent::NTuple{3,Integer}, kind::Symbol, datum=[])
  path = Int[]

  if kind == :rasterup
    for k=1:extent[3], j=1:extent[2], i=1:extent[1]
      push!(path, sub2ind(extent, i,j,k))
    end
  end

  if kind == :rasterdown
    for k=extent[3]:-1:1, j=1:extent[2], i=1:extent[1]
      push!(path, sub2ind(extent, i,j,k))
    end
  end

  if kind == :random
    nelm = prod(extent)
    path = nthperm!(collect(1:nelm), rand(1:factorial(big(nelm))))
  end

  if kind == :dilation
    nelm = prod(extent)
    pivot = rand(1:nelm)

    grid = falses(extent)
    grid[pivot] = true
    push!(path, pivot)

    while !all(grid)
      dilated = dilate(grid, [1,2,3])
      append!(path, find(dilated - grid))
      grid = dilated
    end
  end

  if kind == :datum
    @assert !isempty(datum) "datum path cannot be generated without data"

    shuffle!(datum)

    grid = falses(extent)
    for (i,j,k) in datum
      pivot = sub2ind(extent, i,j,k)
      grid[pivot] = true
      push!(path, pivot)
    end

    while !all(grid)
      dilated = dilate(grid, [1,2,3])
      append!(path, find(dilated - grid))
      grid = dilated
    end
  end

  path
end


text:
module Generator

using Pkg
using Genie, Genie.Generator

function newapp(app_name::String; ui::Bool = true, charts::Bool = true, autostart::Bool = true, kwargs...)
  Genie.Generator.newapp(app_name; autostart = false, kwargs...)

  ui && Pkg.add("StippleUI")
  charts && Pkg.add("StippleCharts")

  autostart && Genie.up()
end

function scaffold(app_name::String)
  clean_name = Genie.Generator.validname(app_name)
  type_name = "$(uppercasefirst(clean_name))Model"
  var_name = "$(lowercase(clean_name))_model"

  Pkg.activate(".")
  Pkg.add("Revise")
  Pkg.add("Stipple")
  Pkg.add("StippleUI")
  Pkg.add("StippleCharts")

  open("$clean_name.jl", "w") do io
    write(io,
    """
    __revise_mode__ = :eval
    using Revise
    using Stipple
    using StippleUI
    using StippleCharts

    using GenieAutoReload
    GenieAutoReload.autoreload(pwd())

    #= Data =#

    Base.@kwdef mutable struct $type_name <: ReactiveModel
      # add fields here
      x::R{Int} = 100
      y::R{String} = "Hello"
    end

    #= Stipple setup =#

    Stipple.register_components($type_name, StippleCharts.COMPONENTS)
    $var_name = Stipple.init($type_name())

    #= Event handlers =#

    onany($var_name.x, $var_name.y) do (_...)
      # handle update events
      @show $var_name.x[] $var_name.y[]
    end

    #= UI =#

    function ui(model::$type_name)
    [
      dashboard(
        vm(model), class="container", title="$app_name", head_content=Genie.Assets.favicon_support(),
        [
          heading("$app_name")

          row([
            cell(class="st-module", [
              h6("X")
              span("", @text(:x))
            ])
            cell(class="st-module", [
              h6("Y")
              span("", @text(:y))
            ])
          ])

          footer(class="st-footer q-pa-md", [
            cell([
              span("$app_name")
            ])
          ])
        ]
      )
      GenieAutoReload.assets()
    ]
    end

    #= routing =#

    route("/") do
      ui($var_name) |> html
    end

    #= start server =#
    # up(rand((8000:9000)), open_browser=true)
    """)
  end
end

end


text:
"""
$(TYPEDEF)

Coordinates as row (vertical), columns (horizontal). Represent the data of each top corner of a pixel.
Topleft for top triangle, bottomleft for bottom triangle
vr, vc: velocity vector of that corner.
Warning: posts/fences: One more corner than number of pixels.
"""
struct Vertex3D
    r::Float64
    c::Float64
    ϕ::Float64

    vr::Float64
    vc::Float64

    Vertex3D(r, c, ϕ, vr, vc) = new(r, c, ϕ, vr, vc)
    Vertex3D() = new(0.0, 0.0, 0.0, 0.0, 0.0)
end

"""
$(SIGNATURES)
"""
function dist(p1::Tuple{Float64,Float64,Float64}, p2::Tuple{Float64,Float64,Float64})
    dr = p2[1] - p1[1]
    dc = p2[2] - p1[2]
    return sqrt(dr^2 + dc^2)
end

"""
$(SIGNATURES)
"""
function dist(p1::Vertex3D, p2::Vertex3D)
    dr = p2.r - p1.r
    dc = p2.c - p1.c
    return sqrt(dr^2 + dc^2)
end

"""
$(SIGNATURES)
"""
function area(v1::Vertex3D, v2::Vertex3D, v3::Vertex3D)
    a = dist(v1, v2)
    b = dist(v2, v3)
    c = dist(v3, v1)
    s = (a + b + c) / 2.0

    # surface_sq = (a + b + c) * (- a + b + c) * (a - b + c) * (a + b - c) / 16.0
    surface_sq = s * (s - a) * (s - b) * (s - c)
    return surface_sq <= 1e-100 ? 0.0 : sqrt(surface_sq)
end

"""
$(SIGNATURES)

Centroid of three points.
"""
centroid(p1::Vertex3D, p2::Vertex3D, p3::Vertex3D) =
    Vertex3D((p1.r + p2.r + p3.r) / 3.0, (p1.c + p2.c + p3.c) / 3.0, 0.0, 0.0, 0.0)

"""
$(SIGNATURES)
"""
function fill_borders!(mat::AbstractMatrix{T}, val::T) where {T}
    mat[begin:end, begin] .= val
    mat[begin:end, end] .= val
    mat[begin, begin:end] .= val
    mat[end, begin:end] .= val
end

"""
$(TYPEDEF)

## Coordinates

Origin is top left, going right and down .`x` goes horizontal and follows columns. x is within 1 and width
`y` goes vertical and follows rows. y is within 1 and height.

# Velocity

Velocity vector along `row` and `col` (velocity along h not used).

Implementation is a struc of arrays. Easier to vectorise.

"""
mutable struct FieldVertex3D
    size::Tuple{Int,Int}

    r::AbstractMatrix{Float64}
    c::AbstractMatrix{Float64}
    vr::AbstractMatrix{Float64}
    vc::AbstractMatrix{Float64}

    # Velocity potential at each corner
    ϕ_top::AbstractMatrix{Float64}
    ϕ_bot::AbstractMatrix{Float64}

    rows_numbers::AbstractMatrix{Float64}
    cols_numbers::AbstractMatrix{Float64}
end

# Information about each corner surrounding a pixel. Posts&fences warning!!!
function FieldVertex3D(height, width)
    # Two matrices comtaining the row and column number of each element of the matrix
    rows_numbers = repeat(Float64.(1:height+1), 1, width + 1)
    cols_numbers = repeat(Float64.(1:width+1)', height + 1, 1)

    mr = rows_numbers
    mc = cols_numbers
    mvr = zeros(Float64, height + 1, width + 1)
    mvc = zeros(Float64, height + 1, width + 1)

    mϕ_top = zeros(Float64, height , width )
    mϕ_bot = zeros(Float64, height , width )

    fv = FieldVertex3D((height, width), mr, mc, mvr, mvc, mϕ_top, mϕ_bot, rows_numbers, cols_numbers)
    reset_border_values!(fv)
    return fv
end

"""
$(SIGNATURES)

Reset the values of the borders corners: all positions at integer values, all velocities at 0.
"""
function reset_border_values!(corners::FieldVertex3D)
    corners.r[1, :] .= corners.rows_numbers[1, :]
    corners.r[end, :] .= corners.rows_numbers[end, :]

    corners.r[:, 1] .= corners.rows_numbers[:, 1]
    corners.r[:, end] .= corners.rows_numbers[:, end]

    corners.c[1, :] .= corners.cols_numbers[1, :]
    corners.c[end, :] .= corners.cols_numbers[end, :]

    corners.c[:, 1] .= corners.cols_numbers[:, 1]
    corners.c[:, end] .= corners.cols_numbers[:, end]

    fill_borders!(corners.vr, 0.0)
    fill_borders!(corners.vc, 0.0)
end

Base.size(fv::FieldVertex3D) = fv.size

"""
$(SIGNATURES)
"""
function Vertex3D(fv::FieldVertex3D, row, col)
    height, width = size(fv)

    if (1 <= row <= height ) && (1 <= col <= width)
        Vertex3D(
            fv.r[row, col],
            fv.c[row, col],
            fv.vr[row, col],
            fv.vc[row, col],
            fv.ϕ_top[min(height, row), col(width, col)],
            fv.ϕ_bot[min(height, row), col(width, col)],
        )
    else
        missing
    end
end

"""
$(TYPEDEF)

The algorithm works by allocating a rectangle from the exit mesh (the face where the light exits the block
of acrylate). Initially, each rectangle representsa single pixel on that face and faces s single pixel on
the caustics picture.

Each rectangle is split into 2 trangles: one upper-left, the other bottom-right like:

```
*------*
|    / |
|   /  |
|  /   |
| /    |
*------*
```
Although the corners of each triangle match corners of the rectangles, it is sometimes easier to think in terms
of triangles rather than corners coming from different rectangles.
"""
struct FaceMesh
    size::Tuple{Int,Int}
    corners::FieldVertex3D

    toptriangles::AbstractMatrix{Tuple{Tuple{Int,Int},Tuple{Int,Int},Tuple{Int,Int}}}
    bottriangles::AbstractMatrix{Tuple{Tuple{Int,Int},Tuple{Int,Int},Tuple{Int,Int}}}

    FaceMesh(s, tl, tt, bt) = new(s, tl, tt, bt)

    function FaceMesh(height::Int, width::Int)
        m_corner = FieldVertex3D(height, width)
        t_triangles = Matrix{Tuple{Tuple{Int,Int},Tuple{Int,Int},Tuple{Int,Int}}}(
            undef,
            height,
            width,
        )
        b_triangles = Matrix{Tuple{Tuple{Int,Int},Tuple{Int,Int},Tuple{Int,Int}}}(
            undef,
            height,
            width,
        )

        for row ∈ 1:height, col ∈ 1:width
            t_triangles[row, col] = ((row, col), (row + 1, col), (row, col + 1))
            b_triangles[row, col] = ((row, col + 1), (row + 1, col), (row + 1, col + 1))
        end

        return FaceMesh((height, width), m_corner, t_triangles, b_triangles)
    end
end

"""
$(SIGNATURES)

Return the dimension of the mesh as the number of rectangles. It does not return the number of corner points.
"""
Base.size(mesh::FaceMesh) = mesh.size

"""
$(SIGNATURES)

Converts a triangle as a triplet of references to mesh vertices to a triplet of 3D coordinates.
"""
function triangle3D(mesh::FaceMesh, row::Int, col::Int, side = Union{:top,:bottom})
    height, width = size(mesh)
    max_distance_squared = (height + 1)^2 + (width + 1)^2

    if 1 <= row <= height && 1 <= col <= width
        if side == :top
            t = mesh.toptriangles[row, col]
        elseif side == :bottom
            t = mesh.bottriangles[row, col]
        end

        t1_row, t1_col = t[1]
        p1 = Vertex3D(mesh.corners, t1_row, t1_col)

        t2_row, t2_col = t[2]
        p2 = Vertex3D(mesh.corners, t2_row, t2_col)

        t3_row, t3_col = t[3]
        p3 = Vertex3D(mesh.corners, t3_row, t3_col)

        return (p1, p2, p3)
    else
        # return (missing, missing, missing)
        @assert false "Coordinates out of bounds at $(row), $(col)"
    end
end

triangle3D(mesh::FaceMesh, ci::CartesianIndex{2}, side = Union{:top,:bottom}) =
    triangle3D(mesh, ci[1], ci[2], side)

"""
$(SIGNATURES)

Centroid of a specific mesh triangle.
"""
function centroid(mesh::FaceMesh, ci::CartesianIndex; side = Union{:top,:bottom})
    if side == :top
        return centroid(top_triangle3D(mesh, ci)...)
    elseif side == :bottom
        return centroid(bot_triangle3D(mesh, ci)...)
    else
        return missing
    end
end

"""
$(SIGNATURES)

Centroid of a specific mesh triangle.
"""
function centroid(mesh::FaceMesh, row::Int, col::Int; side = Union{:top,:bottom})
    if side == :top
        return centroid(top_triangle3D(mesh, row, col)...)
    elseif side == :bottom
        return centroid(bot_triangle3D(mesh, row, col)...)
    else
        return missing
    end
end

"""
$(SIGNATURES)
"""
function area(mesh::FaceMesh, ci::CartesianIndex; side = Union{:top,:bottom})
    if side == :top
        return area(top_triangle3D(mesh, ci)...)
    elseif side == :bottom
        return area(top_triangle3D(mesh, ci)...)
    else
        return missing
    end
end

"""
$(SIGNATURES)
"""
function area(mesh::FaceMesh, row::Int, col::Int; side = Union{:top,:bottom})
    if side == :top
        return area(top_triangle3D(mesh, row, col)...)
    elseif side == :bottom
        return area(top_triangle3D(mesh, row, col)...)
    else
        return missing
    end
end

"""
$(SIGNATURES)
"""
area(t::Tuple{Vertex3D,Vertex3D,Vertex3D}) = area(t...)

"""
$(SIGNATURES)

A Mesh is a collection of triangles. The brightness flowing through a given triangle is just proportional to its
area in the r, r plane. h is ignored.

The function returns a matrix with the quantity of light coming from each 'rectangle'  around a corner. That 'rectangle'
has been shifted and flexed around.
"""
function get_lens_pixels_area(mesh::FaceMesh)
    height, width = size(mesh)

    top_tri_area = get_lens_pixels_area(mesh, :top)
    bot_tri_area =get_lens_pixels_area(mesh, :bottom)

    total_area = top_tri_area + bot_tri_area

    average_energy_per_pixel = average(total_area)
    luminosity_pixels = total_area / average_energy_per_pixel

    @assert abs(sum(luminosity_pixels) - height * width) < 1.0 "Total lens luminosity is incorrect"

    # return luminosity_pixels
    return top_tri_area, bot_tri_area
end

function get_lens_pixels_area(mesh::FaceMesh, side::Union{:top, :bottom})
    height, width = size(mesh)

    # return luminosity_pixels
    return area.([triangle3D(mesh, row, col, side) for row ∈ 1:height, col ∈ 1:width])
end

"""
$(SIGNATURES)

"""
function field_summary(field, fieldname)
    s = round(sum(field), sigdigits = 6)
    max = round(maximum(field), sigdigits = 4)
    min = round(minimum(field), sigdigits = 4)
    avg = round(average(field), sigdigits = 4)
    abs = round(average_absolute(field), sigdigits = 4)

    return """$(fieldname) (Sum/Min/Max/Avg/Avg abs.): $(s) / $(min) / $(max) / $(avg) / $(abs)"""
end

"""
$(SIGNATURES)

"""
field_summary(field) = field_summary(field, "")


text:
using SIMDTypes
using Test

@testset "SIMDTypes.jl" begin
  @test SIMDTypes.Bit(true).data
  @test !SIMDTypes.Bit(false).data
  @test SIMDTypes.Bit(true) isa SIMDTypes.NativeTypes
  @test !(SIMDTypes.Bit(true) isa SIMDTypes.NativeTypesExceptBit)
  @test SIMDTypes.Bit(true) isa SIMDTypes.NativeTypesExceptFloat16
  @test !(SIMDTypes.Bit(true) isa SIMDTypes.NativeTypesExceptBitandFloat16)
  @test Float16(1) isa SIMDTypes.NativeTypes
  @test Float16(1) isa SIMDTypes.NativeTypesExceptBit
  @test !(Float16(1) isa SIMDTypes.NativeTypesExceptFloat16)
  @test !(Float16(1) isa SIMDTypes.NativeTypesExceptBitandFloat16)
  @test 1f0 isa SIMDTypes.NativeTypesExceptBitandFloat16

  @test 1.3 isa SIMDTypes.FloatingTypes
  # @test !(1.3 isa SIMDTypes.IntegerTypes)
  @test !(1.3 isa SIMDTypes.IntegerTypesHW)
  # @test 1 isa SIMDTypes.IntegerTypes
  @test 1 isa SIMDTypes.SignedHW
  @test !(1 isa SIMDTypes.UnsignedHW)
  @test 1 isa SIMDTypes.IntegerTypesHW
  # @test one(UInt) isa SIMDTypes.IntegerTypes
  @test !(one(UInt) isa SIMDTypes.SignedHW)
  @test one(UInt) isa SIMDTypes.UnsignedHW
  @test one(UInt) isa SIMDTypes.IntegerTypesHW

  @test ntuple(Core.VecElement, 2) isa SIMDTypes._Vec
end


text:
export GuildEmbed

"""
A [`Guild`](@ref) embed.
More details [here](https://discordapp.com/developers/docs/resources/guild#guild-embed-object).
"""
struct GuildEmbed
    enabled::Bool
    channel_id::Union{Snowflake, Nothing}
end
@boilerplate GuildEmbed :dict :docs :lower :merge