text:
__precompile__(true)
module DataStreams

module Data

using Missings, WeakRefStrings

import Core.Compiler: return_type

# Data.Schema
"""
A `Data.Schema` describes a tabular dataset, i.e. a set of named, typed columns with records as rows

`Data.Schema` allow `Data.Source` and `Data.Sink` to talk to each other and prepare to provide/receive data through streaming.
`Data.Schema` provides the following accessible properties:

 * `Data.header(schema)` to return the header/column names in a `Data.Schema`
 * `Data.types(schema)` to return the column types in a `Data.Schema`; `Union{T, Missing}` indicates columns that may contain missing data (`missing` values)
 * `Data.size(schema)` to return the (# of rows, # of columns) in a `Data.Schema`; note that # of rows may be `missing`, meaning unknown

`Data.Schema` has the following constructors:

 * `Data.Schema()`: create an "empty" schema with no rows, no columns, and no column names
 * `Data.Schema(types[, header, rows, meta::Dict])`: column element types are provided as a tuple or vector; column names provided as an iterable; # of rows can be an Int or `missing` to indicate unknown # of rows

`Data.Schema` are indexable via column names to get the number of that column in the `Data.Schema`

```julia
julia> sch = Data.Schema([Int], ["column1"], 10)
Data.Schema:
rows: 10	cols: 1
Columns:
 "column1"  Int64

julia> sch["column1"]
1
```

**Developer note**: the full type definition is `Data.Schema{R, T}` where the `R` type parameter will be `true` or `false`, indicating
whether the # of rows are known (i.e not `missing`), respectively. The `T` type parameter is a `Tuple{A, B, ...}` representing the column element types
in the `Data.Schema`. Both of these type parameters provide valuable information that may be useful when constructing `Sink`s or streaming data.
"""
mutable struct Schema{R, T}
    # types::T               # Julia types of columns
    header::Vector{String}   # column names
    rows::Union{Int, Missing}# number of rows in the dataset
    cols::Int                # number of columns in a dataset
    metadata::Dict           # for any other metadata we'd like to keep around (not used for '==' operation)
    index::Dict{String, Int} # maps column names as Strings to their index # in `header` and `types`
end

function Schema(types=(), header=["Column$i" for i = 1:length(types)], rows::Union{Integer,Missing}=0, metadata::Dict=Dict())
    !ismissing(rows) && rows < 0 && throw(ArgumentError("Invalid # of rows for Data.Schema; use `missing` to indicate an unknown # of rows"))
    types2 = Tuple(types)
    header2 = String[string(x) for x in header]
    cols = length(header2)
    cols != length(types2) && throw(ArgumentError("length(header): $(length(header2)) must == length(types): $(length(types2))"))
    return Schema{!ismissing(rows), Tuple{types2...}}(header2, rows, cols, metadata, Dict(n=>i for (i, n) in enumerate(header2)))
end
Schema(types, rows::Union{Integer,Missing}, metadata::Dict=Dict()) = Schema(types, ["Column$i" for i = 1:length(types)], rows, metadata)

header(sch::Schema) = sch.header
types(sch::Schema{R, T}) where {R, T} = Tuple(T.parameters)

function anytypes(sch::Schema{R, T}, weakref) where {R, T}
    types = T.parameters
    if !weakref
        types = map(x->x >: Missing ? ifelse(Missings.T(x) <: WeakRefString, Union{String, Missing}, x) : ifelse(x <: WeakRefString, String, x), types)
    end
    return collect(Any, types)
end

metadata(sch::Schema) = sch.metadata
Base.size(sch::Schema) = (sch.rows, sch.cols)
Base.size(sch::Schema, i::Int) = ifelse(i == 1, sch.rows, ifelse(i == 2, sch.cols, 0))

setrows!(source, rows) = isdefined(source, :schema) ? (source.schema.rows = rows; nothing) : nothing
setrows!(source::Array, rows) = nothing
Base.getindex(sch::Schema, col::String) = sch.index[col]

function Base.show(io::IO, schema::Schema)
    println(io, "Data.Schema:")
    println(io, "rows: $(schema.rows)  cols: $(schema.cols)")
    if schema.cols > 0
        println(io, "Columns:")
        Base.print_matrix(io, hcat(schema.header, collect(types(schema))))
    end
end

function transform(sch::Data.Schema{R, T}, transforms::Dict{Int, <:Base.Callable}, weakref) where {R, T}
    types = Data.types(sch)
    transforms2 = ((get(transforms, x, identity) for x = 1:length(types))...,)
    # NOTE: `return_type` is imported, see above
    newtypes = ((return_type(transforms2[x], (types[x],)) for x = 1:length(types))...,)
    if !weakref
        newtypes = map(x->x >: Missing ? ifelse(Missings.T(x) <: WeakRefString, Union{String, Missing}, x) : ifelse(x <: WeakRefString, String, x), newtypes)
    end
    return Schema(newtypes, Data.header(sch), size(sch, 1), sch.metadata), transforms2
end
transform(sch::Data.Schema, transforms::Dict{String, F}, s) where {F<:Base.Callable} =
    transform(sch, Dict{Int, F}(sch[x]=>f for (x, f) in transforms), s)

# Data.StreamTypes
abstract type StreamType end
struct Field  <: StreamType end
struct Row    <: StreamType end
struct Column <: StreamType end

# Data.Source Interface
abstract type Source end

# Required methods
"""
`Data.schema(s::Source) => Data.Schema`

Return the `Data.Schema` of a source, which describes the # of rows & columns, as well as the column types of a dataset.
Some sources like `CSV.Source` or `SQLite.Source` store the `Data.Schema` directly in the type, whereas
others like `DataFrame` compute the schema on the fly.

The `Data.Schema` of a source is used in various ways during the streaming process:

- The # of rows and if they are known are used to generate the inner streaming loop
- The # of columns determine if the innermost streaming loop can be unrolled automatically or not
- The types of the columns are used in loop unrolling to generate efficient and type-stable streaming

See `?Data.Schema` for more details on how to work with the type.
"""
function schema end
"""
`Data.isdone(source, row, col) => Bool`

Checks whether a source can stream additional fields/columns for a desired
`row` and `col` intersection. Used during the streaming process, especially for sources
that have an unknown # of rows, to detect when a source has been exhausted of data.

Data.Source types must at least implement:

`Data.isdone(source::S, row::Int, col::Int)`

If more convenient/performant, they can also implement:

`Data.isdone(source::S, row::Int, col::Int, rows::Union{Int, Missing}, cols::Int)`

where `rows` and `cols` are the size of the `source`'s schema when streaming.

A simple example of how a DataFrame implements this is:
```julia
Data.isdone(df::DataFrame, row, col, rows, cols) = row > rows || col > cols
```
"""
function isdone end
isdone(source, row, col, rows, cols) = isdone(source, row, col)

"""
`Data.streamtype{T<:Data.Source, S<:Data.StreamType}(::Type{T}, ::Type{S}) => Bool`

Indicates whether the source `T` supports streaming of type `S`. To be overloaded by individual sources according to supported `Data.StreamType`s.
This is used in the streaming process to determine the compatability of streaming from a specific source to a specific sink.
It also helps in determining the preferred streaming method, when matched up with the results of `Data.streamtypes(s::Sink)`.

For example, if `MyPkg.Source` supported `Data.Field` streaming, I would define:

```julia
Data.streamtype(::Type{MyPkg.Source}, ::Type{Data.Field}) = true
```

and/or for `Data.Column` streaming:

```julia
Data.streamtype(::Type{MyPkg.Source}, ::Type{Data.Column}) = true
```
"""
function streamtype end

"""
`Data.reset!(source)`

Resets a source into a state that allows streaming its data again.
For example, for `CSV.Source`, the internal buffer is "seek"ed back to the start position of
the csv data (after the column name headers). For `SQLite.Source`, the source SQL query is re-executed.
"""
function reset! end

"""
Data.Source types must implement one of the following:

`Data.streamfrom(source, ::Type{Data.Field}, ::Type{T}, row::Int, col::Int) where {T}`

`Data.streamfrom(source, ::Type{Data.Column}, ::Type{T}, col::Int) where {T}`

Performs the actually streaming of data "out" of a data source. For `Data.Field` streaming, the single field value of type `T`
at the intersection of `row` and `col` is returned. For `Data.Column` streaming, the column # `col` with element type `T` is returned.

For `Data.Column`, a source can also implement:
```julia
Data.streamfrom(source, ::Type{Data.Field}, ::Type{T}, row::Int, col::Int) where {T}
```
where `row` indicates the # of rows that have already been streamed from the source.
"""
function streamfrom end
Data.streamfrom(source, ::Type{Data.Column}, T, row, col) = Data.streamfrom(source, Data.Column, T, col)
Data.streamfrom(source, ::Type{Data.Column}, T, r::AbstractRange, col) = Data.streamfrom(source, Data.Column, T, first(r), col)[r]

# Generic fallbacks
Data.streamtype(source, ::Type{Data.Row}) = Data.streamtype(source, Data.Field)
Data.streamtype(source, ::Type{<:StreamType}) = false
Data.reset!(source) = nothing

struct RandomAccess end
struct Sequential end

"""
`Data.accesspattern(source) => Data.RandomAccess | Data.Sequential`

returns the data access pattern for a Data.Source.

`RandomAccess` indicates that a source supports streaming data (via calls to `Data.streamfrom`) with arbitrary row/column values in any particular order.

`Sequential` indicates that the source only supports streaming data sequentially, starting w/ row 1, then accessing each column from 1:N, then row 2, and each column from 1:N again, etc.

For example, a `DataFrame` holds all data in-memory, and thus supports easy random access in any order.
A `CSV.Source` however, is streaming the contents of a file, where rows must be read sequentially, and each column sequentially within each rows.

By default, sources are assumed to have a `Sequential` access pattern.
"""
function accesspattern end

accesspattern(x) = Sequential()

const EMPTY_REFERENCE = UInt8[]
"""
Data.Source types can optionally implement

`Data.reference(x::Source) => Vector{UInt8}`

where the type retruns a `Vector{UInt8}` that represents a memory block that should be kept in reference for WeakRefStringArrays.

In many streaming situations, the minimizing of data copying/movement is ideal. Some sources can provide in-memory access to their data
in the form of a `Vector{UInt8}`, i.e. a single byte vector, that sinks can "point to" when streaming, instead of needing to always
copy all the data. In particular, the `WeakRefStrings` package provides utilities for creating "string types" that don't actually hold their
own data, but instead just "point" to data that lives elsewhere. As Strings can be some of the most expensive data structures to copy
and move around, this provides excellent performance gains in some cases when the sink is able to leverage this alternative structure.
"""
function reference end
reference(x) = EMPTY_REFERENCE

# Data.Sink Interface
"""
Represents a type that can have data streamed to it from `Data.Source`s.

To satisfy the `Data.Sink` interface, it must provide two constructors with the following signatures:

```
[Sink](sch::Data.Schema, S::Type{StreamType}, append::Bool, args...; reference::Vector{UInt8}=UInt8[], kwargs...) => [Sink]
[Sink](sink, sch::Data.Schema, S::Type{StreamType}, append::Bool; reference::Vector{UInt8}=UInt8[]) => [Sink]
```

Let's break these down, piece by piece:

- `[Sink]`: this is your sink type, i.e. `CSV.Sink`, `DataFrame`, etc. You're defining a constructor for your sink type.
- `sch::Data.Schema`: in the streaming process, the schema of a `Data.Source` is provided to the sink in order to allow the sink to "initialize" properly in order to receive data according to the format in `sch`. This might mean pre-allocating space according to the # of rows/columns in the source, managing the sink's own schema to match `sch`, etc.
- `S::Type{StreamType}`: `S` represents the type of streaming that will occur from the `Data.Source`, either `Data.Field` or `Data.Column`
- `append::Bool`: a boolean flag indicating whether the data should be appended to a sink's existing data store, or, if `false`, if the sink's data should be fully replaced by the incoming `Data.Source`'s data
- `args...`: In the 1st constructor form, `args...` represents a catchall for any additional arguments your sink may need to construct. For example, `SQLite.jl` defines `Sink(sch, S, append, db, table_name)`, meaning that the `db` and `table_name` are additional required arguments in order to properly create an `SQLite.Sink`.
- `reference::Vector{UInt8}`: if your sink defined `Data.weakrefstrings(sink::MySink) = true`, then it also needs to be able to accept the `reference` keyword argument, where a source's memory block will be passed, to be held onto appropriately by the sink when streaming WeakRefStrings. If a sink does not support streaming WeakRefStrings (the default), the sink constructor doesn't need to support any keyword arguments.
- `kwargs...`: Similar to `args...`, `kwargs...` is a catchall for any additional keyword arguments you'd like to expose for your sink constructor, typically matching supported keyword arguments that are provided through the normal sink constructor
- `sink`: in the 2nd form, an already-constructed sink is passed in as the 1st argument. This allows efficient sink re-use. This constructor needs to ensure the existing sink is modified (enlarged, shrunk, schema changes, etc) to be ready to accept the incoming source data as described by `sch`.

Now let's look at an example implementation from CSV.jl:

```julia
function CSV.Sink(fullpath::AbstractString; append::Bool=false, headers::Bool=true, colnames::Vector{String}=String[], kwargs...)
    io = IOBuffer()
    options = CSV.Options(kwargs...)
    !append && header && !isempty(colnames) && writeheaders(io, colnames, options)
    return CSV.Sink(options, io, fullpath, position(io), !append && header && !isempty(colnames), colnames, length(colnames), append)
end

function CSV.Sink(sch::Data.Schema, T, append, file::AbstractString; reference::Vector{UInt8}=UInt8[], kwargs...)
    sink = CSV.Sink(file; append=append, colnames=Data.header(sch), kwargs...)
    return sink
end

function CSV.Sink(sink, sch::Data.Schema, T, append; reference::Vector{UInt8}=UInt8[])
    sink.append = append
    sink.cols = size(sch, 2)
    !sink.header && !append && writeheaders(sink.io, Data.header(sch), sink.options, sink.quotefields)
    return sink
end
```

In this case, CSV.jl defined an initial constructor that just takes the filename with a few keyword arguments.
The two required Data.Sink constructors are then defined. The first constructs a new Sink, requiring a `file::AbstractString` argument.
We also see that `CSV.Sink` supports WeakRefString streaming by accepting a `reference` keyword argument (which is trivially implemented for CSV, since all data is simply written out to disk as text).

For the 2nd (last) constructor in the definitions above, we see the case where an existing `sink` is passed to `CSV.Sink`.
The sink updates a few of its fields (`sink.append = append`), and some logic is computed to determine if
the column headers should be written.
"""
abstract type Sink end

"""
`Data.streamtypes(::Type{[Sink]}) => Vector{StreamType}`

Returns a vector of `Data.StreamType`s that the sink is able to receive; the order of elements indicates the sink's streaming preference

For example, if my sink only supports `Data.Field` streaming, I would simply define:
```julia
Data.streamtypes(::Type{MyPkg.Sink}) = [Data.Field]
```

If, on the other hand, my sink also supported `Data.Column` streaming, and `Data.Column` streaming happend to be more efficient, I could define:
```julia
Data.streamtypes(::Type{MyPkg.Sink}) = [Data.Column, Data.Field] # put Data.Column first to indicate preference
```

A third option is a sink that operates on entire rows at a time, in which case I could define:
```julia
Data.streamtypes(::Type{MyPkg.Sink}) = [Data.Row]
```
The subsequent `Data.streamto!` method would then require the signature `Data.streamto!(sink::MyPkg.Sink, ::Type{Data.Row}, vals::NamedTuple, row, col, knownrows`
"""
function streamtypes end

"""
`Data.streamto!(sink, S::Type{StreamType}, val, row, col)`

`Data.streamto!(sink, S::Type{StreamType}, val, row, col, knownrows)`

Streams data to a sink. `S` is the type of streaming (`Data.Field`, `Data.Row`, or `Data.Column`). `val` is the value or values (single field, row as a NamedTuple, or column, respectively)
to be streamed to the sink. `row` and `col` indicate where the data should be streamed/stored.

A sink may optionally define the method that also accepts the `knownrows` argument, which will be `Val{true}` or `Val{false}`,
indicating whether the source streaming has a known # of rows or not. This can be useful for sinks that
may know how to pre-allocate space in the cases where the source can tell the # of rows, or in the case
of unknown # of rows, may need to stream the data in differently.
"""
function streamto! end
Data.streamto!(sink, S, val, row, col, knownrows) = Data.streamto!(sink, S, val, row, col)
Data.streamto!(sink, S, val, row, col) = Data.streamto!(sink, S, val, col)

# Optional methods
"""
`Data.cleanup!(sink)`

Sometimes errors occur during the streaming of data from source to sink. Some sinks may be left in
an undesired state if an error were to occur mid-streaming. `Data.cleanup!` allows a sink to "clean up"
any necessary resources in the case of a streaming error. `SQLite.jl`, for example, defines:
```julia
function Data.cleanup!(sink::SQLite.Sink)
    rollback(sink.db, sink.transaction)
    return
end
```
Since a database transaction is initiated at the start of streaming, it must be rolled back in the case of streaming error.

The default definition is: `Data.cleanup!(sink) = nothing`
"""
function cleanup! end

"""
`Data.close!(sink) => sink`

A function to "close" a sink to streaming. Some sinks require a definitive time where data can be "committed",
`Data.close!` allows a sink to perform any necessary resource management or commits to ensure all data that has been
streamed is stored appropriately. For example, the `SQLite` package defines:
```julia
function Data.close!(sink::SQLite.Sink)
    commit(sink.db, sink.transaction)
    return sink
end
```
Which commits a database transaction that was started when the sink was initially "opened".
"""
function close! end

# Generic fallbacks
cleanup!(sink) = nothing
close!(sink) = sink

"""
`Data.weakrefstrings(sink) => Bool`

If a sink is able to appropriately handle `WeakRefString` objects, it can define:
```julia
Data.weakrefstrings(::Type{[Sink]}) = true
```
to indicate that a source may stream those kinds of values to it. By default, sinks do
not support WeakRefString streaming. Supporting WeakRefStrings corresponds to accepting the
`reference` keyword argument in the required sink constructor method, see `?Data.Sink`.
"""
function weakrefstrings end
weakrefstrings(x) = false

# Data.stream!
"""
`Data.stream!(source, sink; append::Bool=false, transforms=Dict())`

`Data.stream!(source, ::Type{Sink}, args...; append::Bool=false, transforms=Dict(), kwargs...)`

Stream data from source to sink. The 1st definition assumes already constructed source & sink and takes two optional keyword arguments:

- `append::Bool=false`: whether the data from `source` should be appended to `sink`
- `transforms::Dict`: A dict with mappings between column # or name (Int or String) to a "transform" function. For `Data.Field` streaming, the transform function should be of the form `f(x::T) => y::S`, i.e. takes a single input of type `T` and returns a single value of type `S`. For `Data.Column` streaming, it should be of the form `f(x::AbstractVector{T}) => y::AbstractVector{S}`, i.e. take an AbstractVector with eltype `T` and return another AbstractVector with eltype `S`.

For the 2nd definition, the Sink type itself is passed as the 2nd argument (`::Type{Sink}`) and is constructed "on-the-fly", being passed `args...` and `kwargs...` like `Sink(args...; kwargs...)`.

While users are free to call `Data.stream!` themselves, oftentimes, packages want to provide even higher-level convenience functions.

An example of of these higher-level convenience functions are from CSV.jl:

```julia
function CSV.read(fullpath::Union{AbstractString,IO}, sink::Type=DataFrame, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}(), kwargs...)
    source = CSV.Source(fullpath; kwargs...)
    sink = Data.stream!(source, sink, args...; append=append, transforms=transforms, kwargs...)
    return Data.close!(sink)
end

function CSV.read(fullpath::Union{AbstractString,IO}, sink::T; append::Bool=false, transforms::Dict=Dict{Int,Function}(), kwargs...) where T
    source = CSV.Source(fullpath; kwargs...)
    sink = Data.stream!(source, sink; append=append, transforms=transforms)
    return Data.close!(sink)
end
```
In this example, CSV.jl defines it's own high-level function for reading from a `CSV.Source`. In these examples, a `CSV.Source` is constructed using the `fullpath` argument, along w/ any extra `kwargs...`.
The sink can be provided as a type with `args...` and `kwargs...` that will be passed to its DataStreams constructor, like `Sink(sch, streamtype, append, args...; kwargs...)`; otherwise, an already-constructed
Sink can be provided directly (2nd example).

Once the `source` is constructed, the data is streamed via the call to `Data.stream(source, sink; append=append, transforms=transforms)`, with the sink being returned.

And finally, to "finish" the streaming process, `Data.close!(sink)` is closed, which returns the finalized sink. Note that `Data.stream!(source, sink)` could be called multiple times with different sources and the same sink,
most likely with `append=true` being passed, to enable the accumulation of several sources into a single sink. A single `Data.close!(sink)` method should be called to officially close or commit the final sink.

Two "builtin" Source/Sink types that are included with the DataStreams package are the `Data.Table` and `Data.RowTable` types. `Data.Table` is a NamedTuple of AbstractVectors, with column names as NamedTuple fieldnames.
This type supports both `Data.Field` and `Data.Column` streaming. `Data.RowTable` is just a Vector of NamedTuples, and as such, only supports `Data.Field` streaming.

In addition, any `Data.Source` can be iterated via the `Data.rows(source)` function, which returns a NamedTuple-iterator over the rows of a source.
"""
function stream! end

# skipfield! and skiprow! only apply to Data.Field/Data.Row streaming
skipfield!(source, S, T, row, col) = Data.accesspattern(source) == Data.RandomAccess() ? nothing : Data.streamfrom(source, S, T, row, col)
function skiprow!(source, S, row, col)
    Data.accesspattern(source) == Data.RandomAccess() && return
    sch = Data.schema(source)
    cols = size(sch, 2)
    types = Data.types(sch)
    for i = col:cols
        Data.streamfrom(source, S, types[i], row, i)
    end
    return
end
function skiprows!(source, S, from, to)
    Data.accesspattern(source) == Data.RandomAccess() && return
    sch = Data.schema(source)
    cols = size(sch, 2)
    types = Data.types(sch)
    for row = from:to
        for col = 1:cols
            Data.streamfrom(source, S, types[col], row, col)
        end
    end
end

datatype(T) = Core.eval(parentmodule(Base.unwrap_unionall(T)), nameof(T))

include("namedtuples.jl")
include("query.jl")

end # module Data

using .Data
export Data

end # module DataStreams


text:
using ActorInterfaces.Classic
using ActorInterfacesTests
using QuickActors
import ActorInterfaces.Implementation.Tick

struct QuickActorsWrapper <: ActorInterfacesTests.ActorLib
    s::QuickScheduler
    QuickActorsWrapper() = new(QuickScheduler())
end

function ActorInterfacesTests.ex_spawn!(lib::QuickActorsWrapper, bhv, args...)
    return spawn!(lib.s, bhv, args...)
end

function ActorInterfacesTests.ex_send!(lib::QuickActorsWrapper, addr::Addr, msg...)
    return send!(lib.s, addr, msg...)
end

function ActorInterfacesTests.ex_actorcount(lib::QuickActorsWrapper)
    return length(lib.s.actorcache)
end

function ActorInterfacesTests.ex_runtofinish(lib::QuickActorsWrapper)
    while Tick.tick!(lib.s) end
    return nothing
end

run_suite(QuickActorsWrapper)

ext: jl


text:
using DReal

## Function Approximation
## ======================

f(x) = x*2 / sin(x)
θ = [Var(Float64, -10.0, 10.0) for i = 1:7]
# g(x) = θ[1] + θ[2]*x + θ[3]*x^2 + θ[4]*x^3

g(x) = θ[1] + θ[2]*x^θ[3] + θ[4]*x^θ[5] + θ[6]*x^θ[7]

x = ForallVar(Float64, 0.0, 10.0)
constraint = g(x) == f(x)
constraint = (g(x) * f(x)) < 10
constraint = abs(g(x) - f(x)) < 1.0 # However, this (incorrectly) is UNSAT
add!(constraint)
@show is_satisfiable()
# result = model(θ...)
# @show result


text:
module lol
export tria6,N,Nl,J,Nxy,ke,me,ε,Ematrix,BC
function tria6(l,xs,ys)
    (l1,l2,l3) = l
    N1 = (2*l1-1)*l1; N2 = (2*l2-1)*l2; N3 = (2*l3-1)*l3;
    N4 = 4*l1*l2; N5 = 4*l2*l3; N6 = 4*l3*l1;
    x = dot(xs,[N1 N2 N3 N4 N5 N6])
    y = dot(ys,[N1 N2 N3 N4 N5 N6])
    return x,y
end
function N(l)
    (l1,l2,l3) = l;
    N1 = (2*l1-1)*l1; N2 = (2*l2-1)*l2; N3 = (2*l3-1)*l3;
    N4 = 4*l1*l2; N5 = 4*l2*l3; N6 = 4*l3*l1;

    N = [N1 0 N2 0 N3 0 N4 0 N5 0 N6 0;
         0 N1 0 N2 0 N3 0 N4 0 N5 0 N6];
    return N
end
function J(l,xs,ys)
    (l1,l2,l3) = l
    Nl1 = [4*l1-1 0 0 4*l2 0 4*l3]
    Nl2 = [0 4*l2-1 0 4*l1 4*l3 0]
    Nl3 = [0 0 4*l3-1 0 4*l2 4*l1]
    Nl = [Nl1; Nl2; Nl3]
    J = Nl*[xs' ys']
end
function Nl(l)
    (l1,l2,l3) = l
    Nl1 = [4*l1-1 0 0 4*l2 0 4*l3]
    Nl2 = [0 4*l2-1 0 4*l1 4*l3 0]
    Nl3 = [0 0 4*l3-1 0 4*l2 4*l1]
    Nl = [Nl1; Nl2; Nl3]
end
function Nxy(l,xs,ys)
    (l1,l2,l3) = l
    j = J(l,xs,ys)
    nl = Nl(l)
    nxy = j\nl
    return nxy
end

function me(xs,ys)
    #neeeds multiplication by rho*h
    me  = zeros(12,12)
    a1 = 0.0597158717; b1 = 0.4701420641; a2 = 0.7974269853; b2 = 0.1012865073
    w1 = 0.2250000000; w2 =0.1323941527; w3=0.1323941527; w4=0.1323941527
    w5 =0.1259391805; w6 = 0.1259391805; w7 = 0.1259391805;
    g1 = [1 1 1]./3; g2 = [a1 b1 b1]; g3 = [b1 a1 b1]; g4 = [b1 b1 a1];
    g5 = [a2 b2 b2]; g6 = [b2 a2 b2]; g7 = [b2 b2 a2]
    w = [w1 w2 w3 w4 w5 w6 w7]
    g = [[g1] [g2] [g3] [g4] [g5] [g6] [g7]]

    for i=1:7
        n = N(g[i])
        me = me + w[i]*n'*n
    end
    return me
end
function Ematrix(e,v,cond)
    if cond=="stress"
        E = (e/(1-v^2))*[1 v 0;
                        v 1 0;
                        0 0 (1-v)/2]
    elseif cond=="strain"
        E = (e*(1-v)/((1+v)*(1-2*v)))*[1 v/(v-1) 0;
                                       v/(v-1) 1 0;
                                        0 0 (1-2*v)/(2-2*v)]
    else
        println("motherfucker no planestress or planestrain condition")
        exit()
    end
    return E
end
function ke(E,xs,ys)
    #neeeds multiplication by h
    me  = zeros(12,12)
    a1 = 0.0597158717; b1 = 0.4701420641; a2 = 0.7974269853; b2 = 0.1012865073
    w1 = 0.2250000000; w2 =0.1323941527; w3=0.1323941527; w4=0.1323941527
    w5 =0.1259391805; w6 = 0.1259391805; w7 = 0.1259391805;
    g1 = [1 1 1]./3; g2 = [a1 b1 b1]; g3 = [b1 a1 b1]; g4 = [b1 b1 a1];
    g5 = [a2 b2 b2]; g6 = [b2 a2 b2]; g7 = [b2 b2 a2]
    w = [w1 w2 w3 w4 w5 w6 w7]
    g = [[g1] [g2] [g3] [g4] [g5] [g6] [g7]]

    ke = zeros(12,12)
    for i=1:7
        nxy = Nxy(g[i],xs,ys)
        nx = nxy[1,:]; ny = nxy[2,:];
        nx = [1 0 0 0 0 0;0 0 0 0 0 0;0 1 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 1 0;0 0 0 0 0 0;0 0 0 0 0 1;0 0 0 0 0 0]*nx
        ny = [0 0 0 0 0 0;1 0 0 0 0 0;0 0 0 0 0 0;0 1 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 1 0;0 0 0 0 0 0;0 0 0 0 0 1]*ny;
        b = [1 0;0 0;0 1]*[nx';ny'] + [0 0;0 1;1 0]*[nx';ny']

        ke = ke + w[i]*b'*E*b
    end
    return ke
end
function BC(f,xedge,yedge,map)
    #f is vertical array
    #needs to multiplied by h
    #(fx,fy)
    if map[1]==0
        l1 = zeros(7,1)
        l2 = [-0.949108  -0.741531  -0.405845  0.405845  0.741531  0.949108  0.0]
        w = [0.129485  0.279705  0.38183  0.38183  0.279705  0.129485  0.417959]
        #transform to 0~1
        l2 = 0.5 .*l2 .+ 0.5
        l3 = 1 .- l2
        w = 0.5 .* w
    elseif map[2]==0
        l2 = zeros(7,1)
        l3 = [-0.949108  -0.741531  -0.405845  0.405845  0.741531  0.949108  0.0]
        w = [0.129485  0.279705  0.38183  0.38183  0.279705  0.129485  0.417959]
        #transform to 0~1
        l3 = 0.5 .*l3 .+ 0.5
        l1 = 1 .- l3
        w = 0.5 .* w
    elseif map[3]==0
        l3 = zeros(7,1)
        l1 = [-0.949108  -0.741531  -0.405845  0.405845  0.741531  0.949108  0.0]
        w = [0.129485  0.279705  0.38183  0.38183  0.279705  0.129485  0.417959]
        #transform to 0~1
        l1 = 0.5 .*l1 .+ 0.5
        l2 = 1 .- l1
        w = 0.5 .* w
    else
        println("this is not a valid edge map")
        exit()
    end
    bc = zeros(12,1)
    for i=1:7
        l = (l1[i],l2[i],l3[i])
        (x,y) = tria6(l,xedge,yedge)
        bc = bc + N(l)'*f([x,y])
    end
    return bc
end
function ε(l,xs,ys)
    nxy = Nxy(l,xs,ys)
    nx = nxy[1,:]; ny = nxy[2,:];
    nx = [1 0 0 0 0 0;0 0 0 0 0 0;0 1 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 1 0;0 0 0 0 0 0;0 0 0 0 0 1;0 0 0 0 0 0]*nx
    ny = [0 0 0 0 0 0;1 0 0 0 0 0;0 0 0 0 0 0;0 1 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 1 0;0 0 0 0 0 0;0 0 0 0 0 1]*ny;
    nxy = [nx';ny']
    xy = [collect(zip(xs,ys))[i][k] for i=1:6 for k=1:2]
    b = [1 0;0 0;0 1]*[nx';ny'] + [0 0;0 1;1 0]*[nx';ny']
    epsilon = b*xy
    return epsilon
end
end

end

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
# f(x)=diagA*X
function FdiagAmultX(A::Array,X::Array)
    if size(A)==(1,1)
        return (A[1].*X,nothing)
    elseif size(X)==(1,1)
        return (diagm(vec(A)).*X[1],nothing)
    else
        return (diagm(vec(A))*X,nothing)
    end
end

function FdiagAmultX_inplace(value,auxvalue,A::Array,X::Array)
    if size(A)==(1,1)
        copy!(value,A[1]*X)
    elseif size(X)==(1,1)
        copy!(value,diagm(vec(A))*X[1])
    else
        copy!(value,diagm(vec(A))*X) # BLAS has no dgmm routine
    end
end

function DdiagAmultX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X)
    if derivativeIDX==1
        if size(A)==(1,1)
            axpy!(1.0,[sum(X.*grad_c)],grad_n)
        elseif size(X)==(1,1)
            axpy!(X[1],sum(grad_c,2),grad_n)
        else
            axpy!(1.0,sum(X.*grad_c,2),grad_n)
        end
    elseif derivativeIDX==2
        if size(A)==(1,1)
            axpy!(A[1],grad_c,grad_n)
        elseif size(X)==(1,1)
            #axpy!(1.0,[sum(A.*diag(grad_c))],grad_n)
            axpy!(1.0,[sum(sum(A,2).*diag(grad_c))],grad_n)
        else
            axpy!(1.0,diagm(vec(A))*grad_c,grad_n)
        end
    end
end

if PROC=="GPU"
    function FdiagAmultX_inplace(value::CudaArray,auxvalue,A::CudaArray,X::CudaArray)
        if size(A)==(1,1)
            gax!(A,X,value)
        elseif size(X)==(1,1)
            diagm!(A,value)
            scale!(X,value) # nb argument converse of Base.scale!
        else
            CUBLAS.dgmm!('L',X,vec(A),value)
        end
    end

    function DdiagAmultX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,X::CudaArray)
        if derivativeIDX==1
            if size(A)==(1,1)
                tmp=CudaArray(Float64,size(X))
                vmult!(1.0,X,grad_c,tmp)
                sum_update!(1.0,tmp,1.0,grad_n)
                free(tmp)
            elseif size(X)==(1,1)
                tmp=CudaArray(Float64,size(grad_n))
                sum!(grad_c,2,0.0,tmp) # tmp=sum(grad_c,2)
                alphaaxpy!(1.0,X,tmp,grad_n)
                free(tmp)
            else
                tmp=CudaArray(Float64,size(X))
                vmult!(1.0,X,grad_c,tmp) #tmp=X.*grad_c
                sum!(tmp,2,1.0,grad_n)
            end
        elseif derivativeIDX==2
            if size(A)==(1,1)
                alphaaxpy!(1.0,A,grad_c,grad_n)
            elseif size(X)==(1,1)
              tmp1=CudaArray(Float64,(size(A,1),1))
              sum!(A,2,0.0,tmp1) # sum(A,2)
              tmp2=CudaArray(Float64,(size(A,1),1))
              diag!(grad_c,tmp2)
              vmult!(1.0,tmp1,tmp2,tmp1)
              sum_update!(1.0,tmp1,1.0,grad_n)
              free(tmp1); free(tmp2)
#            axpy!(1.0,[sum(sum(A,2).*diag(grad_c))],grad_n)
            else
                tmp=CudaArray(Float64,size(grad_c))
                             CUBLAS.dgmm!('L',grad_c,vec(A),tmp)
                axpy!(1.0,tmp,grad_n)
                free(tmp)
            end
        end
    end
end

Derivative[FdiagAmultX]=DdiagAmultX
export FdiagAmultX

Inplace[FdiagAmultX]=FdiagAmultX_inplace

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# Rotation utilities

export RotationLinOp, RotationParametericLinOp, RotationParametericDelayedEval, JacobianRotationEvaluated
export rotation_linop, rotation, Jacobian, ∂

## Linear operator

struct RotationLinOp{T}<:AbstractLinearOperator{T,3,3}
    R::AbstractArray{T,3}
    c::Union{Nothing,AbstractArray{T,2}}
    s::Union{Nothing,AbstractArray{T,2}}
end

rotation_linop(R::AbstractArray{T,3}; c::Union{Nothing,AbstractArray{T,2}}=nothing, s::Union{Nothing,AbstractArray{T,2}}=nothing) where {T<:Real} = RotationLinOp{T}(R, c, s)

function rotation_linop(φ::AbstractArray{T,2}; derivative::Bool=false) where {T<:Real}
    c = cos.(φ)
    s = sin.(φ)
    Rout = rotation_matrix(c, s; derivative=derivative)
    ~derivative ? (return RotationLinOp{T}(Rout, c, s)) : (return (RotationLinOp{T}(Rout[1], c, s), Rout[2]))
end

function rotation_matrix(c::AbstractArray{T,2}, s::AbstractArray{T,2}; derivative::Bool=false) where {T<:Real}
    c1 = c[:,1]; s1 = s[:,1]
    c2 = c[:,2]; s2 = s[:,2]
    c3 = c[:,3]; s3 = s[:,3]

    R = similar(c1, length(c1), 3, 3)
    R[:,1,1] .= c1.*c2;            R[:,1,2] .= -c2.*s1;           R[:,1,3] .= -s2;
    R[:,2,1] .= c3.*s1-c1.*s2.*s3; R[:,2,2] .= c1.*c3+s1.*s2.*s3; R[:,2,3] .= -c2.*s3;
    R[:,3,1] .= c1.*c3.*s2+s1.*s3; R[:,3,2] .= c1.*s3-c3.*s1.*s2; R[:,3,3] .= c2.*c3;
    R = permutedims(R, (1, 3, 2))

    if derivative
        ∂R = similar(c1, length(c1), 3, 3, 3)

        ∂R[:,1,1,1] .= -s1.*c2;            ∂R[:,1,2,1] .= -c2.*c1;            ∂R[:,1,3,1] .= 0;
        ∂R[:,2,1,1] .= c3.*c1+s1.*s2.*s3;  ∂R[:,2,2,1] .= -s1.*c3+c1.*s2.*s3; ∂R[:,2,3,1] .= 0;
        ∂R[:,3,1,1] .= -s1.*c3.*s2+c1.*s3; ∂R[:,3,2,1] .= -s1.*s3-c3.*c1.*s2; ∂R[:,3,3,1] .= 0;

        ∂R[:,1,1,2] .= -c1.*s2;     ∂R[:,1,2,2] .= s2.*s1;      ∂R[:,1,3,2] .= -c2;
        ∂R[:,2,1,2] .= -c1.*c2.*s3; ∂R[:,2,2,2] .= s1.*c2.*s3;  ∂R[:,2,3,2] .= s2.*s3;
        ∂R[:,3,1,2] .= c1.*c3.*c2;  ∂R[:,3,2,2] .= -c3.*s1.*c2; ∂R[:,3,3,2] .= -s2.*c3;

        ∂R[:,1,1,3] .= 0;                  ∂R[:,1,2,3] .= 0;                  ∂R[:,1,3,3] .= 0;
        ∂R[:,2,1,3] .= -s3.*s1-c1.*s2.*c3; ∂R[:,2,2,3] .= -c1.*s3+s1.*s2.*c3; ∂R[:,2,3,3] .= -c2.*c3;
        ∂R[:,3,1,3] .= -c1.*s3.*s2+s1.*c3; ∂R[:,3,2,3] .= c1.*c3+s3.*s1.*s2;  ∂R[:,3,3,3] .= -c2.*s3;

        ∂R = permutedims(∂R, (1, 3, 2, 4))
    end

    ~derivative ? (return R) : (return (R, ∂R))
end

AbstractLinearOperators.domain_size(R::RotationLinOp) = (size(R.R,1), "nk", 3)
AbstractLinearOperators.range_size(R::RotationLinOp) = domain_size(R)

function AbstractLinearOperators.matvecprod(R::RotationLinOp{T}, K::AbstractArray{T,3}) where {T<:Real}
    nt, nk, _ = size(K)
    return sum(reshape(R.R, nt, 1, 3, 3).*reshape(K, nt, nk, 1, 3); dims=4)[:,:,:,1]
end

function AbstractLinearOperators.matvecprod_adj(R::RotationLinOp{T}, K::AbstractArray{T,3}) where {T<:Real}
    nt, nk, _ = size(K)
    return sum(reshape(R.R, nt, 1, 3, 3).*reshape(K, nt, nk, 3, 1); dims=3)[:,:,1,:]
end

## Parameteric linear operator

struct RotationParametericLinOp end

rotation() = RotationParametericLinOp()

(::RotationParametericLinOp)(φ::AbstractArray{T,2}) where {T<:Real} = rotation_linop(φ)

(R::RotationParametericLinOp)() = R

## Parameteric delayed evaluation

struct RotationParametericDelayedEval{T}
    K::AbstractKSpaceSampling{T}
end

Base.:*(::RotationParametericLinOp, K::AbstractKSpaceSampling{T}) where {T<:Real} =  RotationParametericDelayedEval{T}(K)

(RK::RotationParametericDelayedEval{T})(φ::AbstractArray{T,2}) where {T<:Real} = rotation_linop(φ)*coord(RK.K)

## Jacobian of rotation evaluated

struct JacobianRotationEvaluated{T}<:AbstractLinearOperator{Complex{T},2,3}
    K::AbstractKSpaceSampling{T}
    ∂R::AbstractArray{T,4}
end

function Jacobian(RK::RotationParametericDelayedEval{T}, φ::AbstractArray{T,2}) where {T<:Real}
    Rφ, ∂Rφ = rotation_linop(φ; derivative=true)
    return Rφ*coord(RK.K), JacobianRotationEvaluated{T}(RK.K, ∂Rφ)
end

∂(RK::RotationParametericDelayedEval{T}, φ::AbstractArray{T,2}) where {T<:Real} = Jacobian(RK, φ)

AbstractLinearOperators.domain_size(∂RK::JacobianRotationEvaluated) = (size(∂RK.∂R,1),3)
AbstractLinearOperators.range_size(∂RK::JacobianRotationEvaluated) = (size(∂RK.K)...,3)

function AbstractLinearOperators.matvecprod(∂RK::JacobianRotationEvaluated{T}, Δφ::AbstractArray{CT,2}) where {T<:Real,CT<:RealOrComplex{T}}
    ΔR = sum(∂RK.∂R.*reshape(Δφ, :, 1, 1, 3); dims=4)[:,:,:,1]
    return rotation_linop(ΔR)*coord(∂RK.K)
end
Base.:*(∂RK::JacobianRotationEvaluated{T}, Δφ::AbstractArray{T,2}) where {T<:Real} = AbstractLinearOperators.matvecprod(∂RK, Δφ)

function AbstractLinearOperators.matvecprod_adj(∂RK::JacobianRotationEvaluated{T}, ΔRK::AbstractArray{CT,3}) where {T<:Real,CT<:RealOrComplex{T}}
    nt, nk, _ = size(ΔRK)
    KΔRK = sum(ΔRK.*reshape(coord(∂RK.K), nt, nk, 1, 3); dims=2)[:,1,:,:]
    return sum(∂RK.∂R.*reshape(KΔRK, nt, 3, 3, 1); dims=2:3)[:,1,1,:]
end

function LinearAlgebra.dot(∇u::AbstractArray{Complex{T},3}, ∂RK::JacobianRotationEvaluated{T}) where {T<:Real}
    nt, nk, _ = size(∇u)
    ∂R = ∂RK.∂R
    K = coord(∂RK.K)
    J = similar(∇u, nt, nk, 3)
    ∂jRK = similar(∇u, nt, nk, 3)
    @inbounds for j = 1:3
        ∂jRK .= sum(reshape(∂R[:,:,:,j], nt, 1, 3, 3).*reshape(K, nt, nk, 1, 3); dims=4)[:,:,:,1]
        J[:,:,j] .= sum(∇u.*∂jRK; dims=3)[:,:,1]
    end
    return J
end

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
# Ordinary Maximum Likelihood Estimation

############################################################################################
### Types
############################################################################################
"""
Maximum likelihood estimation.

# Constructor

    SemML(;observed, meanstructure = false, approximate_hessian = false, kwargs...)

# Arguments
- `observed::SemObserved`: the observed part of the model
- `meanstructure::Bool`: does the model have a meanstructure?
- `approximate_hessian::Bool`: if hessian-based optimization is used, should the hessian be swapped for an approximation

# Examples
```julia
my_ml = SemML(observed = my_observed)
```

# Interfaces
Analytic gradients are available, and for models without a meanstructure, also analytic hessians.

# Extended help
## Implementation
Subtype of `SemLossFunction`.
"""
struct SemML{INV,M,M2,B, V} <: SemLossFunction
    Σ⁻¹::INV
    Σ⁻¹Σₒ::M
    meandiff::M2
    approximate_hessian::B
    has_meanstructure::V
end

############################################################################################
### Constructors
############################################################################################

function SemML(;observed, meanstructure = false, approximate_hessian = false, kwargs...)
    isnothing(obs_mean(observed)) ?
        meandiff = nothing :
        meandiff = copy(obs_mean(observed))
    return SemML(
        similar(obs_cov(observed)),
        similar(obs_cov(observed)),
        meandiff,
        approximate_hessian,
        Val(meanstructure)
        )
end

############################################################################################
### objective, gradient, hessian methods
############################################################################################

# first, dispatch for meanstructure
objective!(semml::SemML, par, model::AbstractSemSingle) =
    objective!(semml::SemML, par, model, semml.has_meanstructure, imply(model))
gradient!(semml::SemML, par, model::AbstractSemSingle) =
    gradient!(semml::SemML, par, model, semml.has_meanstructure, imply(model))
hessian!(semml::SemML, par, model::AbstractSemSingle) =
    hessian!(semml::SemML, par, model, semml.has_meanstructure, imply(model))
objective_gradient!(semml::SemML, par, model::AbstractSemSingle) =
    objective_gradient!(semml::SemML, par, model, semml.has_meanstructure, imply(model))
objective_hessian!(semml::SemML, par, model::AbstractSemSingle) =
    objective_hessian!(semml::SemML, par, model, semml.has_meanstructure, imply(model))
gradient_hessian!(semml::SemML, par, model::AbstractSemSingle) =
    gradient_hessian!(semml::SemML, par, model, semml.has_meanstructure, imply(model))
objective_gradient_hessian!(semml::SemML, par, model::AbstractSemSingle) =
    objective_gradient_hessian!(semml::SemML, par, model, semml.has_meanstructure, imply(model))

############################################################################################
### Symbolic Imply Types

function objective!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{T},
        imp::SemImplySymbolic) where T

    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml),
        Σ⁻¹ = Σ⁻¹(semml), μ = μ(imply(model)), μₒ = obs_mean(observed(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol) return non_posdef_return(par) end

        ld = logdet(Σ_chol)
        Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
        #mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)

        if T
            μ₋ = μₒ - μ
            return ld + dot(Σ⁻¹, Σₒ) + dot(μ₋, Σ⁻¹, μ₋)
        else
            return ld + dot(Σ⁻¹, Σₒ)
        end
    end
end

function gradient!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{T},
        imp::SemImplySymbolic) where T

    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml),
        Σ⁻¹ = Σ⁻¹(semml), ∇Σ = ∇Σ(imply(model)),
        μ = μ(imply(model)), ∇μ = ∇μ(imply(model)), μₒ = obs_mean(observed(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol)
            return ones(eltype(par), size(par))
        end

        Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
        mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)

        if T
            μ₋ = μₒ - μ
            μ₋ᵀΣ⁻¹ = μ₋'*Σ⁻¹
            gradient = vec(Σ⁻¹*(I - Σₒ*Σ⁻¹ - μ₋*μ₋ᵀΣ⁻¹))'*∇Σ - 2*μ₋ᵀΣ⁻¹*∇μ
            return gradient'
        else
            gradient = (vec(Σ⁻¹)-vec(Σ⁻¹Σₒ*Σ⁻¹))'*∇Σ
            return gradient'
        end
    end
end

function hessian!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{false},
        imp::SemImplySymbolic)

    let Σ = Σ(imply(model)), ∇Σ = ∇Σ(imply(model)), Σₒ = obs_cov(observed(model)),
        Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml), Σ⁻¹ = Σ⁻¹(semml),
        ∇²Σ_function! = ∇²Σ_function(imply(model)), ∇²Σ = ∇²Σ(imply(model))

    copyto!(Σ⁻¹, Σ)
    Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

    if !isposdef(Σ_chol)
        return diagm(fill(one(eltype(par)), length(par)))
    end

    Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)

    if semml.approximate_hessian
        hessian = 2*∇Σ'*kron(Σ⁻¹, Σ⁻¹)*∇Σ
    else
        mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)
        Σ⁻¹ΣₒΣ⁻¹ = Σ⁻¹Σₒ*Σ⁻¹
        # inner
        J = vec(Σ⁻¹ - Σ⁻¹ΣₒΣ⁻¹)'
        ∇²Σ_function!(∇²Σ, J, par)
        # outer
        H_outer = 2*kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹)
        hessian = ∇Σ'*H_outer*∇Σ
        hessian += ∇²Σ
    end

    return hessian
    end
end

function hessian!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{true},
        imp::SemImplySymbolic)
    throw(DomainError(H, "hessian of ML + meanstructure is not available"))
end

function objective_gradient!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{T},
        imp::SemImplySymbolic) where T

    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml),
        Σ⁻¹ = Σ⁻¹(semml), μ = μ(imply(model)), μₒ = obs_mean(observed(model)),
        ∇Σ = ∇Σ(imply(model)), ∇μ = ∇μ(imply(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol)
            return non_posdef_return(par), ones(eltype(par), size(par))
        else
            ld = logdet(Σ_chol)
            Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
            mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)

            if T
                μ₋ = μₒ - μ
                μ₋ᵀΣ⁻¹ = μ₋'*Σ⁻¹

                objective = ld + tr(Σ⁻¹Σₒ) + dot(μ₋, Σ⁻¹, μ₋)
                gradient = vec(Σ⁻¹*(I - Σₒ*Σ⁻¹ - μ₋*μ₋ᵀΣ⁻¹))'*∇Σ - 2*μ₋ᵀΣ⁻¹*∇μ
                return objective, gradient'
            else
                objective = ld + tr(Σ⁻¹Σₒ)
                gradient = (vec(Σ⁻¹)-vec(Σ⁻¹Σₒ*Σ⁻¹))'*∇Σ
                return objective, gradient'
            end
        end
    end
end

function objective_hessian!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{T},
        imp::SemImplySymbolic) where T

    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml),
        Σ⁻¹ = Σ⁻¹(semml), ∇Σ = ∇Σ(imply(model)), ∇μ = ∇μ(imply(model)),
        ∇²Σ_function! = ∇²Σ_function(imply(model)), ∇²Σ = ∇²Σ(imply(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol)
            return non_posdef_return(par), diagm(fill(one(eltype(par)), length(par)))
        else
            ld = logdet(Σ_chol)
            Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
            mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)
            objective = ld + tr(Σ⁻¹Σₒ)

            if semml.approximate_hessian
                hessian = 2*∇Σ'*kron(Σ⁻¹, Σ⁻¹)*∇Σ
            else
                Σ⁻¹ΣₒΣ⁻¹ = Σ⁻¹Σₒ*Σ⁻¹
                # inner
                J = vec(Σ⁻¹ - Σ⁻¹ΣₒΣ⁻¹)'
                ∇²Σ_function!(∇²Σ, J, par)
                # outer
                H_outer = 2*kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹)
                hessian = ∇Σ'*H_outer*∇Σ
                hessian += ∇²Σ
            end

            return objective, hessian
        end
    end
end

function objective_hessian!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{true},
        imp::SemImplySymbolic)
    throw(DomainError(H, "hessian of ML + meanstructure is not available"))
end

function gradient_hessian!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{false},
        imp::SemImplySymbolic)

    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml),
        Σ⁻¹ = Σ⁻¹(semml), ∇Σ = ∇Σ(imply(model)), ∇μ = ∇μ(imply(model)),
        ∇²Σ_function! = ∇²Σ_function(imply(model)), ∇²Σ = ∇²Σ(imply(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol)
            return ones(eltype(par), size(par)), diagm(fill(one(eltype(par)), length(par)))
        end

        Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
        mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)

        Σ⁻¹ΣₒΣ⁻¹ = Σ⁻¹Σₒ*Σ⁻¹

        J = vec(Σ⁻¹ - Σ⁻¹ΣₒΣ⁻¹)'
        gradient = J*∇Σ

        if semml.approximate_hessian
            hessian = 2*∇Σ'*kron(Σ⁻¹, Σ⁻¹)*∇Σ
        else
            # inner
            ∇²Σ_function!(∇²Σ, J, par)
            # outer
            H_outer = 2*kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹)
            hessian = ∇Σ'*H_outer*∇Σ
            hessian += ∇²Σ
        end

        return gradient', hessian
    end
end

function gradient_hessian!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{true},
        imp::SemImplySymbolic)
    throw(DomainError(H, "hessian of ML + meanstructure is not available"))
end

function objective_gradient_hessian!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{false},
        imp::SemImplySymbolic)

    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml),
        Σ⁻¹ = Σ⁻¹(semml), ∇Σ = ∇Σ(imply(model)),
        ∇²Σ_function! = ∇²Σ_function(imply(model)), ∇²Σ = ∇²Σ(imply(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol)
            objective = non_posdef_return(par)
            gradient = ones(eltype(par), size(par))
            hessian = diagm(fill(one(eltype(par)), length(par)))
            return objective, gradient, hessian
        end

        ld = logdet(Σ_chol)
        Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
        mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)
        objective = ld + tr(Σ⁻¹Σₒ)

        Σ⁻¹ΣₒΣ⁻¹ = Σ⁻¹Σₒ*Σ⁻¹

        J = vec(Σ⁻¹ - Σ⁻¹ΣₒΣ⁻¹)'
        gradient = J*∇Σ

        if semml.approximate_hessian
            hessian = 2*∇Σ'*kron(Σ⁻¹, Σ⁻¹)*∇Σ
        else
            Σ⁻¹ΣₒΣ⁻¹ = Σ⁻¹Σₒ*Σ⁻¹
            # inner
            ∇²Σ_function!(∇²Σ, J, par)
            # outer
            H_outer = 2*kron(Σ⁻¹ΣₒΣ⁻¹, Σ⁻¹) - kron(Σ⁻¹, Σ⁻¹)
            hessian = ∇Σ'*H_outer*∇Σ
            hessian += ∇²Σ
        end

        return objective, gradient', hessian
    end
end

function objective_gradient_hessian!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{true},
        imp::SemImplySymbolic)
    throw(DomainError(H, "hessian of ML + meanstructure is not available"))
end

############################################################################################
### Non-Symbolic Imply Types

# no hessians ------------------------------------------------------------------------------

function hessian!(semml::SemML, par, model::AbstractSemSingle, has_meanstructure, imp::RAM)
    throw(DomainError(H, "hessian of ML + non-symbolic imply type is not available"))
end

function objective_hessian!(semml::SemML, par, model::AbstractSemSingle, has_meanstructure, imp::RAM)
    throw(DomainError(H, "hessian of ML + non-symbolic imply type is not available"))
end

function gradient_hessian!(semml::SemML, par, model::AbstractSemSingle, has_meanstructure, imp::RAM)
    throw(DomainError(H, "hessian of ML + non-symbolic imply type is not available"))
end

function objective_gradient_hessian!(semml::SemML, par, model::AbstractSemSingle, has_meanstructure, imp::RAM)
    throw(DomainError(H, "hessian of ML + non-symbolic imply type is not available"))
end

# objective, gradient ----------------------------------------------------------------------

function objective!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{T},
        imp::RAM) where T
    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml),
        Σ⁻¹ = Σ⁻¹(semml), μ = μ(imply(model)), μₒ = obs_mean(observed(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol) return non_posdef_return(par) end

        ld = logdet(Σ_chol)
        Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
        mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)

        if T
            μ₋ = μₒ - μ
            return ld + tr(Σ⁻¹Σₒ) + dot(μ₋, Σ⁻¹, μ₋)
        else
            return ld + tr(Σ⁻¹Σₒ)
        end
    end
end

function gradient!(semml::SemML, par, model::AbstractSemSingle, has_meanstructure::Val{T}, imp::RAM) where T

    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml),
        Σ⁻¹ = Σ⁻¹(semml), S = S(imply(model)), M = M(imply(model)),
        F⨉I_A⁻¹ = F⨉I_A⁻¹(imply(model)), I_A⁻¹ = I_A⁻¹(imply(model)),
        ∇A = ∇A(imply(model)), ∇S = ∇S(imply(model)), ∇M = ∇M(imply(model)),
        μ = μ(imply(model)), μₒ = obs_mean(observed(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol) return ones(eltype(par), size(par)) end

        Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
        #mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)

        C = F⨉I_A⁻¹'*(I-Σₒ*Σ⁻¹)'*Σ⁻¹*F⨉I_A⁻¹
        gradient = 2vec(C*S*I_A⁻¹')'∇A + vec(C)'∇S

        if T
            μ₋ = μₒ - μ
            μ₋ᵀΣ⁻¹ = μ₋'*Σ⁻¹
            k = μ₋ᵀΣ⁻¹*F⨉I_A⁻¹

            gradient += -2k*∇M - 2vec(k'M'I_A⁻¹')'∇A - 2vec(k'k*S*I_A⁻¹')'∇A - vec(k'k)'∇S
        end

        return gradient'
    end
end

function objective_gradient!(
        semml::SemML,
        par,
        model::AbstractSemSingle,
        has_meanstructure::Val{T},
        imp::RAM) where T

    let Σ = Σ(imply(model)), Σₒ = obs_cov(observed(model)), Σ⁻¹Σₒ =  Σ⁻¹Σₒ(semml), Σ⁻¹ = Σ⁻¹(semml),
        S = S(imply(model)), M = M(imply(model)), F⨉I_A⁻¹ = F⨉I_A⁻¹(imply(model)), I_A⁻¹ = I_A⁻¹(imply(model)),
        ∇A = ∇A(imply(model)), ∇S = ∇S(imply(model)), ∇M = ∇M(imply(model)),
        μ = μ(imply(model)), μₒ = obs_mean(observed(model))

        copyto!(Σ⁻¹, Σ)
        Σ_chol = cholesky!(Symmetric(Σ⁻¹); check = false)

        if !isposdef(Σ_chol)
            objective = non_posdef_return(par)
            gradient = ones(eltype(par), size(par))
            return objective, gradient
        else
            ld = logdet(Σ_chol)
            Σ⁻¹ .= LinearAlgebra.inv!(Σ_chol)
            mul!(Σ⁻¹Σₒ, Σ⁻¹, Σₒ)
            objective = ld + tr(Σ⁻¹Σₒ)

            C = F⨉I_A⁻¹'*(I-Σₒ*Σ⁻¹)'*Σ⁻¹*F⨉I_A⁻¹
            gradient = 2vec(C*S*I_A⁻¹')'∇A + vec(C)'∇S

            if T
                μ₋ = μₒ - μ
                objective += dot(μ₋, Σ⁻¹, μ₋)

                μ₋ᵀΣ⁻¹ = μ₋'*Σ⁻¹
                k = μ₋ᵀΣ⁻¹*F⨉I_A⁻¹
                gradient += -2k*∇M - 2vec(k'M'I_A⁻¹')'∇A - 2vec(k'k*S*I_A⁻¹')'∇A - vec(k'k)'∇S
            end

            return objective, gradient'
        end
    end
end

############################################################################################
### additional functions
############################################################################################

function non_posdef_return(par)
    if eltype(par) <: AbstractFloat
        return floatmax(eltype(par))
    else
        return typemax(eltype(par))
    end
end

############################################################################################
### recommended methods
############################################################################################

update_observed(lossfun::SemML, observed::SemObservedMissing; kwargs...) =
    throw(ArgumentError("ML estimation does not work with missing data - use FIML instead"))

function update_observed(lossfun::SemML, observed::SemObserved; kwargs...)
    if size(lossfun.Σ⁻¹) == size(obs_cov(observed))
        return lossfun
    else
        return SemML(;observed = observed, kwargs...)
    end
end

############################################################################################
### additional methods
############################################################################################

Σ⁻¹(semml::SemML) = semml.Σ⁻¹
Σ⁻¹Σₒ(semml::SemML) = semml.Σ⁻¹Σₒ


text:
using AQFED, Test, StatsBase
import AQFED.Black

@testset "JaeckelAccuracyFixedStrike" begin
  strike = 200.0
  forward = 100.0
  tte = 1.0
  isCall = true
  intr = forward - strike
  if !isCall
    intr = -intr
  end
  vols =0.02:0.01:4.0
  mvols = zeros(length(vols))
  mprices = zeros(length(vols))
  prices = zeros(length(vols))

  for (j,vol) in enumerate(vols)
    price = Black.blackScholesFormula(isCall, strike, forward, vol*vol * tte, 1.0, 1.0)
    prices[j] = price
    mvols[j]= Black.impliedVolatilityJaeckel(isCall, price, forward, strike, tte, 1.0)
    mprices[j] = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
    # println(mprice-price, " ",price)
  end
  println("JAECKEL RMSD ",rmsd(mvols,vols)," MAXAD ",maxad(mvols,vols)," PriceMAXAD ",maxad(mprices,prices), " ",maxad(mprices ./ prices,prices ./ prices))
  @test isless(maxad(mvols,vols)/vols[end], 3e-15)
  for (j,vol) in enumerate(vols)
    price = prices[j]
    mvols[j] = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.Halley())
    mprices[j] = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
  end
  println("HALLEY RMSD ",rmsd(mvols,vols)," MAXAD ",maxad(mvols,vols)," PriceMAXAD ",maxad(mprices,prices), " ",maxad(mprices ./ prices,prices ./ prices))

  for (j,vol) in enumerate(vols)
    price = prices[j]
    mvols[j] = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.CMethod())
    mprices[j] = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
  end
  println("CMETHOD RMSD ",rmsd(mvols,vols)," MAXAD ",maxad(mvols,vols)," PriceMAXAD ",maxad(mprices,prices), " ",maxad(mprices ./prices,prices ./prices))

  for (j,vol) in enumerate(vols)
    price = prices[j]
    mvols[j] = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64,Black.Householder())
    mprices[j] = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
    # mprice = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
    # println(mprice-price, " ",price)
  end
  println("HOUSEHOLDER RMSD ",rmsd(mvols,vols)," MAXAD ",maxad(mvols,vols)," PriceMAXAD ",maxad(mprices,prices)," ",maxad(mprices ./ prices,prices ./ prices))

end

@testset "JaeckelAccuracyFixedVol" begin
  forward = 100.0
  tte = 1.0
  isCall = true
  vol =0.1
  strikes = 100.0:1.0:500.0
  mvols = zeros(length(strikes))
  mprices = zeros(length(strikes))
  prices = zeros(length(strikes))
  vols = ones(length(strikes))*vol
  for (j,strike) in enumerate(strikes)
    price = Black.blackScholesFormula(isCall, strike, forward, vol*vol * tte, 1.0, 1.0)
    prices[j] = price
    mvols[j]= Black.impliedVolatilityJaeckel(isCall, price, forward, strike, tte, 1.0)
    mprices[j] = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
    # println(mprice-price, " ",price)
  end
  println("JAECKEL RMSD ",rmsd(mvols,vols)," MAXAD ",maxad(mvols,vols)," PriceMAXAD ",maxad(mprices,prices), " ",maxad(mprices ./ prices,prices ./ prices))
  @test isless(maxad(mvols,vols), 4e-15)
  for (j,strike) in enumerate(strikes)
    price = prices[j]
    mvols[j] = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.Halley())
    mprices[j] = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
  end
  println("HALLEY RMSD ",rmsd(mvols,vols)," MAXAD ",maxad(mvols,vols)," PriceMAXAD ",maxad(mprices,prices), " ",maxad(mprices ./ prices,prices ./ prices))
  @test isless(maxad(mvols,vols), 4e-15)

  for (j,strike) in enumerate(strikes)
    price = prices[j]
    mvols[j] = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.CMethod())
    mprices[j] = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
  end
  println("CMETHOD RMSD ",rmsd(mvols,vols)," MAXAD ",maxad(mvols,vols)," PriceMAXAD ",maxad(mprices,prices), " ",maxad(mprices ./prices,prices ./prices))
  @test isless(maxad(mvols,vols), 4e-15)

  for (j,strike) in enumerate(strikes)
    price = prices[j]
    mvols[j] = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64,Black.Householder())
    mprices[j] = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
    # mprice = Black.blackScholesFormula(isCall, strike, forward, mvols[j]^2 * tte, 1.0, 1.0)
    # println(mprice-price, " ",price)
  end
  println("HOUSEHOLDER RMSD ",rmsd(mvols,vols)," MAXAD ",maxad(mvols,vols)," PriceMAXAD ",maxad(mprices,prices)," ",maxad(mprices ./ prices,prices ./ prices))
  @test isless(maxad(mvols,vols), 4e-15)

end

@testset "JaeckelSet6" begin
  price = 271.43234885190117
  strike = 275.0
  forward = 356.73063159822254
  tte = 1.5917808219178082
  isCall = false
  intr = forward - strike
  if !isCall
    intr = -intr
  end
  ivj = Black.impliedVolatilityJaeckel(isCall, price, forward, strike, tte, 1.0)
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println("Jaeckel ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-25)
  @test isapprox(e2, 0, atol = 1e-25)

  guess = 0.0
  ivj = Black.impliedVolatilityLiSOR(isCall, price, forward, strike, tte, 1.0, guess, 0e-14, 64, Black.SORTS()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SOR-Li ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.Halley())
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-Halley ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.SuperHalley()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-SHalley ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.CMethod()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-CMethod ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)
end

@testset "JaeckelSet5" begin
  price = 355.12907714987386
  forward = 2838.9998434067384
  strike = 2485.0
  tte = 0.08493150684931507

  isCall = true
  intr = forward - strike
  if !isCall
    intr = -intr
  end
  ivj = Black.impliedVolatilityJaeckel(isCall, price, forward, strike, tte, 1.0)
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-25)
  @test isapprox(e2, 0, atol = 1e-25)

  guess = 0.0
  ivj = Black.impliedVolatilityLiSOR(isCall, price, forward, strike, tte, 1.0, guess, 0e-14, 64, Black.SORTS())
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SOR-Li ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.Halley()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-Halley ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.SuperHalley()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-SHalley ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.CMethod()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-CMethod ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

end

@testset "JaeckelSet1" begin
  price = 0.35000000053458924
  forward = 1.0
  strike = 0.65
  tte = 0.0191780822

  isCall = true
  intr = forward - strike
  if !isCall
    intr = -intr
  end
  ivj = Black.impliedVolatilityJaeckel(isCall, price, forward, strike, tte, 1.0)
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-25)
  @test isapprox(e2, 0, atol = 1e-25)

  guess = 0.0
  ivj = Black.impliedVolatilityLiSOR(isCall, price, forward, strike, tte, 1.0, guess, 0e-14, 64, Black.SORTS())
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SOR-Li ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.Halley()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-Halley ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.SuperHalley()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-SHalley ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.CMethod()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-CMethod ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)
end

@testset "JaeckelSetATM" begin
  forward = 100.0
  strike = 100.0
  vol = 0.08
  tte = 1.0
  df = 1.0
  isCall = true
  price = Black.blackScholesFormula(isCall, strike, forward, vol * vol * tte, 1.0, df)
  intr = forward - strike
  if !isCall
    intr = -intr
  end
  ivj = Black.impliedVolatilityJaeckel(isCall, price, forward, strike, tte, 1.0)
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-25)
  @test isapprox(e2, 0, atol = 1e-25)
  strike = 101.0
  price = Black.blackScholesFormula(isCall, strike, forward, vol * vol * tte, 1.0, df)
  intr = forward - strike
  if !isCall
    intr = -intr
  end
  ivj = Black.impliedVolatilityJaeckel(isCall, price, forward, strike, tte, 1.0)
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-25)
  @test isapprox(e2, 0, atol = 1e-25)

  guess = 0.0
  ivj = Black.impliedVolatilityLiSOR(isCall, price, forward, strike, tte, 1.0, guess, 0e-14, 64, Black.SORTS())
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SOR-Li ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.Halley()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-Halley ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.SuperHalley()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-SHalley ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

  ivj = Black.impliedVolatilitySRHalley(isCall, price, forward, strike, tte, 1.0, 0e-14, 64, Black.CMethod()) #1e-14=> 0 but 0=> 1e-13 error; likely because blackfor erfc
  e1 = Black.blackScholesFormula(isCall, strike, forward, ivj * ivj * tte, 1.0, 1.0) - price
  e2 = Black.blackScholesFormula(isCall, forward, strike, ivj * ivj * tte, 1.0, 1.0) - price + intr
  println(" SR-CMethod ", ivj, " ", e1, " ", e2)
  @test isapprox(e1, 0, atol = 1e-12)
  @test isapprox(e2, 0, atol = 1e-12)

end


text:
# The set-up was designed for the
# 8th International Cloud Modelling Workshop
# ([Muhlbauer2013](@cite))
#
# See chapter 2 in [Arabas2015](@cite) for setup details:

using Dates
using DocStringExtensions
using LinearAlgebra
using Logging
using MPI
using Printf
using StaticArrays
using Test

using ClimateMachine
ClimateMachine.init(diagnostics = "default")
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.Diagnostics
using ClimateMachine.Grids
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Thermodynamics:
    gas_constants,
    PhaseEquil,
    PhasePartition_equil,
    PhasePartition,
    internal_energy,
    q_vap_saturation,
    relative_humidity,
    PhaseEquil_ρTq,
    PhaseNonEquil_ρTq,
    air_temperature,
    latent_heat_fusion,
    Liquid,
    Ice,
    supersaturation,
    vapor_specific_humidity

using ClimateMachine.Microphysics
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.VTK

using CLIMAParameters
using CLIMAParameters.Planet:
    R_d, cp_d, cv_d, cv_v, cv_l, cv_i, T_0, T_freeze, e_int_v0, e_int_i0, grav

using CLIMAParameters.Atmos.Microphysics

struct LiquidParameterSet <: AbstractLiquidParameterSet end
struct IceParameterSet <: AbstractIceParameterSet end
struct RainParameterSet <: AbstractRainParameterSet end
struct SnowParameterSet <: AbstractSnowParameterSet end
struct MicropysicsParameterSet{L, I, R, S} <: AbstractMicrophysicsParameterSet
    liquid::L
    ice::I
    rain::R
    snow::S
end
struct EarthParameterSet{M} <: AbstractEarthParameterSet
    microphys_param_set::M
end

const microphys_param_set = MicropysicsParameterSet(
    LiquidParameterSet(),
    IceParameterSet(),
    RainParameterSet(),
    SnowParameterSet(),
)
const param_set = EarthParameterSet(microphys_param_set)
const liquid_param_set = param_set.microphys_param_set.liquid
const ice_param_set = param_set.microphys_param_set.ice
const rain_param_set = param_set.microphys_param_set.rain
const snow_param_set = param_set.microphys_param_set.snow

using ClimateMachine.BalanceLaws:
    BalanceLaw, Prognostic, Auxiliary, Gradient, GradientFlux, Hyperdiffusive

import ClimateMachine.BalanceLaws:
    vars_state,
    init_state_prognostic!,
    init_state_auxiliary!,
    nodal_init_state_auxiliary!,
    nodal_update_auxiliary_state!,
    flux_first_order!,
    flux_second_order!,
    wavespeed,
    parameter_set,
    boundary_conditions,
    boundary_state!,
    source!

import ClimateMachine.DGMethods: DGModel
using ClimateMachine.Mesh.Geometry: LocalGeometry

struct KinematicModelConfig{FT}
    xmax::FT
    ymax::FT
    zmax::FT
    wmax::FT
    θ_0::FT
    p_0::FT
    p_1000::FT
    qt_0::FT
    z_0::FT
    periodicity_x::Bool
    periodicity_y::Bool
    periodicity_z::Bool
    idx_bc_left::Int
    idx_bc_right::Int
    idx_bc_front::Int
    idx_bc_back::Int
    idx_bc_bottom::Int
    idx_bc_top::Int
end

struct KinematicModel{FT, PS, O, M, P, S, BC, IS, DC} <: BalanceLaw
    param_set::PS
    orientation::O
    moisture::M
    precipitation::P
    source::S
    boundarycondition::BC
    init_state_prognostic::IS
    data_config::DC
end

parameter_set(m::KinematicModel) = m.param_set

function KinematicModel{FT}(
    ::Type{AtmosLESConfigType},
    param_set::AbstractParameterSet;
    orientation::O = FlatOrientation(),
    moisture::M = nothing,
    precipitation::P = nothing,
    source::S = nothing,
    boundarycondition::BC = nothing,
    init_state_prognostic::IS = nothing,
    data_config::DC = nothing,
) where {FT <: AbstractFloat, O, M, P, S, BC, IS, DC}

    @assert param_set ≠ nothing
    @assert init_state_prognostic ≠ nothing

    atmos = (
        param_set,
        orientation,
        moisture,
        precipitation,
        source,
        boundarycondition,
        init_state_prognostic,
        data_config,
    )

    return KinematicModel{FT, typeof.(atmos)...}(atmos...)
end

vars_state(m::KinematicModel, ::Gradient, FT) = @vars()

vars_state(m::KinematicModel, ::GradientFlux, FT) = @vars()

function nodal_init_state_auxiliary!(
    m::KinematicModel,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    FT = eltype(aux)
    x, y, z = geom.coord
    dc = m.data_config
    param_set = parameter_set(m)

    _R_d::FT = R_d(param_set)
    _cp_d::FT = cp_d(param_set)
    _grav::FT = grav(param_set)

    # TODO - should R_d and cp_d here be R_m and cp_m?
    R_m, cp_m, cv_m, γ = gas_constants(param_set, PhasePartition(dc.qt_0))

    # Pressure profile assuming hydrostatic and constant θ and qt profiles.
    # It is done this way to be consistent with Arabas paper.
    # It's not necessarily the best way to initialize with our model variables.
    p =
        dc.p_1000 *
        (
            (dc.p_0 / dc.p_1000)^(_R_d / _cp_d) -
            _R_d / _cp_d * _grav / dc.θ_0 / R_m * (z - dc.z_0)
        )^(_cp_d / _R_d)

    @inbounds begin
        aux.p = p
        aux.x_coord = x
        aux.z_coord = z
    end
end

function init_state_prognostic!(
    m::KinematicModel,
    state::Vars,
    aux::Vars,
    localgeo,
    t,
    args...,
)
    m.init_state_prognostic(m, state, aux, localgeo, t, args...)
end
boundary_conditions(::KinematicModel) = (1, 2, 3, 4, 5, 6)

function boundary_state!(
    ::CentralNumericalFluxSecondOrder,
    bctype,
    m::KinematicModel,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    t,
    args...,
) end

@inline function flux_second_order!(
    m::KinematicModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
) end

@inline function flux_second_order!(
    m::KinematicModel,
    flux::Grad,
    state::Vars,
    τ,
    d_h_tot,
) end

function config_kinematic_eddy(
    FT,
    N,
    resolution,
    xmax,
    ymax,
    zmax,
    wmax,
    θ_0,
    p_0,
    p_1000,
    qt_0,
    z_0,
    periodicity_x,
    periodicity_y,
    periodicity_z,
    idx_bc_left,
    idx_bc_right,
    idx_bc_front,
    idx_bc_back,
    idx_bc_bottom,
    idx_bc_top,
)
    # Choose explicit solver
    ode_solver = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    kmc = KinematicModelConfig(
        FT(xmax),
        FT(ymax),
        FT(zmax),
        FT(wmax),
        FT(θ_0),
        FT(p_0),
        FT(p_1000),
        FT(qt_0),
        FT(z_0),
        Bool(periodicity_x),
        Bool(periodicity_y),
        Bool(periodicity_z),
        Int(idx_bc_left),
        Int(idx_bc_right),
        Int(idx_bc_front),
        Int(idx_bc_back),
        Int(idx_bc_bottom),
        Int(idx_bc_top),
    )

    # Set up the model
    model = KinematicModel{FT}(
        AtmosLESConfigType,
        param_set;
        init_state_prognostic = init_kinematic_eddy!,
        data_config = kmc,
    )

    config = ClimateMachine.AtmosLESConfiguration(
        "KinematicModel",
        N,
        resolution,
        FT(xmax),
        FT(ymax),
        FT(zmax),
        param_set,
        init_kinematic_eddy!,
        boundary = (
            (Int(idx_bc_left), Int(idx_bc_right)),
            (Int(idx_bc_front), Int(idx_bc_back)),
            (Int(idx_bc_bottom), Int(idx_bc_top)),
        ),
        periodicity = (
            Bool(periodicity_x),
            Bool(periodicity_y),
            Bool(periodicity_z),
        ),
        xmin = FT(0),
        ymin = FT(0),
        zmin = FT(0),
        solver_type = ode_solver,
        model = model,
    )

    return config
end


text:
"""
    Periodic{P, S}(f, period::P, [phase_shift::S]) <: Transform

Applies a periodic function `f` with provided `period` and `phase_shift` to the data.

The `period` and `phase_shift` must have the same supertype of `Real` or `Period`, depending on
whether the data is `Real` or `TimeType` respectively.

!!! note
    For `TimeType` data, the result will change depending on the type of `period` given,
    even if the same amount of time is described. Example: `Week(1)` vs `Second(Week(1))`;
    the former starts the period on the most recent Monday, while the latter starts the
    period on the most recent multiple of 604800 seconds since time 0.

# Fields
* `f::Union{typeof(cos), typeof(sin)}`: the periodic function
* `period::Union{Real, Period}`: the function period. Must be strictly positive.
* `phase_shift::Union{Real, Period}` (optional): adjusts the phase of the periodic function,
  measured in the same units as the input. Increasing the value translates the function to
  the right, toward higher/later input values.
"""
struct Periodic{P, S} <: Transform
    f::Union{typeof(cos), typeof(sin)}
    period::P
    phase_shift::S

    function Periodic(f, period::P, phase_shift::S) where {P, S}
        if !((P <: Real && S <: Real) || (P <: Period && S <: Period))
            throw(ArgumentError("period and phase_shift must have the same supertype"))
        end
        period > zero(P) || throw(DomainError(period, "period must be strictly positive."))
        return new{P, S}(f, period, phase_shift)
    end
end

cardinality(::Periodic) = OneToOne()

"""
    Periodic(f, period) -> Periodic

A constructor for [`Periodic`](@ref).
Returns a `Periodic` transform with zero phase shift.
"""
Periodic(f, period::P) where P = Periodic(f, period, zero(P))

function _apply(x, P::Periodic; kwargs...)
    return P.f.(2π .* (x .- P.phase_shift) / P.period)
end

function _apply(x, p::Periodic{P}; kwargs...) where P <: Period
    map(xi -> _periodic(p.f, xi, p.period, p.phase_shift), x)
end

"""
    _periodic(f, instant, period, phase_shift=Day(0))

Computes the value of periodic function `f` at the given instant in time.

# Arguments
* `f`: the periodic function
* `period`: the function period
* `phase_shift`: adjusts the phase of the periodic function. Increasing the value translates
    the function to the right, toward higher/later input values.
"""
function _periodic(f, instant, period, phase_shift=Day(0))
    period_begin = floor(instant, period) + phase_shift

    # Make sure the `instant - period_begin < period` after we add the `phase_shift`.
    # Performing these adjustments isn't necessary since we're dealing with a periodic
    # function but it will reduce floating-point errors.
    if phase_shift > zero(typeof(period))
        while period_begin > instant
            period_begin -= period
        end
    elseif phase_shift < zero(typeof(period))
        while period_begin < instant
            period_begin += period  # Make sure `period_begin <= instant`
        end
    end

    period_end = period_begin + period
    # `period_end - period_begin` converts to the same units so that `/` is defined
    normalized = (instant - period_begin) / (period_end - period_begin)
    return f(2π * normalized)
end


text:
# Identity matrix
eye(m::Int64) = Diagonal(ones(m))

unpack(s, f::Symbol) = reduce(hcat, getproperty.(s, f))

function unpack(s, f::Vector{Symbol})
    A = []
    for i in f
        push!(A, reduce(hcat, getproperty.(s, i)))
    end
    return A
end

unpack_matrix(s, f) = getproperty.(s, f)

function unpack_matrix(s, f::Vector{Symbol})
    A = []
    for i in f
        push!(A, getproperty.(s, i))
    end
    return A
end

_return_type(P::Matrix, C) = Matrix(C)
_return_type(P, C) = Diagonal(C)

_square_root_filter(filter::CholeskyFactorization, P) = _return_type(P, cholesky(P).L)

function _square_root_filter(filter::ModifiedCholeskyFactorization, P)
	F = _modified_cholesky(P)
	_return_type(P, F.U * sqrt.(F.D))
end

function _square_root_filter(filter::SVDFactorization, P)
	_C = svd(Matrix(P))
	return _return_type(P, _C.U * sqrt(Diagonal(_C.S)))
end

"""

	F = _modified_cholesky(P)

Computes modified Cholesky factors F.U and F.D of a symmetric positive definite matrix P,
such that F.U is unit upper triangular and F.D is diagonal. Thus, and P = F.U*F.D*F.U'.

Reference: Grewal - Kalman filtering theory and practice

"""
function _modified_cholesky(P::AbstractMatrix)
	isposdef(P) || throw("input matrix must be positive-definite. Modified Cholesky factorization failed.")
	n = size(P, 1)
	_P = (P .+ P') ./ 2.0 # take symmetric part
	U, D = similar(_P), similar(_P)
	for j in n:-1:1
   		for i in j:-1:1
			s = P[i, j]
			for k in j+1:n
				s -= U[i, k]*D[k, k]*U[j, k]
			end
			if i==j
				D[j, j] = s
				U[j, j] = 1.0
			else
				U[i, j] = s / D[j, j]
			end
		end
	end
	return (U=U, D=D)
end


text:
using LowRankModels
using Test, Random, SparseArrays

using LinearAlgebra: norm

###############################################################
# verify basic functionality works
include("basic_functionality.jl")

###############################################################
# verify most losses, regularizers, and ways of calling glrm work
include("hello_world.jl")

################################################################################
# ScikitLearnBase test
import ScikitLearnBase

# Check that KMeans can correctly separate two non-overlapping Gaussians
Random.seed!(21)
gaussian1 = randn(100, 2) .+ 5.
gaussian2 = randn(50, 2) .- 10.
A = vcat(gaussian1, gaussian2)

model = ScikitLearnBase.fit!(LowRankModels.KMeans(), A)
@test Set(sum(ScikitLearnBase.transform(model, A), dims=1)) == Set([100, 50])

###############################################################
# test sparse and streaming functionality
include("sparse_test.jl")
include("streaming_test.jl")

###############################################################
# verify all examples run
include("../examples/runexamples.jl")

###############################################################
# verify all tests of specific loss functions run
include("prob_tests/runtests.jl")


text:
using ARDESPOT
using Random

mt = MersenneTwister(1)

rng = MemorizingRNG(mt)

rand(rng)
rand(rng, Float64)
rand(rng, Int32)
rand(rng, [:a, :b, :c])
rand(rng, Int)
randn(rng)
rand(rng, 1:1_000_000)
randn(rng, (3,3))

@test isapprox(sum(rand(rng, Bool, 1_000_000))/1_000_000, 0.5, atol=0.001)
@test isapprox(sum(rand(rng, [1, 1, 1, 1, 0, 0, 0, 0], 1_000_000))/1_000_000, 0.5, atol=0.001)


text:
function showcoeff(io::IO, n::Number, show_sign::Bool, show_one::Bool=false)
    isimag(n) = !iszero(n) && isreal(im*n)
    isneg(n) = n < 0

    if isreal(n) && isneg(real(n)) || isimag(n) && isneg(imag(n))
        write(io, "- ")
    else
        show_sign && write(io, "+ ")
    end
    n = round(n, digits=6)
    if !(isone(n) || isone(-n))
        if !(isreal(n) || isimag(n))
            write(io, "($n)")
        elseif isimag(n)
            ii = imag(n)
            if !(isone(ii) || isone(-ii))
                write(io, "$(abs(ii))im")
            else
                write(io, "im")
            end
        else
            write(io, string(abs(n)))
        end
    elseif show_one
        write(io, "1")
    end
end


text:
using Turing
using Distributions
using Base.Test

# Test the @assume and @predict macros on a model without conditioning.

@model test begin
  @assume x ~ Bernoulli(1)
  @assume y ~ Bernoulli(x / 2)
  @predict y
  @predict x
end

s = SMC(10000)
p = PG(10,1000)

res = sample(test, s)

@test reduce(&, res[:x]) == 1  #check that x is always 1
@test_approx_eq_eps mean(res[:y]) 0.5 0.1  #check that the mean of y is between 0.4 and 0.6

res = sample(test, p)

@test reduce(&, res[:x]) == 1  #check that x is always 1
@test_approx_eq_eps mean(res[:y]) 0.5 0.1  #check that the mean of y is between 0.4 and 0.6


text:
@doc """
	ConstGate <: QGate
# Structure
- `kernel::GateKernel` : See also **GateKernel**.
""" ->
ConstGate

@doc """
	ConstGate(data::Array{T} where {T<:Number}, name::String="Anonymous")
Construct a ConstGate.

# Arguments
- `data`: A tensor representation of the gate.
- `name`: Name string of the gate; used in print functions. Default value: "Anonymous".
# Example
```Julia
julia> xdata = [0 1; 1 0]
2×2 Array{Int64,2}:
 0  1
 1  0

julia> Xgate = ConstGate(xdata, "self defined X")
self defined X
```
""" ->
ConstGate(data::Array{T} where {T<:Number}, name::String="Anonymous")

@doc """
	ConstGate(f::Function, name::String="Anonymous")
Construct a ConstGate with given function.
# Arguments
- `f`: A function that returns gate tensor.
- `name`: Name string used in print functions.
""" ->
ConstGate(f::Function, name::String="Anonymous")

# @doc """
# 	name(gate::ConstGate)
# Return name of the gate.

# """ ->
# name(gate::ConstGate)

# @doc """
# 	data(gate::ConstGate)
# Return gate tensor.
# """ ->
# data(gate::ConstGate)


text:
"""
Parse files provided by the Stanford Network Analysis Project (SNAP).
"""
module Snap

import CodecZlib: GzipDecompressorStream
using ..Io: parse_uints
using ...GraphIdx: BiAdjacentIndex

"""
    open_snap(f, fname; as_process = false)

Open file `fname` to perform reading with function `f`.
If `fname` ends with `".bz2"`, start a corresponding decompressor.
If `fname == "-"`, read from `stdin`.
"""
function open_snap(f::Function, fname::String; as_process::Bool=false)
    if endswith(fname, ".gz")
        if as_process
            open(`gunzip -d -c $fname`) do io
                return f(io)
            end
        else
            open(fname) do gz
                return f(GzipDecompressorStream(gz))
            end
        end
    else
        if fname == "-"
            f(stdin)
        else
            @assert ispath(fname)
            open(f, fname, "r")
        end
    end
end

"""
    parse_snap_header(io)

Read lines until lines does not start with comment character "%" and
then return the first two integers (number of nodes, number of edges).
"""
function parse_snap_header(io::IO)::Tuple{Int,Int}
    line = ""
    n, m = -1, -1
    while startswith(line, "#") || line == ""
        if (g = match(r"# Nodes: (\d+) Edges: (\d+)", line)) != nothing
            n, m = [parse(Int, g[i]) for i = 1:2]
        end
        mark(io)
        line = readline(io)
    end
    reset(io)
    unmark(io)
    n, m
end

"""
    parse_snap_idx(io)

Parse a Snap file into a bidirectional index.
"""
function parse_snap_idx(io::IO)::BiAdjacentIndex
    n, m, head, tail = parse_snap_edges(io, Int)
    return BiAdjacentIndex(n, head, tail)
end

"""
    parse_snap_edges(io)

Return the graph as `n, m, head, tail`.
"""

function parse_snap_edges(io::IO, ::Type{T} = Int32) where {T<:Integer}
    n, m = parse_snap_header(io)
    head = Vector{T}()
    tail = Vector{T}()
    sizehint!(head, m)
    sizehint!(tail, m)
    parse_uints(io) do u::UInt, last::Bool
        if last
            push!(tail, T(u))
        else
            push!(head, T(u))
        end
    end
    return n, m, head, tail
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
abstract type AbstractCollisionDynamics end

struct ElasticCollision <: AbstractCollisionDynamics end

"Elastic collision of ball with FixedPlane.
The ball is assumed to be touching the FixedPlane."
function collide!(b::MovableBall{N,T}, Π::FixedPlane{N,T}, ::ElasticCollision) where {N,T}
    v = b.v
    n = Π.n

    b.v = v - 2 * (v⋅n) * n
end

"Assumes b1 and b2 are touching"
function collide!(b1::MovableBall, b2::MovableBall, ::ElasticCollision)
    Δx = b1.x - b2.x

    # @show norm(Δx)

    v1 = b1.v
    v2 = b2.v

    m1 = b1.m
    m2 = b2.m

    # reference: https://en.wikipedia.org/wiki/Elastic_collision#Two-dimensional_collision_with_two_moving_objects

    # Component orthogonal to line joining centers
    # is unchanged by elastic collision (since impulse is applied in direction joining sphere centers)

    n = normalize(Δx)   # vector joining sphere centers

    v1n = (v1 ⋅ n) * n  # velocity components in that direction
    v2n = (v2 ⋅ n) * n

    # treat terms m1 / (m1 + m2)  and m2 / (m1 + m2)
    # by dividing top and bottom by m1
    # This allows one of them to be Inf

    mass_ratio = (1 / ((m1 / m2) + 1))

    v1′ = v1 + 2 * mass_ratio * (v2n - v1n)  # update only components in n direction
    v2′ = v2 + 2 * (1 - mass_ratio) * (v1n - v2n)

    b1.v = v1′
    b2.v = v2′

end

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
function ConjugateGradient(arg0::IterationManager, arg1::jdouble, arg2::jboolean)
    return ConjugateGradient((IterationManager, jdouble, jboolean), arg0, arg1, arg2)
end

function ConjugateGradient(arg0::jint, arg1::jdouble, arg2::jboolean)
    return ConjugateGradient((jint, jdouble, jboolean), arg0, arg1, arg2)
end

function get_iteration_manager(obj::IterativeLinearSolver)
    return jcall(obj, "getIterationManager", IterationManager, ())
end

function should_check(obj::ConjugateGradient)
    return jcall(obj, "shouldCheck", jboolean, ())
end

function solve(obj::PreconditionedIterativeLinearSolver, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealLinearOperator, RealVector), arg0, arg1, arg2)
end

function solve(obj::PreconditionedIterativeLinearSolver, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector, arg3::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2, arg3)
end

function solve(obj::PreconditionedIterativeLinearSolver, arg0::RealLinearOperator, arg1::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealVector), arg0, arg1)
end

function solve(obj::PreconditionedIterativeLinearSolver, arg0::RealLinearOperator, arg1::RealVector, arg2::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2)
end

function solve_in_place(obj::ConjugateGradient, arg0::RealLinearOperator, arg1::RealLinearOperator, arg2::RealVector, arg3::RealVector)
    return jcall(obj, "solveInPlace", RealVector, (RealLinearOperator, RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2, arg3)
end

function solve_in_place(obj::PreconditionedIterativeLinearSolver, arg0::RealLinearOperator, arg1::RealVector, arg2::RealVector)
    return jcall(obj, "solveInPlace", RealVector, (RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2)
end