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

using Markdown
using InteractiveUtils

# ╔═╡ 7dbb1cfc-46c0-11eb-3cde-2dacfdb76910
using Pkg

# ╔═╡ b2bdd276-46c0-11eb-04ea-077ed82cf318
using SpecialFunctions;

# ╔═╡ 98ad6aee-46c3-11eb-2659-1df8839b5675
begin
	module MyModule
	export func2
	a=42
	function func1(x)
	    return x^2
	end

	function func2(x)
	    return func1(x) + a
	end

	end #end of module
	using .MyModule
end

# ╔═╡ 426b72dc-46c0-11eb-08ea-490eb2b4e2a2
md" ## _Modules_
In this lesson we will learn what modules are and how they can be used for code reusability.

## Working with Modules
Libraries in Julia come in the form of module which can be loaded via the using notation. A module is a separate environment with its sets of variables and functions, some of which are exported in the calling scope, which means that you can call exported functions by simply typing their name as if they where defined in the same scope, while others are accessible only through the ModuleName.functionName notation.

In order to use an existing official module, we need first to install it and then import it, you can do it using the package manager. For this example we will use the Special Functions package, which contains functions such as the gamma function and the Bessel functions.
"

# ╔═╡ 84828a4a-46c0-11eb-23a6-3708b6c8629f
Pkg.add("SpecialFunctions")

# ╔═╡ aeff4b1c-46c0-11eb-2c22-dfafe3c2a2c9
md"At line 1 we import the module called Pkg (with the using keyword) and at line 2 we call the add function which is defined inside it. add takes as its argument the name of the package which we want to install and it will download and build it for us. When it is done (it may take a few minutes) we are ready to use the functions available inside the package!"

# ╔═╡ c15fe60c-46c0-11eb-25e3-3703cca3e55b
sinint(5)

# ╔═╡ c7e491aa-46c0-11eb-04e1-9d578a50a92e
md"If we don’t want to import all of the functions available inside SpecialFunctions but only some of them, for example the gamma function and sinint, but not cosint, we can do it in the following way."

# ╔═╡ 0fadc6aa-46c1-11eb-35dd-153ff4ca0078
cosint(5)

# ╔═╡ 17321bc4-46c1-11eb-0f87-c5e661fe114e
md"Sometimes it is useful to call a function taking in consideration the package where it is defined. If a function is not exported by the module (more on export later) or if there are several modules which export a function with the same name and same argument signature, we can specify which module the function belongs to using the following syntax:"

# ╔═╡ ba1df408-46c0-11eb-3db7-47c4e4e54059
gamma(3)

# ╔═╡ 07fb1dae-46c1-11eb-2c8f-a1a393fa332c
gamma(4)

# ╔═╡ 5621266c-46c2-11eb-0b1e-b1a4906e36f6
gamma(3)

# ╔═╡ 6df7227a-46c2-11eb-05a1-5dbaab866f31
SpecialFunctions.gamma(3)

# ╔═╡ 82434edc-46c2-11eb-3bf0-d56282316427
md"As you can see at line 7, a warning is shown to let us know that the gamma function exported by SpecialFunctions is in conflict with the gamma function which we have defined at line 1.

At line 9 we call the gamma function and as we can see the first definition of the function is what is used (i.e. the user defined function), if we want to call the gamma function inside SpecialFunctions we need to specify the module which contains it (as done on line 13).

In the case where a conflict like this may arise, it is better to avoid the using notations and instead use import, which “imports” the desired module all the same but without exporting any function in the calling scope."

# ╔═╡ a38c3be4-46c2-11eb-2555-8b66580eb7fc
code

# ╔═╡ a8458c44-46c2-11eb-02be-bd3aeb05d50a
md" At line 4 we see that gamma is not defined if we simply type import SpecialFunctions and at line 6 we see that it becomes mandatory to call it using SpecialFunctions.gamma(3), after we define the gamma function at line 9-12 it becomes possible to call our “custom” gamma function.

## User Defined Modules

It is possible to think of modules as compact blocks of variables and functions which can be easily imported in another program. One should not think of a module as something similar to what is a class in object oriented languages such as C++ and Python, but instead as a separate global scope with its own set of variables and functions which can be called from another program. One difference with a class is that it is not possible to import a module several times to have different sets of “global variables”, while it is usual in OOP languages to have different instances of the same class.

In Julia the functions inside a module should thus depend only on their input (and eventually some global variables or other functions defined inside the module which will be shared). We will see in the next lesson how all the data which should be passed to a function and which has to be eventually modified by the function can be conveniently wrapped into a structure. As an anticipation, what may look like object.foo(x) in Python/C++, in Julia will look like foo(x, dataStructure). Although it may seem a trivial difference, we will see that this difference is at the base of multiple dispatch, which is one of the main strength of Julia!

In the following example we will define a simple module which performs some operations and exports a function .
"

# ╔═╡ d77a0ed8-485a-11eb-0dac-3ff8dee5168a
md"
A module starts with the module keyword and should end with end. Contrarily to functions and other blocks, one should not add indentation to a module block.

We define a variable, a, and two functions, func1 and func2, but we export only func2 (see line 2), which means that only func2 will be accessible in the outer scope if we don’t specify the module to which it belongs.

At line 17 we import MyModule, notice the . before the module name, this is needed as MyModule is not an “official” package and ultimately because MyModule is defined in the Main scope (the .Name notation is an abbreviation of Main.Name). At line 19 we call func2 and at line 22 func1, notice that an error is thrown when we call func1 as it is not exported. If we want to access it we need to type MyModule.func1(3) as shown at line 25.

## Code inclusion

As modules become bigger and longer, it is a good practice to split code among different files. In Julia it is simple to include code from another file in the current program as if the code were written inside the same file through the include() function.

include() takes as its argument the name of the file which should be “pasted” inside the current program, pay attention to its location as you may need to write the path to the file.

For this example let’s create three files as shown in the picture:
"

# ╔═╡ fa5636fc-485a-11eb-2f34-75e68933eb82

# ╔═╡ fdbbaca2-485a-11eb-2705-9370aa177831
md"The content of `big-module-part-1.jl` should be:"

# ╔═╡ 13e98006-485b-11eb-1d57-ed09c6bde282
function func1big(x)
    return x^2
end

# ╔═╡ 252120ce-485b-11eb-3637-c14cc1d0153c
md"The content of big-module-part-2.jl should be:"

# ╔═╡ 2bace11c-485b-11eb-01a3-630898d751eb
begin
	a = 42

	function func2big(x)
	    return func1big(x) + a
	end
end

# ╔═╡ 35c07c8e-485b-11eb-24f2-9b5b8cececd6
md"And the content of big-module.jl should be:"

# ╔═╡ 52ee9660-485b-11eb-1e3f-9be2725ee183
md"We can now import MyBigModule inside our program (lesson-modules.jl) in the following way:"

# ╔═╡ 92ea00b0-485b-11eb-3db8-7d6ec720cfcb
md"Although one should give meaningful names to the files which make up a module (and not part1, part2, etc.) this was an example of how one can structure a module:

make a main “module file” which contains the module, imports all of the other files using include and exports the desired functions.
make several files with meaningful names which perform a group of operations with a common topic.
This structure lets you easily extend the module (simply add new files) and makes the code more maintainable, if functions which perform similar tasks are grouped in the same file.

Notice that in Julia it is not important the order in which you include the files in the main module, but in my opinion it is a good practice to include the files in a sort of chronological order. In our example the function func2big depends on func1big which is defined inside big-module-part-1.jl, so we should import it before we includebig-module-part-2.jl.

Remember that when we call include() the code gets “pasted” inside the file where include is called, so it is not necessary to call include inside big-module-part-2.jl as the compiler will see part1, part2 and the main module as a unique file (furthermore, using include inside a file which gets included by another file may lead to errors).

## Code reusability

The goal of a module is to write a set of functions, define a series of variables or types which can be easily reused in other programs (your other programs for example), so one should give meaningful names to the functions (not like func1 and func2) and make them as general and stand-alone as possible.

Since you will likely come back to a module you have written once every few months, you may not remember what a specific function does: it is thus a good idea to add annotations to your code through # comment and write a small description of the function.

When you type ? functionName in the REPL you get a description of that function and usually an example of how it can be used. We will now learn how to write such description for our functions.

## Code documentation
One can write the description of a function in the following way:
"

# ╔═╡ 2323afe6-485c-11eb-18db-c7f6864abb07
"""
Description of the function
"""
function foo(x)
    #... function implementation
end

# ╔═╡ 33fd5894-485c-11eb-281d-c1fe27b864c1
md"There is a set of rules for writing documentation listed in the official documentation. Here is a summary:

Start with the function signature (i.e. function name and arguments) indented by 4 spaces
Write a small summary of what the function does
If necessary explain what the arguments mean/do
Optionally give an example of usage of the function
It is advisable to prepend `@doc raw...`, in this way you will be able to write markdown code inside the description string without the need to escape special characters.

Let’s see how we can document, for example, func2big:"

# ╔═╡ 6a165f92-485d-11eb-3c7a-8b6501eccc03
@doc raw"""
    func4big(x)

Computes the square of `x` and add `a` to it.

# Examples
```julia-repl
julia> include("big-module.jl")
julia> using .MyBigModule
julia> func4big(3)
51
```
"""
function func4big(x)
    return func1big(x) + a
end

# ╔═╡ c9ba62ea-485d-11eb-2e1c-a58ff592ca6c
md"After we have imported MyBigModule we can type `? func2big` in the REPL to read the documentation:"

# ╔═╡ f2741294-485d-11eb-15db-29302d1a7f4c

# ╔═╡ 2a69538a-485e-11eb-2d4b-4b29a9bdd635
md"Or, if you are using the Juno IDE, we can look up the documentation for func2big in the documentation explorer (the documentation explorer may take a few minutes to compile the documentation the first time you use it):"

# ╔═╡ 355e34d6-485e-11eb-3a15-939e18fc8978
md"Pretty neat, don’t you think?

It is good practice to write the documentation with at least the function signature and a short description for each function you define: this will make navigating through the code much easier and you can keep open the documentation tab in Juno to look for the signature (i.e. the arguments) of a specific function you wrote.

## Conclusions

We have learned how to import an existing “official” module and how to write our own. We have learned how it is possible to split a piece of code between multiple files and how code reusability can be improved by module usage. Finally, we have learned how to write proper code documentation in order to make it easier to find out and remember what a function does."

# ╔═╡ 221b0d82-46c1-11eb-2605-29fbe2572d2f
begin
	function gamma(x)
		println("I am another 'gamma' function")
		return x^2
	end
	gamma(3)
end

# ╔═╡ d6be33ca-46c0-11eb-3d75-7312125cded8
using SpecialFunctions: gamma, sinint

# ╔═╡ 837a9644-485b-11eb-23af-c7439f1b6ac0
begin
	include("big-module.jl")

	using .MyBigModule
end

# ╔═╡ 3e922ce8-485b-11eb-2a34-57c74d16a76e
module MyBigModule
include("files/big-module-part-1.jl")
include("files/big-module-part-2.jl")

export func2big

end #end of module

# ╔═╡ Cell order:
# ╟─426b72dc-46c0-11eb-08ea-490eb2b4e2a2
# ╠═7dbb1cfc-46c0-11eb-3cde-2dacfdb76910
# ╠═84828a4a-46c0-11eb-23a6-3708b6c8629f
# ╟─aeff4b1c-46c0-11eb-2c22-dfafe3c2a2c9
# ╠═b2bdd276-46c0-11eb-04ea-077ed82cf318
# ╠═ba1df408-46c0-11eb-3db7-47c4e4e54059
# ╠═c15fe60c-46c0-11eb-25e3-3703cca3e55b
# ╟─c7e491aa-46c0-11eb-04e1-9d578a50a92e
# ╠═d6be33ca-46c0-11eb-3d75-7312125cded8
# ╠═07fb1dae-46c1-11eb-2c8f-a1a393fa332c
# ╠═0fadc6aa-46c1-11eb-35dd-153ff4ca0078
# ╟─17321bc4-46c1-11eb-0f87-c5e661fe114e
# ╠═221b0d82-46c1-11eb-2605-29fbe2572d2f
# ╠═5621266c-46c2-11eb-0b1e-b1a4906e36f6
# ╠═6df7227a-46c2-11eb-05a1-5dbaab866f31
# ╟─82434edc-46c2-11eb-3bf0-d56282316427
# ╠═a38c3be4-46c2-11eb-2555-8b66580eb7fc
# ╟─a8458c44-46c2-11eb-02be-bd3aeb05d50a
# ╠═98ad6aee-46c3-11eb-2659-1df8839b5675
# ╠═d77a0ed8-485a-11eb-0dac-3ff8dee5168a
# ╠═fa5636fc-485a-11eb-2f34-75e68933eb82
# ╠═fdbbaca2-485a-11eb-2705-9370aa177831
# ╠═13e98006-485b-11eb-1d57-ed09c6bde282
# ╠═252120ce-485b-11eb-3637-c14cc1d0153c
# ╠═2bace11c-485b-11eb-01a3-630898d751eb
# ╠═35c07c8e-485b-11eb-24f2-9b5b8cececd6
# ╠═3e922ce8-485b-11eb-2a34-57c74d16a76e
# ╠═52ee9660-485b-11eb-1e3f-9be2725ee183
# ╠═837a9644-485b-11eb-23af-c7439f1b6ac0
# ╠═92ea00b0-485b-11eb-3db8-7d6ec720cfcb
# ╠═2323afe6-485c-11eb-18db-c7f6864abb07
# ╠═33fd5894-485c-11eb-281d-c1fe27b864c1
# ╠═6a165f92-485d-11eb-3c7a-8b6501eccc03
# ╠═c9ba62ea-485d-11eb-2e1c-a58ff592ca6c
# ╠═f2741294-485d-11eb-15db-29302d1a7f4c
# ╠═2a69538a-485e-11eb-2d4b-4b29a9bdd635
# ╠═355e34d6-485e-11eb-3a15-939e18fc8978


text:
#= test/id.jl
=#

  @testset "id" begin
    m = 128
    n =  64
    M = matrixlib(:fourier, rand(m), rand(n))
    opts = LRAOptions(maxdet_tol=0., sketch_randn_niter=1)

    @time for (t, s) in ((:none,                 :none ),
                         (:RandomGaussian,       :randn),
                         (:RandomSubset,         :sub  ),
                         (:SRFT,                 :srft ),
                         (:SparseRandomGaussian, :sprn ))
      opts.sketch = s
      for T in (Float32, Float64, ComplexF32, ComplexF64)
        let A
            println("  $t/$T")

            rtol = 5*eps(real(T))
            approx_rtol = 100*rtol
            opts.rtol = rtol

            A = convert(Array{T}, T <: Real ? real(M) : M)
            nrm = norm(A)

            V  = idfact(:n, A, opts)
            F  = ID(:n, A, V)
            Vc = idfact(:c, A, opts)
            Fc = ID(:c, A, Vc)
            @test norm(A  - Matrix(F )) < approx_rtol*nrm
            @test norm(A' - Matrix(Fc)) < approx_rtol*nrm

            for (N, G, ρ, q) in ((Matrix(V), V, V[:k], n), (A, F, m, n))
              xp = rand(T, ρ)
              xq = rand(T, q)
              y = N  *xq; @test norm(y - G  *xq) < approx_rtol*norm(y)
              y = N' *xp; @test norm(y - G' *xp) < approx_rtol*norm(y)
              y = transpose(N)*xp; @test norm(y - transpose(G)*xp) < approx_rtol*norm(y)

              Bp = rand(T, ρ, ρ)
              Bq = rand(T, q, q)
              C = N   *Bq  ; @test norm(C - G   *Bq  ) < approx_rtol*norm(C)
              C = N   *Bq' ; @test norm(C - G   *Bq' ) < approx_rtol*norm(C)
              C = N   *transpose(Bq); @test norm(C - G   *transpose(Bq)) < approx_rtol*norm(C)
              C = N'  *Bp  ; @test norm(C - G'  *Bp  ) < approx_rtol*norm(C)
              C = N'  *Bp' ; @test norm(C - G'  *Bp' ) < approx_rtol*norm(C)
              C = transpose(N) *Bp  ; @test norm(C - transpose(G) *Bp  ) < approx_rtol*norm(C)
              C = transpose(N) *transpose(Bp); @test norm(C - transpose(G) *transpose(Bp)) < approx_rtol*norm(C)
              C = Bp  *N   ; @test norm(C - Bp  *G   ) < approx_rtol*norm(C)
              C = Bq  *N'  ; @test norm(C - Bq  *G'  ) < approx_rtol*norm(C)
              C = Bq  *transpose(N) ; @test norm(C - Bq  *transpose(G) ) < approx_rtol*norm(C)
              C = Bp' *N   ; @test norm(C - Bp' *G   ) < approx_rtol*norm(C)
              C = Bq' *N'  ; @test norm(C - Bq' *G'  ) < approx_rtol*norm(C)
              C = transpose(Bp)*N   ; @test norm(C - transpose(Bp)*G   ) < approx_rtol*norm(C)
              C = transpose(Bq)*transpose(N) ; @test norm(C - transpose(Bq)*transpose(G) ) < approx_rtol*norm(C)
            end
        end
      end
    end
  end

ext: jl


text:
import QuantumOpticsBase: Operator, SuperOperator, Basis
import LinearAlgebra: ishermitian

ishermitian(A::SuperOperator{B1,B2}) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis}} =
    ishermitian(A.data)

"""
    Commutator(A)

Create commutator in a form of superoperator from a given operator.

``\\text{Commutator}(A) \\:\\cdot\\: = [A, \\:\\cdot\\:]``

"""
function Commutator(A::Operator)::SuperOperator
    return spre(A) - spost(A)
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems


text:
export ActorTrait, InvalidActorTrait, BaseActorTrait, NextActorTrait, ErrorActorTrait, CompletionActorTrait
export AbstractActor, Actor, NextActor, ErrorActor, CompletionActor
export AbstractActorFactory, create_actor
export next!, error!, complete!
export on_next!, on_error!, on_complete!
export as_actor, actor_extract_type

export InvalidActorTraitUsageError, InconsistentSourceActorDataTypesError
export MissingDataArgumentInNextCall, MissingErrorArgumentInErrorCall, ExtraArgumentInCompleteCall
export MissingOnNextImplementationError, MissingOnErrorImplementationError, MissingOnCompleteImplementationError
export MissingCreateActorFactoryImplementationError

import Base: show, showerror
import Base: eltype

"""
Abstract type for all possible actor traits

See also: [`BaseActorTrait`](@ref), [`NextActorTrait`](@ref), [`ErrorActorTrait`](@ref), [`CompletionActorTrait`](@ref), [`InvalidActorTrait`](@ref)
"""
abstract type ActorTrait end

"""
Base actor trait specifies actor to listen for all `next!`, `error!` and `complete!` events.
`BaseActorTrait` is a subtype of `ActorTrait`.

See also: [`Actor`](@ref)
"""
struct BaseActorTrait{T} <: ActorTrait end

"""
Next actor trait specifies actor to listen for `next!` events only.
`NextActorTrait` is a subtype of `ActorTrait`.

See also: [`NextActor`](@ref)
"""
struct NextActorTrait{T} <: ActorTrait end

"""
Error actor trait specifies actor to listen for `error!` events only.
`ErrorActorTrait` is a subtype of `ActorTrait`.

See also: [`ErrorActor`](@ref)
"""
struct ErrorActorTrait{T} <: ActorTrait end

"""
Completion actor trait specifies actor to listen for `complete!` events only.
`CompletionActorTrait` is a subtype of `ActorTrait`.

See also: [`CompletionActor`](@ref)
"""
struct CompletionActorTrait{T} <: ActorTrait end

"""
Default actor trait behavior for any object. Actor with such a trait specificaion cannot be used as a valid actor in `subscribe!` function.
Doing so will raise an error. `InvalidActorTrait` is a subtype of `ActorTrait`.

See also: [`ActorTrait`](@ref)
"""
struct InvalidActorTrait <: ActorTrait end

"""
Supertype type for `Actor`, `NextActor`, `ErrorActor` and `CompletionActor` types.

See also: [`Actor`](@ref), [`NextActor`](@ref), [`ErrorActor`](@ref), [`CompletionActor`](@ref)
"""
abstract type AbstractActor{T} end

"""
Can be used as a super type for common actor. Automatically specifies a `BaseActorTrait` trait behavior.
Every `Actor` must implement its own methods for `on_next!(actor, data)`, `on_error!(actor, err)` and `on_complete!(actor)` functions.
`Actor` is a subtype of `AbstractActor` type.

# Examples
```jldoctest
using Rocket

struct MyActor <: Actor{String} end

Rocket.as_actor(MyActor)

# output

BaseActorTrait{String}()
```

See also: [`AbstractActor`](@ref), [`as_actor`](@ref), [`BaseActorTrait`](@ref), [`ActorTrait`](@ref), [`on_next!`](@ref), [`on_error!`](@ref), [`on_complete!`](@ref)
"""
abstract type Actor{T} <: AbstractActor{T} end

"""
Can be used as a super type for "next-only" actor. Automatically specifies a `NextActorTrait` trait behavior.
Every `NextActor` must implement its own methods for `on_next!(actor, data)` function only.
`NextActor` is a subtype of `AbstractActor` type.

# Examples
```jldoctest
using Rocket

struct MyNextActor <: NextActor{String} end

Rocket.as_actor(MyNextActor)

# output

NextActorTrait{String}()
```

See also: [`AbstractActor`](@ref), [`as_actor`](@ref), [`NextActorTrait`](@ref), [`ActorTrait`](@ref), [`on_next!`](@ref)
"""
abstract type NextActor{T} <: AbstractActor{T} end

"""
Can be used as a super type for "error-only" actor. Automatically specifies a `ErrorActorTrait` trait behavior.
Every `ErrorActor` must implement its own methods for `on_error!(actor, err)` function only.
`ErrorActor` is a subtype of `AbstractActor` type.

# Examples
```jldoctest
using Rocket

struct MyErrorActor <: ErrorActor{String} end

Rocket.as_actor(MyErrorActor)

# output

ErrorActorTrait{String}()
```

See also: [`AbstractActor`](@ref), [`as_actor`](@ref), [`ErrorActorTrait`](@ref), [`ActorTrait`](@ref), [`on_error!`](@ref)
"""
abstract type ErrorActor{T} <: AbstractActor{T} end

"""
Can be used as a super type for "completion-only" actor. Automatically specifies a `CompletionActorTrait` trait behavior.
Every `CompletionActor` must implement its own methods for `on_complete!(actor)` function only.
`CompletionActor` is a subtype of `AbstractActor` type.

# Examples
```jldoctest
using Rocket

struct MyCompletionActor <: CompletionActor{String} end

Rocket.as_actor(MyCompletionActor)

# output

CompletionActorTrait{String}()
```

See also: [`AbstractActor`](@ref), [`as_actor`](@ref), [`CompletionActorTrait`](@ref), [`ActorTrait`](@ref), [`on_complete!`](@ref)
"""
abstract type CompletionActor{T} <: AbstractActor{T} end

"""
    as_actor(any)

This function checks actor trait behavior specification. May be used explicitly to specify actor trait behavior for any object.

See also: [`ActorTrait`](@ref)
"""
as_actor(::Type)                                  = InvalidActorTrait()
as_actor(::Type{ <: Actor{T} })           where T = BaseActorTrait{T}()
as_actor(::Type{ <: NextActor{T} })       where T = NextActorTrait{T}()
as_actor(::Type{ <: ErrorActor{T} })      where T = ErrorActorTrait{T}()
as_actor(::Type{ <: CompletionActor{T} }) where T = CompletionActorTrait{T}()
as_actor(::O)                             where O = as_actor(O)

actor_extract_type(::Type{A}) where A = actor_extract_type(as_actor(A), A)
actor_extract_type(::A)       where A = actor_extract_type(A)

actor_extract_type(::BaseActorTrait{T}, _)       where T = T
actor_extract_type(::NextActorTrait{T}, _)       where T = T
actor_extract_type(::ErrorActorTrait{T}, _)      where T = T
actor_extract_type(::CompletionActorTrait{T}, _) where T = T
actor_extract_type(::InvalidActorTrait,  actor)          = throw(InvalidActorTraitUsageError(actor))

Base.eltype(actor::Actor)           = actor_extract_type(actor)
Base.eltype(actor::NextActor)       = actor_extract_type(actor)
Base.eltype(actor::ErrorActor)      = actor_extract_type(actor)
Base.eltype(actor::CompletionActor) = actor_extract_type(actor)

Base.eltype(actor::Type{ <: Actor })           = actor_extract_type(actor)
Base.eltype(actor::Type{ <: NextActor })       = actor_extract_type(actor)
Base.eltype(actor::Type{ <: ErrorActor })      = actor_extract_type(actor)
Base.eltype(actor::Type{ <: CompletionActor }) = actor_extract_type(actor)

next!(_)  = throw(MissingDataArgumentInNextCall())
error!(_) = throw(MissingErrorArgumentInErrorCall())

"""
    next!(actor, data)
    next!(actor, data, scheduler)

This function is used to deliver a "next" event to an actor with some `data`. Takes optional `scheduler` object to schedule execution of data delivery.

See also: [`AbstractActor`](@ref), [`on_next!`](@ref)
"""
next!(actor::T, data)            where T = actor_on_next!(as_actor(T), actor, data)
next!(actor::T, data, scheduler) where T = actor_on_next!(as_actor(T), actor, data, scheduler)

# Specialized methods for default built-in type of actor (more efficient usage of stack and unnecessary checks bypass)
next!(actor::Actor{T}, data::T)            where T = begin on_next!(actor, data); return nothing end
next!(actor::Actor{T}, data::T, scheduler) where T = begin scheduled_next!(actor, data, scheduler); return nothing end
next!(actor::Actor{T}, data::R)            where { T, R } = throw(InconsistentSourceActorDataTypesError{T, R}(actor))
next!(actor::Actor{T}, data::R, scheduler) where { T, R } = throw(InconsistentSourceActorDataTypesError{T, R}(actor))

"""
    error!(actor, err)
    error!(actor, err, scheduler)

This function is used to deliver a "error" event to an actor with some `err`. Takes optional `scheduler` object to schedule execution of error delivery.

See also: [`AbstractActor`](@ref), [`on_error!`](@ref)
"""
error!(actor::T, err)            where T = actor_on_error!(as_actor(T), actor, err)
error!(actor::T, err, scheduler) where T = actor_on_error!(as_actor(T), actor, err, scheduler)

# Specialized methods for default built-in type of actor (more efficient usage of stack and unnecessary checks bypass)
error!(actor::Actor, err)            = begin on_error!(actor, err); return nothing end
error!(actor::Actor, err, scheduler) = begin scheduled_error!(actor, err, scheduler); return nothing end

"""
    complete!(actor)
    complete!(actor, scheduler)

This function is used to deliver a "complete" event to an actor. Takes optional `scheduler` object to schedule execution of complete event delivery.

See also: [`AbstractActor`](@ref), [`on_complete!`](@ref)
"""
complete!(actor::T)            where T = actor_on_complete!(as_actor(T), actor)
complete!(actor::T, scheduler) where T = actor_on_complete!(as_actor(T), actor, scheduler)

# Specialized methods for default built-in type of actor (more efficient usage of stack and unnecessary checks bypass)
complete!(actor::Actor)            = begin on_complete!(actor); return nothing end
complete!(actor::Actor, scheduler) = begin scheduled_complete!(actor, scheduler); return nothing end

actor_on_next!(::InvalidActorTrait,     actor, data) = throw(InvalidActorTraitUsageError(actor))
actor_on_error!(::InvalidActorTrait,    actor, err)  = throw(InvalidActorTraitUsageError(actor))
actor_on_complete!(::InvalidActorTrait, actor)       = throw(InvalidActorTraitUsageError(actor))

actor_on_next!(::BaseActorTrait{T},       actor, data::R) where { R, T } = throw(InconsistentSourceActorDataTypesError{T, R}(actor))
actor_on_next!(::NextActorTrait{T},       actor, data::R) where { R, T } = throw(InconsistentSourceActorDataTypesError{T, R}(actor))
actor_on_next!(::ErrorActorTrait{T},      actor, data::R) where { R, T } = throw(InconsistentSourceActorDataTypesError{T, R}(actor))
actor_on_next!(::CompletionActorTrait{T}, actor, data::R) where { R, T } = throw(InconsistentSourceActorDataTypesError{T, R}(actor))

actor_on_next!(::BaseActorTrait{T},       actor, data::R) where { T, R <: T } = begin on_next!(actor, data); return nothing end
actor_on_next!(::NextActorTrait{T},       actor, data::R) where { T, R <: T } = begin on_next!(actor, data); return nothing end
actor_on_next!(::ErrorActorTrait{T},      actor, data::R) where { T, R <: T } = nothing
actor_on_next!(::CompletionActorTrait{T}, actor, data::R) where { T, R <: T } = nothing

actor_on_next!(::BaseActorTrait{T},       actor, data::R, scheduler) where { T, R <: T } = begin scheduled_next!(actor, data, scheduler); return nothing end
actor_on_next!(::NextActorTrait{T},       actor, data::R, scheduler) where { T, R <: T } = begin scheduled_next!(actor, data, scheduler); return nothing end
actor_on_next!(::ErrorActorTrait{T},      actor, data::R, scheduler) where { T, R <: T } = nothing
actor_on_next!(::CompletionActorTrait{T}, actor, data::R, scheduler) where { T, R <: T } = nothing

actor_on_error!(::BaseActorTrait,       actor, err) = begin on_error!(actor, err); return nothing end
actor_on_error!(::NextActorTrait,       actor, err) = nothing
actor_on_error!(::ErrorActorTrait,      actor, err) = begin on_error!(actor, err); return nothing end
actor_on_error!(::CompletionActorTrait, actor, err) = nothing

actor_on_error!(::BaseActorTrait,       actor, err, scheduler) = begin scheduled_error!(actor, err, scheduler); return nothing end
actor_on_error!(::NextActorTrait,       actor, err, scheduler) = nothing
actor_on_error!(::ErrorActorTrait,      actor, err, scheduler) = begin scheduled_error!(actor, err, scheduler); return nothing end
actor_on_error!(::CompletionActorTrait, actor, err, scheduler) = nothing

actor_on_complete!(::BaseActorTrait,       actor) = begin on_complete!(actor); return nothing end
actor_on_complete!(::NextActorTrait,       actor) = nothing
actor_on_complete!(::ErrorActorTrait,      actor) = nothing
actor_on_complete!(::CompletionActorTrait, actor) = begin on_complete!(actor); return nothing end

actor_on_complete!(::BaseActorTrait,       actor, scheduler) = begin scheduled_complete!(actor, scheduler); return nothing end
actor_on_complete!(::NextActorTrait,       actor, scheduler) = nothing
actor_on_complete!(::ErrorActorTrait,      actor, scheduler) = nothing
actor_on_complete!(::CompletionActorTrait, actor, scheduler) = begin scheduled_complete!(actor, scheduler); return nothing end

"""
    on_next!(actor, data)

Both Actor and NextActor objects must implement its own method for `on_next!` function which will be called on "next" event.

See also: [`Actor`](@ref), [`NextActor`](@ref)
"""
on_next!(actor, data) = throw(MissingOnNextImplementationError(actor, data))

"""
    on_error!(actor, err)

Both Actor and ErrorActor objects must implement its own method for `on_error!` function which will be called on "error" event.

See also: [`Actor`](@ref), [`ErrorActor`](@ref)
"""
on_error!(actor, err) = throw(MissingOnErrorImplementationError(actor))

"""
    on_complete!(actor)

Both Actor and CompletionActor objects must implement its own method for `on_complete!` function which will be called on "complete" event.

See also: [`Actor`](@ref), [`CompletionActor`](@ref)
"""
on_complete!(actor) = throw(MissingOnCompleteImplementationError(actor))

# -------------------------------- #
# Actor factory                    #
# -------------------------------- #

"""
Abstract type for all possible actor factories

See also: [`Actor`](@ref)
"""
abstract type AbstractActorFactory end

"""
    create_actor(::Type{L}, factory::F) where { L, F <: AbstractActorFactory }

Actor creator function for a given factory `F`. Should be implemented explicitly for any `AbstractActorFactory` object

See also: [`AbstractActorFactory`](@ref), [`MissingCreateActorFactoryImplementationError`](@ref)
"""
create_actor(::Type{L}, factory::F) where { L, F <: AbstractActorFactory } = throw(MissingCreateActorFactoryImplementationError(factory))

# -------------------------------- #
# Errors                           #
# -------------------------------- #

"""
This error will be throw if Julia cannot find specific method of 'create_actor()' function for given actor factory

See also: [`AbstractActorFactory`](@ref), [`create_actor`](@ref)
"""
struct MissingCreateActorFactoryImplementationError
    factory
end

function Base.showerror(io::IO, err::MissingCreateActorFactoryImplementationError)
    print(io, "You probably forgot to implement create_actor(::Type{L}, factory::$(typeof(err.factory))).")
end

"""
This error will be thrown if `next!`, `error!` or `complete!` functions are called with invalid actor object

See also: [`next!`](@ref), [`error!`](@ref), [`complete!`](@ref), [`InvalidActorTrait`](@ref)
"""
struct InvalidActorTraitUsageError
    actor
end

function Base.showerror(io::IO, err::InvalidActorTraitUsageError)
    print(io, "Type $(typeof(err.actor)) is not a valid actor type. \nConsider extending your actor with one of the abstract actor types <: (Actor{T}, NextActor{T}, ErrorActor{T}, CompletionActor{T}) or implement as_actor(::Type{<:$(typeof(err.actor))}).")
end

"""
This error will be thrown if `next!` function is called with inconsistent data type

See also: [`AbstractActor`](@ref), [`Subscribable`](@ref), [`next!`](@ref)
"""
struct InconsistentSourceActorDataTypesError{T, R}
    actor
end

function Base.showerror(io::IO, err::InconsistentSourceActorDataTypesError{T, R}) where T where R
    print(io, "Actor of type $(typeof(err.actor)) expects data to be of type $(T), but data of type $(R) has been found.")
end

"""
This error will be thrown if `next!` function is called without data argument

See also: [`next!`](@ref)
"""
struct MissingDataArgumentInNextCall end

function Base.showerror(io::IO, ::MissingDataArgumentInNextCall)
    print(io, "Missing data argument in next! callback")
end

"""
This error will be thrown if `error!` function is called without err argument

See also: [`error!`](@ref)
"""
struct MissingErrorArgumentInErrorCall end

function Base.showerror(io::IO, ::MissingErrorArgumentInErrorCall)
    print(io, "Missing err argument in error! callback")
end

"""
This error will be thrown if `complete!` function is called with extra data/err argument

See also: [`complete!`](@ref)
"""
struct ExtraArgumentInCompleteCall end

function Base.showerror(io::IO, ::ExtraArgumentInCompleteCall)
    print(io, "Extra argument in complete! callback")
end

"""
This error will be thrown if Julia cannot find specific method of 'on_next!()' function for given actor and data

See also: [`on_next!`](@ref)
"""
struct MissingOnNextImplementationError
    actor
    data
end

function Base.showerror(io::IO, err::MissingOnNextImplementationError)
    print(io, "You probably forgot to implement on_next!(actor::$(typeof(err.actor)), data::$(typeof(err.data))).")
end

"""
This error will be thrown if Julia cannot find specific method of 'on_error!()' function for given actor

See also: [`on_error!`](@ref)
"""
struct MissingOnErrorImplementationError
    actor
end

function Base.showerror(io::IO, err::MissingOnErrorImplementationError)
    print(io, "You probably forgot to implement on_error!(actor::$(typeof(err.actor)), err).")
end

"""
This error will be thrown if Julia cannot find specific method of 'on_complete!()' function for given actor and data

See also: [`on_next!`](@ref)
"""
struct MissingOnCompleteImplementationError
    actor
end

function Base.showerror(io::IO, err::MissingOnCompleteImplementationError)
    print(io, "You probably forgot to implement on_complete!(actor::$(typeof(err.actor))).")
end


text:
module GtkTextUtils

    using Gtk
    using ..JuliaWordsUtils
    import Gtk: GtkTextIter, GLib.mutable, GLib.Mutable, hasselection

    export select_word_double_click, select_word, get_line_text, cursor_position,
    get_text_iter_at_cursor, get_current_line_text,
    get_text_left_of_iter, get_text_right_of_iter,
    move_cursor_to_sentence_start, move_cursor_to_sentence_end,
    hasselection, replace_text, offset

    offset(it::Gtk.TI) = get_gtk_property(it, :offset, Integer)
    line(it::Gtk.TI) = get_gtk_property(it,:line,Integer)+1#Gtk counts from zero

    function hasselection(b::GtkTextBuffer)
        (found, it_start, it_end) = selection_bounds(b)
        found
    end
    function replace_text(buffer::GtkTextBuffer, itstart::T, itend::T, str::AbstractString) where {T<:Gtk.TI}
        pos = offset(itstart)+1
        splice!(buffer, itstart:itend)
        insert!(buffer, GtkTextIter(buffer, pos), str)
    end

    cursor_position(b::GtkTextBuffer) = get_gtk_property(b, :cursor_position, Int)

    get_text_iter_at_cursor(b::GtkTextBuffer) =
        GtkTextIter(b, cursor_position(b)+1) #+1 because there's a -1 in gtk.jl

    function get_current_line_text(buffer::GtkTextBuffer)
        it = get_text_iter_at_cursor(buffer)
        return get_line_text(buffer, it)
    end

    function select_word(it::GtkTextIter, buffer::GtkTextBuffer, stop_at_dot::Bool)

        (txt, line_start, line_end) = get_line_text(buffer, it)

        pos = offset(it) - offset(line_start) +1#not sure about the +1 but it feels better
        if pos <= 0
            return ("", GtkTextIter(buffer, offset(it)),
            GtkTextIter(buffer, offset(it)))
        end

        word, i, j = extend_word(pos, txt, stop_at_dot)

        its = GtkTextIter(buffer, i + offset(line_start) )
        ite = GtkTextIter(buffer, j + offset(line_start) + 1)

        return (word, its, ite)
    end
    select_word(it::GtkTextIter, buffer::GtkTextBuffer) = select_word(it, buffer, true)

    function JuliaWordsUtils.select_word_backward(it::GtkTextIter, buffer::GtkTextBuffer, stop_at_dot::Bool)

        (txt, line_start, line_end) = get_line_text(buffer, it)
        pos = offset(it) - offset(line_start) #position of cursor in txt

        if pos <= 0 || length(txt) == 0
            return ("", GtkTextIter(buffer, offset(it)),
            GtkTextIter(buffer, offset(it)))
        end

        txt = CharArray(txt, pos)
        (i, j) = select_word_backward(pos, txt, stop_at_dot)

        its = GtkTextIter(buffer, i + offset(line_start) )
        ite = GtkTextIter(buffer, offset(it))

        return (txt[i:j], its, it)
    end

    function get_line_text(buffer::GtkTextBuffer, it::GtkTextIter)

        itstart, itend = mutable(it), mutable(it)
        li = get_gtk_property(itstart, :line, Integer)

        skip(itstart, :backward_line) #seems there's no skip to line start
        li != get_gtk_property(itstart, :line, Integer) && skip(itstart, 1, :line)#for fist line
        !get_gtk_property(itend, :ends_line, Bool) && skip(itend, :forward_to_line_end)

        return ( (itstart:itend).text[String], itstart, itend)
    end

    function select_word_double_click(textview::GtkTextView, buffer::GtkTextBuffer, x::Integer, y::Integer)

        x, y = Gtk.window_to_buffer_coords(textview, x, y)
        iter_end = Gtk.text_iter_at_position(textview, x, y)
        #iter_end = mutable( get_text_iter_at_cursor(buffer) ) #not using this because the cursor position is modified somewhere

        (w, iter_start, iter_end) = select_word(iter_end, buffer)
        select_range(buffer, iter_start, iter_end)
    end

    function get_text_right_of_cursor(buffer::GtkTextBuffer)
        it = mutable(get_text_iter_at_cursor(buffer))
        return (it:(it+1)).text[String]
    end
    function get_text_left_of_cursor(buffer::GtkTextBuffer)
        it = mutable(get_text_iter_at_cursor(buffer))
        return ((it-1):it).text[String]
    end

    get_text_left_of_iter(it::Mutable{GtkTextIter}) = ((it-1):it).text[String]
    get_text_right_of_iter(it::Mutable{GtkTextIter}) = (it:(it+1)).text[String]

    get_text_left_of_iter(it::GtkTextIter) = ((mutable(it)-1):mutable(it)).text[String]
    get_text_right_of_iter(it::GtkTextIter) = (mutable(it):(mutable(it)+1)).text[String]

    nlines(it_start, it_end) = abs(line(it_end)-line(it_start))+1

    function move_cursor_to_sentence_start(buffer::GtkTextBuffer)
        it = mutable( get_text_iter_at_cursor(buffer) )
        skip(it, :backward_sentence_start)
        place_cursor(buffer, it)
    end
    function move_cursor_to_sentence_end(buffer::GtkTextBuffer)
        it = mutable( get_text_iter_at_cursor(buffer) )
        skip(it, :forward_sentence_end)
        place_cursor(buffer, it)
    end

    function toggle_wrap_mode(v::GtkTextView)
        wm = get_gtk_property(v, :wrap_mode, Int)
        wm = convert(Bool, wm)
        setproperty!(v, :wrap_mode, !wm)
        nothing
    end

end # module


text:
module GalewskyThermalShallowWaterExplicit

using Gridap
using GridapGeosciences

include("GalewskyInitialConditions.jl")

# Solves the Galewsky test case for the shallow water equations on a sphere
# of physical radius 6371220m. Involves a shear flow instability of a zonal
# jet triggered by an initial gravity wave.
# reference:
#   Galewsky, Scott and Polvani (2004) Tellus, 56A 429-440

# initialise the depth weighted buoyancy field such that the buoyancy is just
# a constant gravity field, such that in the continuous form the thermal shallow
# water equations default to just the shallow water equations
function gh₀(xyz)
  h = h₀(xyz)
  E = g*h
  E
end

order  = 1
degree = 4

n      = 48
nstep  = 20*24*60 # 20 days
dt     = 60.0

model = CubedSphereDiscreteModel(n; radius=rₑ)

hf, uf = thermal_shallow_water_explicit_time_stepper(model, order, degree,
                                    h₀, u₀, gh₀, f,
                                    dt, 0.5*dt, nstep;
                                    write_solution=true,
                                    write_solution_freq=240,
                                    write_diagnostics=true,
                                    write_diagnostics_freq=1,
                                    dump_diagnostics_on_screen=true)

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

# The code here is mostly copied from the flip_sign.jl code for FlipSignBridge and GreaterToLessBridge

"""
    AbstractToIntervalBridge{T, S1, F}

Bridge a `F`-in-`Interval` constraint into an `F`-in-`Interval{T}` constraint where we have either:
* `S1 = MOI.GreaterThan{T}`
* `S1 = MOI.LessThan{T}`

The `F`-in-`Interval{T}` constraint is stored in the `constraint`
field by convention.

!!! warning
    It is required that `T` be a `AbstractFloat` type because otherwise
    `typemin` and `typemax` would either be not implemented (e.g. `BigInt`)
    or would not give infinite value (e.g. `Int`). For this reason,
    this bridge is only added to
    [`MathOptInterface.Bridges.full_bridge_optimizer`](@ref).

    when `T` is a subtype of `AbstractFloat`.
"""
abstract type AbstractToIntervalBridge{
    T<:AbstractFloat,
    S1<:MOI.AbstractScalarSet,
    F<:MOI.AbstractScalarFunction,
} <: SetMapBridge{T,MOI.Interval{T},S1,F,F} end

# The function map is the identity. It is also an involution, symmetric, and a symmetric involution.
MOIB.map_function(::Type{<:AbstractToIntervalBridge{T}}, func) where {T} = func
MOIB.inverse_map_function(::Type{<:AbstractToIntervalBridge}, func) = func
MOIB.adjoint_map_function(::Type{<:AbstractToIntervalBridge}, func) = func
function MOIB.inverse_adjoint_map_function(
    ::Type{<:AbstractToIntervalBridge},
    func,
)
    return func
end

# FIXME are these modify functions necessary?
function MOI.modify(
    model::MOI.ModelLike,
    bridge::AbstractToIntervalBridge,
    change::MOI.ScalarCoefficientChange,
)
    MOI.modify(model, bridge.constraint, change)
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::AbstractToIntervalBridge,
    change::MOI.MultirowChange{T},
) where {T}
    MOI.modify(model, bridge.constraint, change)
    return
end

"""
    GreaterToIntervalBridge{T, F<:MOI.AbstractScalarFunction} <:
        AbstractToIntervalBridge{T, MOI.GreaterThan{T}, F}

Transforms a `F`-in-`GreaterThan{T}` constraint into an `F`-in-`Interval{T}`
constraint.
"""
struct GreaterToIntervalBridge{T,F<:MOI.AbstractScalarFunction} <:
       AbstractToIntervalBridge{T,MOI.GreaterThan{T},F}
    constraint::MOI.ConstraintIndex{F,MOI.Interval{T}}
end

const GreaterToInterval{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{GreaterToIntervalBridge{T},OT}

function MOIB.map_set(::Type{<:GreaterToIntervalBridge}, set::MOI.GreaterThan)
    return MOI.Interval(set.lower, typemax(set.lower))
end

function MOIB.inverse_map_set(
    ::Type{<:GreaterToIntervalBridge},
    set::MOI.Interval,
)
    return MOI.GreaterThan(set.lower)
end

function concrete_bridge_type(
    ::Type{<:GreaterToIntervalBridge{T}},
    F::Type{<:MOI.AbstractScalarFunction},
    ::Type{MOI.GreaterThan{T}},
) where {T}
    return GreaterToIntervalBridge{T,F}
end

"""
    LessToIntervalBridge{T, F<:MOI.AbstractScalarFunction} <:
        AbstractToIntervalBridge{T, MOI.LessThan{T}, F}

Transforms a `F`-in-`LessThan{T}` constraint into an `F`-in-`Interval{T}`
constraint.
"""
struct LessToIntervalBridge{T,F<:MOI.AbstractScalarFunction} <:
       AbstractToIntervalBridge{T,MOI.LessThan{T},F}
    constraint::MOI.ConstraintIndex{F,MOI.Interval{T}}
end

const LessToInterval{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{LessToIntervalBridge{T},OT}

function MOIB.map_set(::Type{<:LessToIntervalBridge}, set::MOI.LessThan)
    return MOI.Interval(typemin(set.upper), set.upper)
end

function MOIB.inverse_map_set(::Type{<:LessToIntervalBridge}, set::MOI.Interval)
    return MOI.LessThan(set.upper)
end

function concrete_bridge_type(
    ::Type{<:LessToIntervalBridge{T}},
    F::Type{<:MOI.AbstractScalarFunction},
    ::Type{MOI.LessThan{T}},
) where {T}
    return LessToIntervalBridge{T,F}
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using Simulacrum
using Symbolics
import SymbolicUtils: Term, Sym

@test Simulacrum.get_lowest_subtype(Term, Term{Real, Nothing}) == Term{Real, Nothing}
@test Simulacrum.get_lowest_subtype(Sym{Real, Nothing}, Sym) == Sym{Real, Nothing}
@test Simulacrum.get_lowest_subtype(Sym{Real, Nothing}, Sym, Sym{Real}) == Sym{Real, Nothing}

@variables t x(t) y
D = Differential(t)
z = t + t^2
diff = D(z)
eq = Equation(diff, x)

# code below is only useful if this is true
# @test typeof(x) == typeof(diff) == typeof(eq.lhs) == typeof(eq.rhs) == Num
# update: Equation types now come unwrapped out of the box; still keeping these
# in case that changes again in the future

# testing unwrapping of Num via get_value
conv_eq = Simulacrum.get_value(eq)
@test eq == conv_eq
@test typeof(conv_eq.lhs) <: Term{Real} && typeof(conv_eq.rhs) <: Term{Real}
eq = Equation(diff, y)
conv_eq = Simulacrum.get_value(eq)
@test typeof(conv_eq.rhs) <: Sym{Real}

@test Simulacrum.get_sym_name(Symbolics.variable(:x, T = Symbolics.FnType)(t)) == :x # Term
@test Simulacrum.get_sym_name(Symbolics.variable(:x)) == :x # Sym
@test Simulacrum.get_sym_name(Num(Symbolics.variable(:x))) == :x # Num

# make_var (& gen_var)
@variables t
make_var(:test1)
make_var(:test2, 1)
make_var(:test3, t)
make_var(:test4, 1, t)
@test_throws ErrorException make_var(:test5, make_var(:somevar, t))
@test_throws MethodError make_var(:test6, t, t)
x = make_var(:x)
@test x === make_var(:x)
y = make_var(:y, t)
@test isequal(y, make_var(:y, t))

# subscript_var: Term, Sym, and Num types
a = Symbolics.variable(:x, T = Symbolics.FnType)(t)
b = Symbolics.variable(:t)
c = (@variables x(t))[1]
typeof(a) == Term{Real, Nothing}
typeof(b) == Sym{Real, Nothing}
typeof(c) == Num
@test typeof(Simulacrum.subscript_var(a, 1)) == typeof(a)
@test typeof(Simulacrum.subscript_var(b, 1)) == typeof(b)
@test typeof(Simulacrum.subscript_var(c, 1)) == typeof(c)

# subscript_var and subscript_vars
a = Simulacrum.subscript_vars([Symbolics.variable(:x)], 1)[1]
b = Simulacrum.subscript_var(Symbolics.variable(:x), 1)
@test isequal(a, b)

# subscript_var: Num type
@variables x(t) y
@test typeof(Simulacrum.subscript_var(x, 1)) == typeof(x) == Num
@test typeof(Simulacrum.get_value(Simulacrum.subscript_var(x, 1))) == typeof(Simulacrum.get_value(x))
@test typeof(Simulacrum.subscript_var(y, 1)) == typeof(y) == Num

# wrapping Symbolic variables into Num type
@variables t X(t) Y(t)
r = Reaction(2, [X], [Y])
che = convert(ChemicalHyperEdge, r)
@test eltype(convert(ChemicalHyperEdge{Num}, che)) == Num
chx = ChemicalHyperGraph(che)
@test eltype(convert(ChemicalHyperGraph{Num}, chx)) == Num

# converting ordinary edge indices to hyperedge indices for patterns
n = rand((1:1000))
idx = indices(AllToAll(n))
@test length(Simulacrum.edge_idx_to_hyperedge_idx(idx)) == n
idx = indices(Cycle(n))
@test length(Simulacrum.edge_idx_to_hyperedge_idx(idx)) == length(idx)
idx = indices(Path(n))
@test length(Simulacrum.edge_idx_to_hyperedge_idx(idx)) == length(idx)


text:
using Random
using StaticArrays
using Test

using CLIMA
using CLIMA.Atmos
using CLIMA.ConfigTypes
using CLIMA.Diagnostics
using CLIMA.GenericCallbacks
using CLIMA.ODESolvers
using CLIMA.Mesh.Filters
using CLIMA.MoistThermodynamics
using CLIMA.VariableTemplates

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# ------------------------ Description ------------------------- #
# 1) Dry Rising Bubble (circular potential temperature perturbation)
# 2) Boundaries - `All Walls` : Impenetrable(FreeSlip())
#                               Laterally periodic
# 3) Domain - 2500m[horizontal] x 2500m[horizontal] x 2500m[vertical]
# 4) Timeend - 1000s
# 5) Mesh Aspect Ratio (Effective resolution) 1:1
# 7) Overrides defaults for
#               `init_on_cpu`
#               `solver_type`
#               `sources`
#               `C_smag`
# 8) Default settings can be found in `src/Driver/Configurations.jl`
# ------------------------ Description ------------------------- #
function init_risingbubble!(bl, state, aux, (x, y, z), t)
    FT = eltype(state)
    R_gas::FT = R_d(bl.param_set)
    c_p::FT = cp_d(bl.param_set)
    c_v::FT = cv_d(bl.param_set)
    p0::FT = MSLP(bl.param_set)
    _grav::FT = grav(bl.param_set)
    γ::FT = c_p / c_v

    xc::FT = 1250
    yc::FT = 1250
    zc::FT = 1000
    r = sqrt((x - xc)^2 + (y - yc)^2 + (z - zc)^2)
    rc::FT = 500
    θ_ref::FT = 300
    Δθ::FT = 0

    if r <= rc
        Δθ = FT(5) * cospi(r / rc / 2)
    end

    #Perturbed state:
    θ = θ_ref + Δθ # potential temperature
    π_exner = FT(1) - _grav / (c_p * θ) * z # exner pressure
    ρ = p0 / (R_gas * θ) * (π_exner)^(c_v / R_gas) # density
    q_tot = FT(0)
    ts = LiquidIcePotTempSHumEquil(θ, ρ, q_tot, bl.param_set)
    q_pt = PhasePartition(ts)

    ρu = SVector(FT(0), FT(0), FT(0))

    #State (prognostic) variable assignment
    e_kin = FT(0)
    e_pot = gravitational_potential(bl.orientation, aux)
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)
    state.ρ = ρ
    state.ρu = ρu
    state.ρe = ρe_tot
    state.moisture.ρq_tot = ρ * q_pt.tot
end

function config_risingbubble(FT, N, resolution, xmax, ymax, zmax)

    # Choose explicit solver
    ode_solver = CLIMA.MultirateSolverType(
        linear_model = AtmosAcousticGravityLinearModel,
        slow_method = LSRK144NiegemannDiehlBusch,
        fast_method = LSRK144NiegemannDiehlBusch,
        timestep_ratio = 10,
    )

    # Set up the model
    C_smag = FT(0.23)
    ref_state =
        HydrostaticState(DryAdiabaticProfile(typemin(FT), FT(300)), FT(0))
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        param_set;
        turbulence = SmagorinskyLilly(C_smag),
        hyperdiffusion = StandardHyperDiffusion(60),
        source = (Gravity(),),
        ref_state = ref_state,
        init_state = init_risingbubble!,
    )

    # Problem configuration
    config = CLIMA.AtmosLESConfiguration(
        "DryRisingBubble",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_risingbubble!,
        solver_type = ode_solver,
        model = model,
    )
    return config
end

function config_diagnostics(driver_config)
    interval = 10000 # in time steps
    dgngrp = setup_atmos_default_diagnostics(interval, driver_config.name)
    return CLIMA.DiagnosticsConfiguration([dgngrp])
end

function main()
    CLIMA.init()

    # Working precision
    FT = Float64
    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δh = FT(50)
    Δv = FT(50)
    resolution = (Δh, Δh, Δv)
    # Domain extents
    xmax = FT(2500)
    ymax = FT(2500)
    zmax = FT(2500)
    # Simulation time
    t0 = FT(0)
    timeend = FT(1000)

    # Courant number
    CFL = FT(20)

    driver_config = config_risingbubble(FT, N, resolution, xmax, ymax, zmax)
    solver_config = CLIMA.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFL,
    )
    dgn_config = config_diagnostics(driver_config)

    # User defined filter (TMAR positivity preserving filter)
    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do (init = false)
        Filters.apply!(solver_config.Q, 6, solver_config.dg.grid, TMARFilter())
        nothing
    end

    # Invoke solver (calls solve! function for time-integrator)
    result = CLIMA.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbtmarfilter,),
        check_euclidean_distance = true,
    )

    @test isapprox(result, FT(1); atol = 1.5e-3)
end

main()

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
#=
  GLM Family Functions
=#

#=
  References:

  1. Generalized Linear Models & Extensions: J. W. Hardin & J. M. Hilbe
=#

#=
  Types for Distribution Functions
=#
abstract type AbstractDistribution end
struct BernoulliDistribution <: AbstractDistribution end
struct BinomialDistribution <: AbstractDistribution end
struct GammaDistribution <: AbstractDistribution end
struct PoissonDistribution <: AbstractDistribution end
struct GaussianDistribution <: AbstractDistribution end
struct InverseGaussianDistribution <: AbstractDistribution end
struct NegativeBernoulliDistribution{T <: AbstractFloat} <: AbstractDistribution
  Alpha::T
end
struct PowerDistribution{T <: AbstractFloat} <: AbstractDistribution
  k::T
end
# struct LogNormalDistribution <: AbstractDistribution end

#====================================================================
  Copy constructors for distribution types
====================================================================#
import Base.copy
copy(::BernoulliDistribution) = BernoulliDistribution()
copy(::BinomialDistribution) = BinomialDistribution()
copy(::GammaDistribution) = GammaDistribution()
copy(::PoissonDistribution) = PoissonDistribution()
copy(::GaussianDistribution) = GaussianDistribution()
copy(::InverseGaussianDistribution) = InverseGaussianDistribution()
function copy(distrib::NegativeBernoulliDistribution{T}) where {T}
  return NegativeBernoulliDistribution(distrib.Alpha)
end
function copy(distrib::PowerDistribution{T}) where {T}
  return PowerDistribution(distrib.k)
end
#===================================================================#

# Type returned from init!() function
InitType{T} = Tuple{Array{T, 1}, Array{T, 1}, Array{T, 1}} # where {T <: AbstractFloat}

#=
  Distribution Implementations
=#

# Default Initializer for all the distributions
@inline function init!(::AbstractDistribution, y::Array{T}, wts::Array{T, 1})::InitType{T} where {T <: AbstractFloat}
  y = y[:, 1]
  return y, y, wts
end

function init!(::AbstractDistribution, y::Array{Array{T, 2}, 1}, wts::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  if length(size(y[1])) == 2
    y = [y[i][:, 1] for i in 1:nBlocks]
  end
  mu = Array{Array{T, 1}, 1}(undef, nBlocks);
  for i in 1:nBlocks
    mu[i] = copy(y[i]);
  end
  return y, mu, wts;
end

function init!(::Block1DParallel, ::AbstractDistribution, y::Array{Array{T, 2}, 1}, wts::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  isMat::Bool = length(size(y[1])) == 2
  mu = Array{Array{T, 1}, 1}(undef, nBlocks);
  yNew = Array{Array{T, 1}, 1}(undef, nBlocks);
  @threads for i in 1:nBlocks
    tmp = isMat ? y[i][:, 1] : y[i]
    mu[i] = copy(tmp);
    yNew[i] = tmp
  end
  return yNew, mu, wts;
end

#=
  BinomialDistribution
=#

# Initializer
@inline function init!(::BinomialDistribution, y::Array{T}, wts::Array{T, 1})::InitType{T} where {T <: AbstractFloat}
  if size(y)[2] == 1
    y = y[:, 1]
    if length(wts) == 0
      mu = (y .+ T(0.5))./2
    else
      mu = (wts .* y .+ T(0.5))./(wts .+ T(1))
    end
  elseif size(y)[2] == 2
    n = size(y)[1]
    events = y[:, 1]
    N = events .+ y[:, 2]
    y = [if N[i] != 0; events[i]/N[i]; else T(0) end for i in 1:n]
    if length(wts) != 0
      wts .*= N
    else
      wts = N
    end
    mu = (N .* y .+ T(0.5))./(N .+ T(1))
  else
    error("There was a problem")
  end
  return y, mu, wts
end

@inline function init!(::BinomialDistribution, y::Array{Array{T, 2}, 1}, wts::Array{Array{T, 1}, 1})::InitType{Array{T, 1}} where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  if size(y[1])[2] == 1
    y = [y[i][:, 1] for i in 1:nBlocks]
    if length(wts) == 0
      mu = [(y[i] .+ T(0.5))./2 for i in 1:nBlocks]
    else
      mu = [(wts[i] .* y .+ T(0.5))./(wts[i] .+ T(1)) for i in 1:nBlocks]
    end
  elseif size(y[1])[2] == 2
    events = [y[i][:, 1] for i in 1:nBlocks]
    N = [events[i] .+ y[i][:, 2] for i in 1:nBlocks]
    yNew = Array{Array{T, 1}, 1}(undef, nBlocks)
    for i in 1:nBlocks
      n = size(y[i])[1]
      tmp = zeros(T, n)
      for j in 1:n
        tmp[j] = if N[i][j] != 0; events[i][j]/N[i][j]; else T(0) end
      end
      yNew[i] = tmp
    end
    y = yNew
    # y = [if N[i] != 0; events[i]/N[i]; else T(0) end for i in 1:n]
    if length(wts) != 0
      [wts[i] .*= N[i] for i in 1:nBlocks]
    else
      wts = N
    end
    mu = [(N[i] .* y[i] .+ T(0.5))./(N[i] .+ T(1)) for i in 1:nBlocks]
  else
    error("There was a problem with the dimensions of y")
  end
  return y, mu, wts
end

@inline function init!(::Block1DParallel, ::BinomialDistribution, y::Array{Array{T, 2}, 1}, wts::Array{Array{T, 1}, 1})::InitType{Array{T, 1}} where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  yNew = Array{Array{T, 1}, 1}(undef, nBlocks)
  mu = Array{Array{T, 1}, 1}(undef, nBlocks)
  events = Array{Array{T, 1}, 1}(undef, nBlocks)
  N = Array{Array{T, 1}, 1}(undef, nBlocks)
  if size(y[1])[2] == 1
    # y = [y[i][:, 1] for i in 1:nBlocks]
    @threads for i in 1:nBlocks
      yNew[i] = y[i][:, 1]
    end
    y = yNew
    if length(wts) == 0
      # mu = [(y[i] .+ T(0.5))./2 for i in 1:nBlocks]
      @threads for i in 1:nBlocks
        mu[i] = (y[i] .+ T(0.5))./2
      end
    else
      # mu = [(wts[i] .* y .+ T(0.5))./(wts[i] .+ T(1)) for i in 1:nBlocks]
      @threads for i in 1:nBlocks
        mu[i] = (wts[i] .* y .+ T(0.5))./(wts[i] .+ T(1))
      end
    end
  elseif size(y[1])[2] == 2
    # events = [y[i][:, 1] for i in 1:nBlocks]
    @threads for i in 1:nBlocks
      events[i] = y[i][:, 1]
    end
    # N = [events[i] .+ y[i][:, 2] for i in 1:nBlocks]
    @threads for i in 1:nBlocks
      N[i] = events[i] .+ y[i][:, 2]
    end
    @threads for i in 1:nBlocks
      n = size(y[i])[1]
      tmp = zeros(T, n)
      for j in 1:n
        tmp[j] = if N[i][j] != 0; events[i][j]/N[i][j]; else T(0) end
      end
      yNew[i] = tmp
    end
    y = yNew
    # y = [if N[i] != 0; events[i]/N[i]; else T(0) end for i in 1:n]
    if length(wts) != 0
      # [wts[i] .*= N[i] for i in 1:nBlocks]
      @threads for i in 1:nBlocks
        wts[i] = wts[i] .* N[i]
      end
    else
      wts = N
    end
    # mu = [(N[i] .* y[i] .+ T(0.5))./(N[i] .+ T(1)) for i in 1:nBlocks]
    @threads for i in 1:nBlocks
      mu[i] = (N[i] .* y[i] .+ T(0.5))./(N[i] .+ T(1))
    end
  else
    error("There was a problem with the dimensions of y")
  end
  return y, mu, wts
end

# Variance Functions
@inline function variance(::BinomialDistribution, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return mu .* (1 .- mu)
end

function variance(distrib::AbstractDistribution, mu::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [variance(distrib, mu[i]) for i in 1:nBlocks]
end

function variance(::Block1DParallel, distrib::AbstractDistribution, mu::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = variance(distrib, mu[i])
  end
  return ret
end

# As in R
function y_log_y(y::T, x::T)::T where {T <: AbstractFloat}
  return y != 0 ? y * log(y/x) : 0
end
# Deviance and Chi-Squared for BernoulliDistribution/BinomialDistribution
@inline function devianceResiduals(::BinomialDistribution, mu::Array{T, 1},
            y::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return 2 .* (y_log_y.(y, mu) .+ y_log_y.(1 .- y, 1 .- mu))
end
@inline function devianceResiduals(::BinomialDistribution, mu::Array{T, 1},
            y::Array{T, 1}, wts::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return 2 .* wts .* (y_log_y.(y, mu) .+ y_log_y.(1 .- y, 1 .- mu))
end
@inline function X2(distrib::BinomialDistribution, mu::Array{T, 1}, y::Array{T, 1})::T where {T <: AbstractFloat}
  return sum(((y .- mu).^2)./variance(distrib, mu))
end

# For 1D block
function devianceResiduals(distrib::AbstractDistribution, mu::Array{Array{T, 1}, 1},
        y::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks = length(y)
  return [devianceResiduals(distrib, mu[i], y[i]) for i in 1:nBlocks]
end

function devianceResiduals(distrib::AbstractDistribution, mu::Array{Array{T, 1}, 1},
        y::Array{Array{T, 1}, 1}, wts::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks = length(y)
  return [devianceResiduals(distrib, mu[i], y[i], wts[i]) for i in 1:nBlocks]
end

# Parallel 1D Block
function devianceResiduals(::Block1DParallel, distrib::AbstractDistribution, mu::Array{Array{T, 1}, 1},
        y::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks = length(y)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = devianceResiduals(distrib, mu[i], y[i])
  end
  return ret
end
function devianceResiduals(::Block1DParallel, distrib::AbstractDistribution, mu::Array{Array{T, 1}, 1},
        y::Array{Array{T, 1}, 1}, wts::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks = length(y)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = devianceResiduals(distrib, mu[i], y[i], wts[i])
  end
  return ret
end

#=
  PoissonDistribution
=#

@inline function init!(::PoissonDistribution, y::Array{T}, wts::Array{T, 1})::InitType{T} where {T <: AbstractFloat}
  y = y[:, 1]
  mu = y .+ T(0.1)
  return y, mu, wts
end

function init!(::PoissonDistribution, y::Array{Array{T, 2}, 1},
            wts::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  tmp = T(0.1)
  y = [y[i][:, 1] for i in 1:nBlocks]
  mu = [y[i][:, 1] .+ tmp for i in 1:nBlocks]
  return y, mu, wts
end

# The Variance Function
@inline function variance(::PoissonDistribution, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return mu
end
# Deviance and Chi-Squared for PoissonDistribution
@inline function devianceResiduals(::PoissonDistribution, mu::Array{T, 1},
    y::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  dev::Array{T, 1} = zeros(T, size(y))
  for i in 1:length(y)
    if y[i] == 0
      dev[i] = 2 * mu[i]
    else y[i] > 0
      dev[i] = 2 * (y[i] * log(y[i]/mu[i]) - (y[i] - mu[i]))
    end
  end
  return dev
end
@inline function devianceResiduals(::PoissonDistribution, mu::Array{T, 1},
    y::Array{T, 1}, wts::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  dev::Array{T, 1} = zeros(T, size(y))
  for i in 1:length(y)
    if y[i] == 0
      dev[i] = 2 * wts[i] * mu[i]
    else y[i] > 0
      dev[i] = 2 * wts[i] * (y[i] * log(y[i]/mu[i]) - (y[i] - mu[i]))
    end
  end
  return dev
end
@inline function X2(distrib::PoissonDistribution, mu::Array{T, 1}, y::Array{T, 1})::T where {T <: AbstractFloat}
  return sum(((y .- mu).^2)./variance(distrib, mu))
end

#=
  GammaDistribution
=#
@inline function variance(::GammaDistribution, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return mu.^2
end
# Deviance and Chi-Squared for GammaDistribution
@inline function devianceResiduals(::GammaDistribution, mu::Array{T, 1},
    y::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return 2 .*(((y .- mu)./mu) .- log.(y./mu))
end
@inline function devianceResiduals(::GammaDistribution, mu::Array{T, 1},
    y::Array{T, 1}, wts::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return 2 .* wts .*(((y .- mu)./mu) .- log.(y./mu))
end
@inline function X2(distrib::GammaDistribution, mu::Array{T, 1}, y::Array{T, 1})::T where {T <: AbstractFloat}
  return sum(((y .- mu).^2)./variance(distrib, mu))
end

#=
  GaussianDistribution
=#
@inline function variance(::GaussianDistribution, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return fill(T(1), size(mu))
end
# Deviance and Chi-Squared for GaussianDistribution
@inline function devianceResiduals(distrib::GaussianDistribution, mu::Array{T, 1},
    y::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return (y .- mu).^2
end
@inline function devianceResiduals(distrib::GaussianDistribution, mu::Array{T, 1},
  y::Array{T, 1}, wts::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return wts .* (y .- mu).^2
end
@inline function X2(distrib::GaussianDistribution, mu::Array{T, 1}, y::Array{T, 1})::T where {T <: AbstractFloat}
  return sum(((y .- mu).^2)./variance(distrib, mu))
end

#=
  InverseGaussianDistribution
=#
@inline function variance(::InverseGaussianDistribution, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return mu.^3
end

# Deviance and Chi-Squared for InverseGaussianDistribution
@inline function devianceResiduals(::InverseGaussianDistribution, mu::Array{T, 1},
            y::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return ((y .- mu).^2)./(y .* (mu.^2))
end
@inline function devianceResiduals(::InverseGaussianDistribution, mu::Array{T, 1},
    y::Array{T, 1}, wts::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return wts .* ((y .- mu).^2)./(y .* (mu.^2))
end
@inline function X2(distrib::InverseGaussianDistribution, mu::Array{T, 1}, y::Array{T, 1})::T where {T <: AbstractFloat}
  sqDiff = (y .- mu).^2
  v = sqrt(sum(sqDiff)/length(y))
  return v*sum((sqDiff)./variance(distrib, mu))
end

#=
  NegativeBernoulliDistribution
=#
@inline function init!(::NegativeBernoulliDistribution{T}, y::Array{T}, wts::Array{T, 1})::InitType{T} where {T <: AbstractFloat}
  y = y[:, 1]
  mu = copy(y)
  tmp = T(1/6)
  for i in 1:length(y)
    if y[i] == 0
      mu[i] += tmp
    end
  end
  return y, mu, wts
end

function init!(::NegativeBernoulliDistribution{T}, y::Array{Array{T, 2}, 1},
            wts::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  tmp = T(1/6)
  y = [y[i][:, 1] for i in 1:nBlocks]
  mu = [y[i][:, 1] .+ tmp for i in 1:nBlocks]
  return y, mu, wts
end

# Variance functions
@inline function variance(distrib::NegativeBernoulliDistribution{T}, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return mu .+ distrib.Alpha .* mu.^2
end

# Deviance and Chi-Squared for NegativeBernoulliDistribution
@inline function devianceResiduals(distrib::NegativeBernoulliDistribution{T}, mu::Array{T, 1},
            y::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  Alpha = distrib.Alpha
  dev::Array{T, 1} = zeros(T, size(y))
  iAlpha = T(1)/Alpha
  for i in 1:length(y)
    if y[i] == 0
      dev[i] = 2 * iAlpha*log(1/(1 + Alpha*mu[i]))
    else y[i] > 0
      dev[i] = 2 * (y[i] * log(y[i]/mu[i]) - (y[i] + iAlpha)*log((1 + Alpha*y[i])/(1 + Alpha*mu[i])))
    end
  end
  return dev
end
@inline function devianceResiduals(distrib::NegativeBernoulliDistribution{T}, mu::Array{T, 1},
            y::Array{T, 1}, wts::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  Alpha = distrib.Alpha
  dev::Array{T, 1} = zeros(T, size(y))
  iAlpha = T(1)/Alpha
  for i in 1:length(y)
    if y[i] == 0
      dev[i] = 2 * wts[i] * iAlpha*log(1/(1 + Alpha*mu[i]))
    else y[i] > 0
      dev[i] = 2 * wts[i] * (y[i] * log(y[i]/mu[i]) - (y[i] + iAlpha)*log((1 + Alpha*y[i])/(1 + Alpha*mu[i])))
    end
  end
  return dev
end
@inline function X2(distrib::NegativeBernoulliDistribution{T}, mu::Array{T, 1}, y::Array{T, 1})::T where {T <: AbstractFloat}
  return sum(((y .- mu).^2)./variance(distrib, mu))
end

#=
  PowerDistribution
=#

# Variance distribution functions
@inline function variance(distrib::PowerDistribution{T}, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  k = distrib.k
  return mu .^k
end

# Deviance and Chi-Squared for PowerDistribution
@inline function devianceResiduals(distrib::PowerDistribution{T}, mu::Array{T, 1},
            y::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  k = distrib.k
  dev::Array{T, 1} = zeros(T, size(y))
  ok = 1 - k
  tk = 2 - k
  for i in 1:length(y)
    dev[i] = ((2*y[i]/(ok*((y[i]^(ok)) - (mu[i]^ok)))) - (2/(tk * ((y[i]^(tk)) - (mu[i]^tk)))))
  end
  return dev
end
@inline function devianceResiduals(distrib::PowerDistribution{T}, mu::Array{T, 1},
            y::Array{T, 1}, wts::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  k = distrib.k
  dev::Array{T, 1} = zeros(T, size(y))
  ok = 1 - k
  tk = 2 - k
  for i in 1:length(y)
    dev[i] = wts[i] * ((2*y[i]/(ok*((y[i]^(ok)) - (mu[i]^ok)))) - (2/(tk * ((y[i]^(tk)) - (mu[i]^tk)))))
  end
  return dev
end
@inline function deviance(distrib::AbstractDistribution, mu::Array{T, 1}, y::Array{T, 1})::T where {T <: AbstractFloat}
  return sum(devianceResiduals(distrib, mu, y))
end

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
function copy_sp_data_range!(newR, newstrd, newoffsets, ind, R, rg)
  for i ∈ rg
    push!(newR.args, R[i])
    push!(newstrd.args, :($getfield(strd, $i, false)))
    push!(newoffsets.args, :($getfield(offs, $i, false)))
    push!(ind.args, :($getfield(offs, $i, false)))
  end
  nothing
end
@generated function subsetview(
  ptr::AbstractStridedPointer{T,N,C,B,R,X,O},
  ::StaticInt{I},
  i::Integer,
) where {T,N,C,B,R,X,O,I}
  I > N && return :ptr
  @assert B ≤ 0 "Batched dims not currently supported."
  newC = C == I ? -1 : (C < I ? C : C - 1)
  newR = Expr(:tuple)
  newstrd = Expr(:tuple)
  newoffsets = Expr(:tuple)
  ind = Expr(:tuple)
  copy_sp_data_range!(newR, newstrd, newoffsets, ind, R, 1:I-1)
  push!(ind.args, :i)
  copy_sp_data_range!(newR, newstrd, newoffsets, ind, R, I+1:N)
  gptr = Expr(:call, :gep, :ptr, ind)
  quote
    $(Expr(:meta, :inline))
    strd = strides(ptr)
    offs = offsets(ptr)
    si = ArrayInterface.StrideIndex{$(N - 1),$newR,$newC}($newstrd, $newoffsets)
    stridedpointer($gptr, si, StaticInt{$B}())
  end
end
@inline _subsetview(ptr::AbstractStridedPointer, ::StaticInt{I}, J::Tuple{}) where {I} = ptr
@inline _subsetview(
  ptr::AbstractStridedPointer,
  ::StaticInt{I},
  J::Tuple{J1},
) where {I,J1} = subsetview(ptr, StaticInt{I}(), first(J))
@inline _subsetview(
  ptr::AbstractStridedPointer,
  ::StaticInt{I},
  J::Tuple{J1,J2,Vararg},
) where {I,J1,J2} =
  _subsetview(subsetview(ptr, StaticInt{I}(), first(J)), StaticInt{I}(), Base.tail(J))
@inline subsetview(
  ptr::AbstractStridedPointer,
  ::StaticInt{I},
  J::CartesianIndex,
) where {I} = _subsetview(ptr, StaticInt{I}(), Tuple(J))

@inline _gesp(sp::VectorizationBase.FastRange, ::StaticInt{1}, i, ::StaticInt{1}) =
  gesp(sp, (i,))
@generated function _gesp(
  sp::AbstractStridedPointer{T,N},
  ::StaticInt{I},
  i::Integer,
  ::StaticInt{D},
) where {I,N,T,D}
  t = Expr(:tuple)
  for j ∈ 1:I-1
    # push!(t.args, staticexpr(0))
    push!(t.args, VectorizationBase.NullStep())
  end
  push!(t.args, :i)
  if D > 1
    for j ∈ I+1:N
      # push!(t.args, staticexpr(0))
      push!(t.args, VectorizationBase.NullStep())
    end
  end
  quote
    $(Expr(:meta, :inline))
    gesp(sp, $t)
  end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Boundary{P,Q,T}

The boundary relation from rank `P` to smaller rank `Q` for
topology of type `T`.
"""
struct Boundary{P,Q,T<:Topology} <: TopologicalRelation
  topology::T
end

Boundary{P,Q}(topology::T) where {P,Q,T} = Boundary{P,Q,T}(topology)

# -------------------
# HALF-EDGE TOPOLOGY
# -------------------

function (∂::Boundary{2,1,T})(elem::Integer) where {T<:HalfEdgeTopology}
  t = ∂.topology
  l = loop(half4elem(elem, t))
  v = CircularVector(l)
  [edge4pair((v[i], v[i+1]), t) for i in 1:length(v)]
end

function (∂::Boundary{2,0,T})(elem::Integer) where {T<:HalfEdgeTopology}
  loop(half4elem(elem, ∂.topology))
end

function (∂::Boundary{1,0,T})(edge::Integer) where {T<:HalfEdgeTopology}
  e = half4edge(edge, ∂.topology)
  [e.head, e.half.head]
end


text:
# generate the dyadic tree, suitable for output to DOT
using SurrealNumbers
# var_args = "color=red"
var_args = ""
n=3; x_scale = 2.25
# n=4; x_scale = 5.00
out_dir = "Data/"

######## FIRST PICTURE

file = "$(out_dir)dyadic_tree_$n.dot"
io = open(file, "w")
println(io, "digraph dyadic_tree {")

name_map = Dict{SurrealFinite, Integer}()

z = zero(SurrealFinite)
m = 0
name_map[z] = m
var_args = "pos=\"0,0!\""
surreal2node(io, z, 0, extra_args=var_args)
println(io, "node_day_0 [shape=none, pos=\"$(x_scale*(n+0.5)),$(0)!\", label=\"Day 0\", fontsize=18]")

for k=1:n
    println(io, "node_day_$k [shape=none, pos=\"$(x_scale*(n+0.5)),$(k)!\", label=\"Day $k\", fontsize=18]")

    var_args = "pos=\"-$(x_scale*k),$(k)!\""
    x = convert(SurrealFinite, -k)
    m += 1
    surreal2node(io, x, m, extra_args=var_args)
    name_map[x] = m
    println(" k = $k, x = $x")
    tmp = name_map[x.R[1]]
    println(io, "   node_$m:R -> node_$tmp [color=\"blue3\"];")

    var_args = "pos=\"$(x_scale*k),$(k)!\""
    x = convert(SurrealFinite, k)
    m += 1
    surreal2node(io, x, m, extra_args=var_args)
    name_map[x] = m
    println(" k = $k, x = $x")
    tmp = name_map[x.L[1]]
    println(io, "   node_$m:L -> node_$tmp [color=\"red3\"];")

    for f=2:k
        for j=(f-2)*2^(k-f+1) + (1:2:2^(k-f+1))
            x = convert(SurrealFinite, j // 2^(k-f+1))
            var_args = "pos=\"$(x_scale*convert(Float64, x)),$(k)!\""
            m += 1
            surreal2node(io, x, m, extra_args=var_args)
            name_map[x] = m
            println(" k = $k, f = $f, x = $x")
            tmp = name_map[x.L[1]]
            println(io, "   node_$m:L -> node_$tmp [color=\"red3\"];")
            tmp = name_map[x.R[1]]
            println(io, "   node_$m:R -> node_$tmp [color=\"blue3\"];")

            x = convert(SurrealFinite, -j // 2^(k-f+1))
            var_args = "pos=\"$(x_scale*convert(Float64, x)),$(k)!\""
            m += 1
            surreal2node(io, x, m, extra_args=var_args)
            name_map[x] = m
            println(" k = $k, f = $f, x = $x")
            tmp = name_map[x.L[1]]
            println(io, "   node_$m:L -> node_$tmp [color=\"red3\"];")
            tmp = name_map[x.R[1]]
            println(io, "   node_$m:R -> node_$tmp [color=\"blue3\"];")
        end
    end
end

println(io, "}")
close(io)
# run(`dot -n -Kfdp -Tpdf -O $file`)
run(`neato -Tpdf -O $file`)
run(`neato -Tsvg -O $file`)

######## Colouring

n=4; x_scale = 5.00
# n=3; x_scale = 2.25

file2 = "$(out_dir)dyadic_tree_chi_$n.dot"
io = open(file2, "w")
println(io, "digraph dyadic_tree {")

name_map2 = Dict{SurrealFinite, Integer}()
color = Dict{SurrealFinite, Integer}()
color_map = Dict(1=>"red3", 2=>"blue3", 3=>"green3")
C = Dict(1=>"R", 2=>"B", 3=>"G")
fs = 24

z = zero(SurrealFinite)
m = 0
name_map2[z] = m
var_args = "pos=\"0,0!\""
cn = 1
println(io, "   node_$m [label=0, fontsize=$fs, color=\"$(color_map[cn])\", $var_args];")
x = z
color[x] = cn

for k=1:n
    cn = mod(k,2) + 1

    var_args = "pos=\"-$(x_scale*k),$(k)!\""
    x = convert(SurrealFinite, -k)
    m += 1
    println(io, "   node_$m [label=$(-k), fontsize=$fs, color=\"$(color_map[cn])\", $var_args];")
    name_map2[x] = m
    color[x] = cn
    println(" k = $k, x = $x")
    tmp = name_map2[x.R[1]]
    println(io, "   node_$m -> node_$tmp;")

    var_args = "pos=\"$(x_scale*k),$(k)!\""
    x = convert(SurrealFinite, k)
    m += 1
    println(io, "   node_$m [label=$(k), fontsize=$fs, color=\"$(color_map[cn])\", $var_args];")
    name_map2[x] = m
    color[x] = cn
    println(" k = $k, x = $x")
    tmp = name_map2[x.L[1]]
    println(io, "   node_$m -> node_$tmp;")

    for f=2:k
        for j=(f-2)*2^(k-f+1) + (1:2:2^(k-f+1))
            v = j // 2^(k-f+1)
            s = string(v.num) * "/" * string(v.den)
            x = convert(SurrealFinite, v)
            cn = setdiff([1,2,3], [color[x.L[1]], color[x.R[1]]] )[1]
            color[x] = cn

            var_args = "pos=\"$(x_scale*convert(Float64, x)),$(k)!\""
            m += 1
            println(io, "   node_$m [label=\"$(s)\", fontsize=$fs, color=\"$(color_map[cn])\", $var_args];")
            name_map[x] = m
            println(" k = $k, f = $f, x = $x")
            tmp = name_map[x.L[1]]
            println(io, "   node_$m -> node_$tmp;")
            tmp = name_map[x.R[1]]
            println(io, "   node_$m -> node_$tmp;")

            v = -j // 2^(k-f+1)
            s = string(v.num) * "/" * string(v.den)
            x = convert(SurrealFinite, v)
            cn = setdiff([1,2,3], [color[x.L[1]], color[x.R[1]]] )[1]
            color[x] = cn
            var_args = "pos=\"$(x_scale*convert(Float64, x)),$(k)!\""
            m += 1
            println(io, "   node_$m [label=\"$(s)\", fontsize=$fs, color=\"$(color_map[cn])\", $var_args];")
            name_map[x] = m
            println(" k = $k, f = $f, x = $x")
            tmp = name_map[x.L[1]]
            println(io, "   node_$m -> node_$tmp;")
            tmp = name_map[x.R[1]]
            println(io, "   node_$m -> node_$tmp;")
        end
    end
end

println(io, "}")
close(io)
# run(`dot -n -Kfdp -Tpdf -O $file`)
run(`neato -Tpdf -O $file2`)
run(`neato -Tsvg -O $file2`)

X = sort( collect(keys(color)))
for x in X
    print(" , $(C[color[x]])")
end
println()


text:
using PyCall

@pyimport selenium.webdriver as wd

immutable Driver
    o::PyObject
end
PyObject(x::Driver) = x.o

for attr in [:session_id, :capabilities, :command_executor, :error_handler]
    @eval $attr(w::Driver) = w.o[:$attr]
end

immutable WebElement
    o::PyObject
end
PyObject(x::WebElement) = x.o

export init_phantomjs, init_firefox, init_chrome, init_ie

init_phantomjs(args...) = Driver(wd.PhantomJS(args...))
init_firefox(args...) = Driver(wd.Firefox(args...))
init_chrome(args...) = Driver(wd.Chrome(args...))
init_ie(args...) = Driver(wd.Ie(args...))

immutable ActionChain
    o::PyObject
end

ActionChain(x::Driver) = ActionChain(wd.ActionChains(x.o))

PyObject(x::ActionChain) = x.o


text:
@theme GitHubTheme Dict(
    :name => "GitHub",
    :description => "Based on the Rouge 'GitHub' theme.",
    :style => S"bg: f8f8f8",
    :tokens => Dict(
        COMMENT             => S"fg: 999988; italic",
        COMMENT_MULTILINE   => S"fg: 999988; italic",
        COMMENT_PREPROC     => S"fg: 999999; bold",
        COMMENT_SINGLE      => S"fg: 999988; italic",
        COMMENT_SPECIAL     => S"fg: 999999; italic; bold",

        ERROR               => S"fg: a61717; bg: e3d2d2",

        GENERIC             => S"fg: 000000",
        GENERIC_DELETED     => S"fg: 000000; bg: ffdddd",
        GENERIC_EMPH        => S"fg: 000000; italic",
        GENERIC_ERROR       => S"fg: aa0000",
        GENERIC_HEADING     => S"fg: 999999",
        GENERIC_INSERTED    => S"fg: 000000; bg: ddffdd",
        GENERIC_OUTPUT      => S"fg: 888888",
        GENERIC_PROMPT      => S"fg: 555555",
        GENERIC_STRONG      => S"bold",
        GENERIC_SUBHEADING  => S"fg: aaaaaa",
        GENERIC_TRACEBACK   => S"fg: aa0000",

        KEYWORD             => S"fg: 000000; bold",
        KEYWORD_CONSTANT    => S"fg: 000000; bold",
        KEYWORD_DECLARATION => S"fg: 000000; bold",
        KEYWORD_NAMESPACE   => S"fg: 000000; bold",
        KEYWORD_PSEUDO      => S"fg: 000000; bold",
        KEYWORD_RESERVED    => S"fg: 000000; bold",
        KEYWORD_TYPE        => S"fg: 445588; bold",

        NUMBER              => S"fg: 009999",
        NUMBER_FLOAT        => S"fg: 009999",
        NUMBER_HEX          => S"fg: 009999",
        NUMBER_INTEGER      => S"fg: 009999",
        NUMBER_INTEGER_LONG => S"fg: 009999",
        NUMBER_OCT          => S"fg: 009999",

        STRING              => S"fg: d14",
        STRING_BACKTICK     => S"fg: d14",
        STRING_CHAR         => S"fg: d14",
        STRING_DOC          => S"fg: d14",
        STRING_DOUBLE       => S"fg: d14",
        STRING_ESCAPE       => S"fg: d14",
        STRING_HEREDOC      => S"fg: d14",
        STRING_INTERPOL     => S"fg: d14",
        STRING_OTHER        => S"fg: d14",
        STRING_REGEX        => S"fg: 009926",
        STRING_SINGLE       => S"fg: d14",
        STRING_SYMBOL       => S"fg: 990073",

        NAME                => S"fg: 008080",
        NAME_ATTRIBUTE      => S"fg: 008080",
        NAME_BUILTIN        => S"fg: 0086b3",
        NAME_BUILTIN_PSEUDO => S"fg: 999999",
        NAME_CLASS          => S"fg: 445588; bold",
        NAME_CONSTANT       => S"fg: 008080",
        NAME_DECORATOR      => S"fg: 3c5d5d; bold",
        NAME_ENTITY         => S"fg: 800080",
        NAME_EXCEPTION      => S"fg: 990000; bold",
        NAME_FUNCTION       => S"fg: 990000; bold",
        NAME_LABEL          => S"fg: 990000; bold",
        NAME_NAMESPACE      => S"fg: 555555",
        NAME_TAG            => S"fg: 000080",

        NAME_VARIABLE          => S"fg: 008080",
        NAME_VARIABLE_CLASS    => S"fg: 008080",
        NAME_VARIABLE_GLOBAL   => S"fg: 008080",
        NAME_VARIABLE_INSTANCE => S"fg: 008080",

        OPERATOR            => S"fg: 000000; bold",
        OPERATOR_WORD       => S"fg: 000000; bold",

        TEXT                => S"bg: f8f8f8",
        WHITESPACE     => S"fg: bbbbbb",
    )
)


text:
struct Linear <: Degree{1} end  # boundary conditions not supported

"""
    Linear()

Indicate that the corresponding axis should use linear interpolation.

# Extended help

Assuming uniform knots with spacing 1, the `i`th piece of linear b-spline
implemented here is defined as follows.

    y_i(x) = c p(x) + cp p(1-x)

where

    p(δx) = x

and the values `cX` for `X ∈ {_, p}` are the coefficients.

Linear b-splines are naturally interpolating, and require no prefiltering;
there is therefore no need for boundary conditions to be provided.

Also, although the implementation is slightly different in order to re-use
the framework built for general b-splines, the resulting interpolant is just
a piecewise linear function connecting each pair of neighboring data points.
"""
Linear

function positions(::Linear, ax::AbstractUnitRange{<:Integer}, x)
    f = floor(x)
    # When x == last(ax) we want to use the x-1, x pair
    f = ifelse(x == last(ax), f - oneunit(f), f)
    fi = fast_trunc(Int, f)
    return fi, x-f
end

value_weights(::Linear, δx) = (1-δx, δx)
gradient_weights(::Linear, δx) = (-oneunit(δx), oneunit(δx))
hessian_weights(::Linear, δx) = (zero(δx), zero(δx))

padded_axis(ax::AbstractUnitRange, ::BSpline{Linear}) = ax


text:
using LinearAlgebra
using ProgressMeter

struct HamiltonSystem
    d
    m
    F
    q0
    p0
    Minv
    sol
    U
    HamiltonSystem(d,m,F,q0,p0,Minv,sol,U) =
        length(q0)==d*length(m) && length(p0)==d*length(m) ?
        new(d,m,F,q0,p0,Minv,sol,U) : error("dimension mismatch")
end

function getMinv(m,d)
    m_long=zeros(d*length(m))
    for (i,n) in enumerate(m)
        m_long[(i-1)*d+1:(i)*d] = length(n)==d ? collect(n) : [n for i in 1:d]
    end
    return Diagonal(1 ./ m_long)
end

HamiltonSystem(d,m,q0,p0,F=missing,U=missing,sol=missing) = HamiltonSystem(d,m,F,q0,p0,getMinv(m,d),sol,U)

function getEnergy(HS,q,p)
    return 1/2*p'*HS.Minv*p + HS.U(q)
end

function getEnergyNorm(HS,Q,P)
    QA = [Q[:,i] for i in 1:size(Q,2)]
    PA = [P[:,i] for i in 1:size(P,2)]
    E = map(getEnergy,[HS for q in QA],QA,PA)
    E = E.-E[1]
    return map(abs,E)
end

getEnergyErrors(HS,Qs,Ps) = [getEnergyNorm(HS,Q,P) for (Q,P) in zip(Qs,Ps)]
#getMaximumErrors(HS,ts,Qs)= [[norm((Q-hcat(map(HS.sol,t)...))[:,i]) for i in 1:length(t)] for (t,Q) in zip(ts,Qs)]
function getMaximumErrors(HS,ts,Qs)
    ME = []
    for (t,Q) in zip(ts,Qs)
        QA = QA = [Q[:,i] for i in 1:size(Q,2)]
        E=map(norm,QA-map(HS.sol,t))
        push!(ME,E)
    end
    return ME
end

println("Finished loading HamiltonSystem")

max_stars_repo_name: BoergeSt/HamiltonODE

max_issues_repo_name: BoergeSt/HamiltonODE

max_forks_repo_name: BoergeSt/HamiltonODE


text:
using Adapt
using Statistics
using Oceananigans.Grids
using Oceananigans.Grids: interior_parent_indices

struct AveragedField{X, Y, Z, S, A, D, G, T, N, O} <: AbstractReducedField{X, Y, Z, A, G, T, N}
            data :: D
    architecture :: A
            grid :: G
            dims :: NTuple{N, Int}
         operand :: O
          status :: S

    function AveragedField{X, Y, Z}(data::D, arch::A, grid::G, dims, operand::O;
                                    recompute_safely=true) where {X, Y, Z, D, A, G, O}

        dims = validate_reduced_dims(dims)
        validate_reduced_locations(X, Y, Z, dims)
        validate_field_data(X, Y, Z, data, grid)

        status = recompute_safely ? nothing : FieldStatus(0.0)

        S = typeof(status)
        N = length(dims)
        T = eltype(grid)

        return new{X, Y, Z, S, A, D, G, T, N, O}(data, arch, grid, dims, operand, status)
    end

    function AveragedField{X, Y, Z}(data::D, arch::A, grid::G, dims, operand::O, status::S) where {X, Y, Z, D, A, G, O, S}
        return new{X, Y, Z, S, A, D, G, eltype(grid), length(dims), O}(data, arch, grid, dims, operand, status)
    end
end

"""
    AveragedField(operand::AbstractField; dims, data=nothing, recompute_safely=false)

Returns an AveragedField.
"""
function AveragedField(operand::AbstractField; dims, data=nothing, recompute_safely=true)

    arch = architecture(operand)
    loc = reduced_location(location(operand), dims=dims)
    grid = operand.grid

    if isnothing(data)
        data = new_data(arch, grid, loc)
        recompute_safely = false
    end

    return AveragedField{loc[1], loc[2], loc[3]}(data, arch, grid, dims, operand,
                                                 recompute_safely=recompute_safely)
end

"""
    compute!(avg::AveragedField, time=nothing)

Compute the average of `avg.operand` and store the result in `avg.data`.
"""
function compute!(avg::AveragedField, time=nothing)
    compute_at!(avg.operand, time)
    avg .= 0 # should we have to do this?
    mean!(avg, avg.operand)
    return nothing
end

compute_at!(avg::AveragedField{X, Y, Z, <:FieldStatus}, time) where {X, Y, Z} =
    conditional_compute!(avg, time)

#####
##### Adapt
#####

Adapt.adapt_structure(to, averaged_field::AveragedField{X, Y, Z}) where {X, Y, Z} =
    AveragedField{X, Y, Z}(Adapt.adapt(to, averaged_field.data), nothing,
                           nothing, averaged_field.dims, nothing, nothing)