text:
using TOML
using DocStringExtensions

export AbstractTOMLDict
export ParamDict, AliasParamDict

export float_type,
    get_parameter_values!,
    get_parameter_values,
    write_log_file,
    log_parameter_information,
    create_toml_dict

"""
    AbstractTOMLDict{FT <: AbstractFloat}

Abstract parameter dict. Two subtypes:
 - [`ParamDict`](@ref)
 - [`AliasParamDict`](@ref)
"""
abstract type AbstractTOMLDict{FT <: AbstractFloat} end

const NAMESTYPE =
    Union{AbstractVector{S}, NTuple{N, S} where {N}} where {S <: AbstractString}

"""
    ParamDict(data::Dict, override_dict::Union{Nothing,Dict})

Structure to hold information read-in from TOML
file, as well as a parametrization type `FT`.

Uses the name to search

# Fields

$(DocStringExtensions.FIELDS)
"""
struct ParamDict{FT} <: AbstractTOMLDict{FT}
    "dictionary representing a default/merged parameter TOML file"
    data::Dict
    "either a nothing, or a dictionary representing an override parameter TOML file"
    override_dict::Union{Nothing, Dict}
end

"""
    AliasParamDict(data::Dict, override_dict::Union{Nothing,Dict})

Structure to hold information read-in from TOML
file, as well as a parametrization type `FT`.

Uses the alias to search

# Fields

$(DocStringExtensions.FIELDS)
"""
struct AliasParamDict{FT} <: AbstractTOMLDict{FT}
    "dictionary representing a default/merged parameter TOML file"
    data::Dict
    "either a nothing, or a dictionary representing an override parameter TOML file"
    override_dict::Union{Nothing, Dict}
    "Alias->name map"
    alias_to_name_map::Dict
    function AliasParamDict{FT}(
        data::Dict,
        override_dict::Union{Nothing, Dict},
    ) where {FT}
        alias_to_name_map = Dict(map(collect(keys(data))) do key
            Pair(data[key]["alias"], key)
        end)
        return new{FT}(data, override_dict, alias_to_name_map)
    end
end

"""
    float_type(::AbstractTOMLDict)

The float type from the parameter dict.
"""
float_type(::AbstractTOMLDict{FT}) where {FT} = FT

function Base.iterate(pd::AliasParamDict)
    it = iterate(pd.data)
    if it !== nothing
        return (Pair(it[1].second["alias"], it[1].second), it[2])
    else
        return nothing
    end
end

function Base.iterate(pd::AliasParamDict, state)
    it = iterate(pd.data, state)
    if it !== nothing
        return (Pair(it[1].second["alias"], it[1].second), it[2])
    else
        return nothing
    end
end

Base.iterate(pd::ParamDict, state) = Base.iterate(pd.data, state)
Base.iterate(pd::ParamDict) = Base.iterate(pd.data)

Base.getindex(pd::ParamDict, i) = getindex(pd.data, i)
Base.getindex(pd::AliasParamDict, i) =
    getindex(pd.data, pd.alias_to_name_map[i])

"""
    log_component!(pd::AbstractTOMLDict, names, component)

Adds a new key,val pair: `("used_in",component)` to each
named parameter in `pd`.
Appends a new val: `component` if "used_in" key exists.
"""
function log_component!(
    pd::AliasParamDict,
    names::NAMESTYPE,
    component::AbstractString,
)
    component_key = "used_in"
    data = pd.data
    for name in names
        for (key, val) in data
            name ≠ val["alias"] && continue
            data[key][component_key] = if component_key in keys(data[key])
                unique([data[key][component_key]..., component])
            else
                [component]
            end
        end
    end
end

function log_component!(
    pd::ParamDict,
    names::NAMESTYPE,
    component::AbstractString,
)
    component_key = "used_in"
    data = pd.data
    for name in names
        for (key, val) in data
            name ≠ key && continue
            data[key][component_key] = if component_key in keys(data[key])
                unique([data[key][component_key]..., component])
            else
                [component]
            end
        end
    end
end

"""
    _get_typed_value(pd, val, valname, valtype)

enforces `val` to be of type as specified in the toml file
- `float_type(pd)` if type=\"float\"
- `Int` if type=\"integer\"
- `String` if  type=\"string\"
Default type of `String` is used if no type is provided.
"""
function _get_typed_value(
    pd::AbstractTOMLDict,
    val,
    valname::AbstractString,
    valtype,
)

    if valtype == "float"
        return float_type(pd)(val)
    elseif valtype == "integer"
        return Int(val)
    elseif valtype == "string"
        return String(val)
    else
        error(
            "For parameter with identifier: \"",
            valname,
            "\", the attribute: type = \"",
            valtype,
            "\", is not recognised, ",
            "\n please select from: type = \"string\", \"float\", or \"integer\" ",
        )
    end
end

"""
    get_values(pd::AbstractTOMLDict, names)

gets the `value` of the named parameters.
"""
function get_values(pd::AliasParamDict, aliases::NAMESTYPE)
    data = pd.data
    # TODO: use map
    ret_values = []
    for alias in aliases
        for (key, val) in data
            alias ≠ val["alias"] && continue
            param_value = val["value"]
            param_type = haskey(val, "type") ? val["type"] : "string"
            elem =
                isa(param_value, AbstractVector) ?
                map(
                    x -> _get_typed_value(pd, x, alias, param_type),
                    param_value,
                ) : _get_typed_value(pd, param_value, alias, param_type)
            push!(ret_values, Pair(Symbol(alias), elem))
        end
    end
    return ret_values
end

function get_values(pd::ParamDict, names::NAMESTYPE)
    data = pd.data
    ret_values = map(names) do name
        param_value = data[name]["value"]
        param_type =
            haskey(data[name], "type") ? data[name]["type"] : "string"
        elem =
            isa(param_value, AbstractVector) ?
            map(x -> _get_typed_value(pd, x, name, param_type), param_value) : _get_typed_value(pd, param_value, name, param_type)
        Pair(Symbol(name), elem)
    end
    return ret_values
end

"""
    get_parameter_values!(
        pd::AbstractTOMLDict,
        names::Union{String,Vector{String}},
        component::String
    )

(Note the `!`) Gets the parameter values, and logs
the component (if given) where parameters are used.
"""
function get_parameter_values!(
    pd::AbstractTOMLDict,
    names::NAMESTYPE,
    component::Union{AbstractString, Nothing} = nothing,
)
    if !isnothing(component)
        log_component!(pd, names, component)
    end
    return get_values(pd, names)
end

get_parameter_values!(
    pd::AbstractTOMLDict,
    names::AbstractString,
    args...;
    kwargs...,
) = first(get_parameter_values!(pd, [names], args..., kwargs...))

"""
    get_parameter_values(pd::AbstractTOMLDict, names)

Gets the parameter values only.
"""
get_parameter_values(
    pd::AbstractTOMLDict,
    names::Union{NAMESTYPE, AbstractString},
) = get_parameter_values!(pd, names, nothing)

"""
    check_override_parameter_usage(pd::ParamDict, strict)

Checks if parameters in the ParamDict.override_dict have the
key "used_in" (i.e. were these parameters used within the model run).
Throws warnings in each where parameters are not used. Also throws
an error if `strict == true` .
"""
check_override_parameter_usage(pd::AbstractTOMLDict, strict::Bool) =
    check_override_parameter_usage(pd, strict, pd.override_dict)

check_override_parameter_usage(pd::AbstractTOMLDict, strict::Bool, ::Nothing) =
    nothing

function check_override_parameter_usage(
    pd::AbstractTOMLDict,
    strict::Bool,
    override_dict,
)
    unused_override = Dict()
    for (key, val) in override_dict
        logged_val = pd.data[key]
        unused_override[key] = !("used_in" in keys(logged_val))
    end
    if any(values(unused_override))
        unused_override_keys = collect(keys(unused_override))
        filter!(key -> unused_override[key], unused_override_keys)
        @warn(
            string(
                "Keys are present in parameter file but not used",
                "in the simulation. \n Typically this is due to",
                "a mismatch in parameter name in toml and in source.",
                "Offending keys: $(unused_override_keys)",
            )
        )
        if strict
            @error(
                "At least one override parameter set and not used in simulation"
            )
            error(
                "Halting simulation due to unused parameters." *
                "\n Typically this is due to a typo in the parameter name." *
                "\n change `strict` flag to `true` to prevent this causing an exception",
            )
        end
    end
    return nothing
end

"""
    write_log_file(pd::AbstractTOMLDict, filepath)

Writes a log file of all used parameters of `pd` at
the `filepath`. This file can be used to rerun the experiment.
"""
function write_log_file(pd::AbstractTOMLDict, filepath::AbstractString)
    used_parameters = Dict()
    for (key, val) in pd.data
        if "used_in" in keys(val)
            used_parameters[key] = val
        end
    end
    open(filepath, "w") do io
        TOML.print(io, used_parameters)
    end
end

"""
    log_parameter_information(
        pd::AbstractTOMLDict,
        filepath;
        strict::Bool = false
    )

Writes the parameter log file at `filepath`; checks that
override parameters are all used.

If `strict = true`, errors if override parameters are unused.
"""
function log_parameter_information(
    pd::AbstractTOMLDict,
    filepath::AbstractString;
    strict::Bool = false,
)
    #[1.] write the parameters to log file
    write_log_file(pd, filepath)
    #[2.] send warnings or errors if parameters were not used
    check_override_parameter_usage(pd, strict)
end

"""
    merge_override_default_values(
        override_toml_dict::AbstractTOMLDict{FT},
        default_toml_dict::AbstractTOMLDict{FT}
    ) where {FT}

Combines the `default_toml_dict` with the `override_toml_dict`,
precedence is given to override information.
"""
function merge_override_default_values(
    override_toml_dict::PDT,
    default_toml_dict::PDT,
) where {FT, PDT <: AbstractTOMLDict{FT}}
    data = default_toml_dict.data
    override_dict = override_toml_dict.override_dict
    for (key, val) in override_toml_dict.data
        if !(key in keys(data))
            data[key] = val
        else
            for (kkey, vval) in val # as val is a Dict too
                data[key][kkey] = vval
            end
        end
    end
    return PDT(data, override_dict)
end

"""
    create_toml_dict(FT;
        override_file,
        default_file,
        dict_type="alias"
    )

Creates a `ParamDict{FT}` struct, by reading and merging upto
two TOML files with override information taking precedence over
default information.
"""
function create_toml_dict(
    ::Type{FT};
    override_file::Union{Nothing, String} = nothing,
    default_file::String = joinpath(@__DIR__, "parameters.toml"),
    dict_type = "alias",
) where {FT <: AbstractFloat}
    @assert dict_type in ("alias", "name")
    PDT = _toml_dict(dict_type, FT)
    if isnothing(override_file)
        return PDT(TOML.parsefile(default_file), nothing)
    end
    override_toml_dict =
        PDT(TOML.parsefile(override_file), TOML.parsefile(override_file))
    default_toml_dict = PDT(TOML.parsefile(default_file), nothing)

    #overrides the defaults where they clash
    return merge_override_default_values(override_toml_dict, default_toml_dict)
end

function _toml_dict(s::String, ::Type{FT}) where {FT}
    if s == "alias"
        return AliasParamDict{FT}
    elseif s == "name"
        return ParamDict{FT}
    else
        error("Bad string option given")
    end
end

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
push!(LOAD_PATH, "./")

using Printf
using Statistics
using stiffness

const S = StiffnessData()
const x = Float64[0., 1., 0.]
const y = Float64[0., 0., 1.]
m = zeros(21)

nb_runs = 200

times = zeros(nb_runs)

for irun in 1:nb_runs
    times[irun] = @elapsed stiffness.op!(S, x, y, m)
end

@printf("%.3f µs\n", median(times) * 1e6)

ext: jl


text:
function get_ticks(f::Function, values; l = length(values))
    idxs = unique(floor.(Int, range(1, length(values); length = l)))
    values = values[idxs]
    values, f.(values)
end


text:
function FootprintOverlapDetector(arg0::FieldOfView, arg1::OneAxisEllipsoid, arg2::SphericalPolygonsSet, arg3::jdouble)
    return FootprintOverlapDetector((FieldOfView, OneAxisEllipsoid, SphericalPolygonsSet, jdouble), arg0, arg1, arg2, arg3)
end

function g(obj::FootprintOverlapDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_body(obj::FootprintOverlapDetector)
    return jcall(obj, "getBody", BodyShape, ())
end

function get_fov(obj::FootprintOverlapDetector)
    return jcall(obj, "getFOV", FieldOfView, ())
end

function get_zone(obj::FootprintOverlapDetector)
    return jcall(obj, "getZone", SphericalPolygonsSet, ())
end


text:
cd(@__DIR__)
using DifferentialEvolutionMCMC, Random, Parameters, Distributions
Random.seed!(88484)

prior_loglike(θ) = logpdf(Beta(1, 1), θ)

sample_prior() = rand(Beta(1, 1))

bounds = ((0,1),)
names = (:θ,)

N = 10
k = rand(Binomial(N, .5))
data = (N = N,k = k)

function loglike(data, θ)
    @unpack N,k = data
    n_sim = 10^4
    counter(_) = rand(Binomial(N, θ)) == k ? 1 : 0
    cnt = mapreduce(counter, +, 1:n_sim)
    return log(cnt / n_sim)
end

# loglike(θ, data) = logpdf(Binomial(data.N, θ), data.k)

model = DEModel(;
    sample_prior,
    prior_loglike,
    loglike,
    data,
    names
)

de = DE(;sample_prior, bounds, burnin=1000, Np=3, σ=.01)
n_iter = 2000
chains = sample(model, de, MCMCThreads(), n_iter, progress=true)


text:
"""
Generate a physically valid random density integrating to the given number of electrons.
"""
function random_density(basis::PlaneWaveBasis; n_electrons=basis.model.n_electrons)
    T = eltype(basis)
    ρtot  = rand(T, basis.fft_size)
    ρtot  = ρtot .* n_electrons ./ (sum(ρtot) * basis.dvol)  # Integration to n_electrons
    ρspin = nothing
    if basis.model.n_spin_components > 1
        ρspin = rand((-1, 1), basis.fft_size ) .* rand(T, basis.fft_size) .* ρtot
        @assert all(abs.(ρspin) .≤ ρtot)
    end
    ρ_from_total_and_spin(ρtot, ρspin)
end

@doc raw"""
    guess_density(basis, magnetic_moments)

Build a superposition of atomic densities (SAD) guess density.

We take for the guess density a gaussian centered around the atom, of
length specified by `atom_decay_length`, normalized to get the right number of electrons
```math
\hat{ρ}(G) = Z \exp\left(-(2π \text{length} |G|)^2\right)

When magnetic moments are provided, construct a symmetry-broken density guess.
The magnetic moments should be specified in units of ``μ_B``.
"""
function guess_density(basis::PlaneWaveBasis, magnetic_moments=[])
    guess_density(basis, basis.model.atoms, magnetic_moments)
end
@timing function guess_density(basis::PlaneWaveBasis{T}, atoms, magnetic_moments) where {T}
    ρtot = _guess_total_density(basis, atoms)
    if basis.model.n_spin_components == 1
        ρspin = nothing
    else
        ρspin = _guess_spin_density(basis, atoms, magnetic_moments)
    end

    ρ_from_total_and_spin(ρtot, ρspin)
end

function _guess_total_density(basis::PlaneWaveBasis{T}, atoms) where {T}
    # build ρtot
    gaussians_tot = [(T(n_elec_valence(spec)), T(atom_decay_length(spec)), pos)
                     for (spec, positions) in atoms for pos in positions]
    ρtot = gaussian_superposition(basis, gaussians_tot)
end

function _guess_spin_density(basis::PlaneWaveBasis{T}, atoms, magnetic_moments) where {T}
    model = basis.model
    if model.spin_polarization in (:none, :spinless)
        isempty(magnetic_moments) && return nothing
        error("Initial magnetic moments can only be used with collinear models.")
    end

    # If no magnetic moments start with a zero spin density
    all_magmoms = (normalize_magnetic_moment(magmom) for (_, magmoms) in magnetic_moments
                   for magmom in magmoms)
    if all(iszero, all_magmoms)
        @warn("Returning zero spin density guess, because no initial magnetization has " *
              "been specified in any of the given elements / atoms. Your SCF will likely " *
              "not converge to a spin-broken solution.")
        return zeros(T, basis.fft_size)
    end

    gaussians = Tuple{T, T, Vec3{T}}[]
    @assert length(magnetic_moments) == length(atoms)
    for (ispec, (spec, magmoms)) in enumerate(magnetic_moments)
        positions = atoms[ispec][2]
        @assert charge_nuclear(spec) == charge_nuclear(atoms[ispec][1])
        @assert length(magmoms) == length(positions)
        for (ipos, r) in enumerate(positions)
            magmom = Vec3{T}(normalize_magnetic_moment(magmoms[ipos]))
            iszero(magmom) && continue
            iszero(magmom[1:2]) || error("Non-collinear magnetization not yet implemented")

            magmom[3] ≤ n_elec_valence(spec) || error(
                "Magnetic moment $(magmom[3]) too large for element $(spec.symbol) with " *
                "only $(n_elec_valence(spec)) valence electrons."
            )
            push!(gaussians, (magmom[3], atom_decay_length(spec), r))
        end
    end
    gaussian_superposition(basis, gaussians)
end

@doc raw"""
Build a superposition of Gaussians as a guess for the density and magnetisation.
Expects a list of tuples `(coefficient, length, position)` for each of the Gaussian,
which follow the functional form
```math
\hat{ρ}(G) = \text{coefficient} \exp\left(-(2π \text{length} |G|)^2\right)
```
and are placed at `position` (in fractional coordinates).
"""
function gaussian_superposition(basis::PlaneWaveBasis{T}, gaussians) where {T}
    ρ = zeros(complex(T), basis.fft_size)
    isempty(gaussians) && return G_to_r(basis, ρ)

    # Fill ρ with the (unnormalized) Fourier transform, i.e. ∫ e^{-iGx} f(x) dx,
    # where f(x) is a weighted gaussian
    #
    # is formed from a superposition of atomic densities, each scaled by a prefactor
    for (iG, G) in enumerate(G_vectors(basis))
        Gsq = sum(abs2, basis.model.recip_lattice * G)
        for (coeff, decay_length, r) in gaussians
            form_factor::T = exp(-Gsq * T(decay_length)^2)
            ρ[iG] += T(coeff) * form_factor * cis(-2T(π) * dot(G, r))
        end
    end

    # projection in the normalized plane wave basis
    G_to_r(basis, ρ / sqrt(basis.model.unit_cell_volume))
end

@doc raw"""
Get the lengthscale of the valence density for an atom with `n_elec_core` core
and `n_elec_valence` valence electrons.
```
"""
function atom_decay_length(n_elec_core, n_elec_valence)
    # Adapted from ABINIT/src/32_util/m_atomdata.F90,
    # from which also the data has been taken.

    n_elec_valence = round(Int, n_elec_valence)
    if n_elec_valence == 0
        return 0.0
    end

    data = if n_elec_core < 0.5
        # Bare ions: Adjusted on 1H and 2He only
        [0.6, 0.4, 0.3, 0.25, 0.2]
    elseif n_elec_core < 2.5
        # 1s2 core: Adjusted on 3Li, 6C, 7N, and 8O
        [1.8, 1.4, 1.0, 0.7, 0.6, 0.5, 0.4, 0.35, 0.3]
    elseif n_elec_core < 10.5
        # Ne core (1s2 2s2 2p6): Adjusted on 11na, 13al, 14si and 17cl
        [2.0, 1.6, 1.25, 1.1, 1.0, 0.9, 0.8, 0.7 , 0.7, 0.7, 0.6]
    elseif n_elec_core < 12.5
        # Mg core (1s2 2s2 2p6 3s2): Adjusted on 19k, and on n_elec_core==10
        [1.9, 1.5, 1.15, 1.0, 0.9, 0.8, 0.7, 0.6 , 0.6, 0.6, 0.5]
    elseif n_elec_core < 18.5
        # Ar core (Ne + 3s2 3p6): Adjusted on 20ca, 25mn and 30zn
        [2.0, 1.8, 1.5, 1.2, 1.0, 0.9, 0.85, 0.8, 0.75, 0.7, 0.65, 0.65, 0.6]
    elseif n_elec_core < 28.5
        # Full 3rd shell core (Ar + 3d10): Adjusted on 31ga, 34se and 38sr
        [1.5, 1.25, 1.15, 1.05, 1.00, 0.95, 0.95, 0.9, 0.9, 0.85, 0.85, 0.80,
         0.8 , 0.75, 0.7]
    elseif n_elec_core < 36.5
        # Krypton core (Ar + 3d10 4s2 4p6): Adjusted on 39y, 42mo and 48cd
        [2.0, 2.00, 1.60, 1.40, 1.25, 1.10, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.7]
    else
        # For the remaining elements, consider a function of n_elec_valence only
        [2.0 , 2.00, 1.55, 1.25, 1.15, 1.10, 1.05, 1.0 , 0.95, 0.9, 0.85, 0.85, 0.8]
    end
    data[min(n_elec_valence, length(data))]
end
atom_decay_length(sp::Element) = atom_decay_length(n_elec_core(sp), n_elec_valence(sp))


text:
using Test
using LinearMaps
using SparseArrays
using LinearAlgebra

# adopted from: https://discourse.julialang.org/t/way-to-return-the-number-of-allocations/5167/10
macro numalloc(expr)
    return quote
        let
            local f
            function f()
                n1 = Base.gc_num()
                $(expr)
                n2 = Base.gc_num()
                diff = Base.GC_Diff(n2, n1)
                Base.gc_alloc_count(diff)
            end
            f()
        end
    end
end

A = 2 * rand(ComplexF64, (20, 10)) .- 1
v = rand(ComplexF64, 10)
w = rand(ComplexF64, 20)
V = rand(ComplexF64, 10, 3)
W = rand(ComplexF64, 20, 3)
α = rand()
β = rand()
M = LinearMap(A)
N = LinearMap(M)

@testset "LinearMaps.jl" begin
    @test eltype(M) == eltype(A)
    @test size(M) == size(A)
    @test size(N) == size(A)
    @test !isreal(M)
    @test ndims(M) == 2
    @test_throws ErrorException size(M, 3)
    @test length(M) == length(A)
    # matrix generation/conversion
    @test Matrix(M) == A
    @test Array(M) == A
    @test convert(Matrix, M) == A
    @test convert(Array, M) == A
    @test Matrix(M') == A'
    @test Matrix(transpose(M)) == copy(transpose(A))
    # sparse matrix generation/conversion
    @test sparse(M) == sparse(Array(M))
    @test convert(SparseMatrixCSC, M) == sparse(Array(M))

    B = copy(A)
    B[rand(1:length(A), 30)] .= 0
    MS = LinearMap(B)
    @test sparse(MS) == sparse(Array(MS)) == sparse(B)
end

Av = A * v
AV = A * V
@testset "mul! and *" begin
    @test M * v == Av
    @test N * v == Av
    @test @inferred mul!(copy(w), M, v) == mul!(copy(w), A, v)
    @test ((@allocated mul!(w, M, v)) == 0)
    @test @inferred mul!(copy(w), N, v) == Av

    # mat-vec-mul
    @test @inferred mul!(copy(w), M, v, 0, 0) == zero(w)
    @test @inferred mul!(copy(w), M, v, 0, 1) == w
    @test @inferred mul!(copy(w), M, v, 0, β) == β * w
    @test @inferred mul!(copy(w), M, v, 1, 1) ≈ Av + w
    @test @inferred mul!(copy(w), M, v, 1, β) ≈ Av + β * w
    @test @inferred mul!(copy(w), M, v, α, 1) ≈ α * Av + w
    @test @inferred mul!(copy(w), M, v, α, β) ≈ α * Av + β * w
    @test @inferred mul!(copy(w), M, v, α)    ≈ α * Av

    # test mat-mat-mul!
    @test @inferred mul!(copy(W), M, V, α, β) ≈ α * AV + β * W
    @test @inferred mul!(copy(W), M, V, α) ≈ α * AV
    @test @inferred mul!(copy(W), M, V) ≈ AV
    @test typeof(M * V) <: LinearMap
end

@testset "transpose/adjoint" begin
    @test @inferred M' * w == A' * w
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W
    @test @inferred transpose(M) * w == transpose(A) * w
    @test @inferred transpose(LinearMap(transpose(M))) * v == A * v
    @test @inferred LinearMap(M')' * v == A * v
    @test @inferred transpose(transpose(M)) == M
    @test (M')' == M
    Mherm = LinearMap(A'A)
    @test @inferred ishermitian(Mherm)
    @test @inferred !issymmetric(Mherm)
    @test @inferred !issymmetric(transpose(Mherm))
    @test @inferred ishermitian(transpose(Mherm))
    @test @inferred ishermitian(Mherm')
    @test @inferred isposdef(Mherm)
    @test @inferred isposdef(transpose(Mherm))
    @test @inferred isposdef(adjoint(Mherm))
    @test @inferred !(transpose(M) == adjoint(M))
    @test @inferred !(adjoint(M) == transpose(M))
    @test @inferred transpose(M') * v ≈ transpose(A') * v
    @test @inferred transpose(LinearMap(M')) * v ≈ transpose(A') * v
    @test @inferred LinearMap(transpose(M))' * v ≈ transpose(A)' * v
    @test @inferred transpose(LinearMap(transpose(M))) * v ≈ Av
    @test @inferred adjoint(LinearMap(adjoint(M))) * v ≈ Av

    @test @inferred mul!(copy(w), transpose(LinearMap(M')), v) ≈ transpose(A') * v
    @test @inferred mul!(copy(w), LinearMap(transpose(M))', v) ≈ transpose(A)' * v
    @test @inferred mul!(copy(w), transpose(LinearMap(transpose(M))), v) ≈ Av
    @test @inferred mul!(copy(w), adjoint(LinearMap(adjoint(M))), v) ≈ Av
    @test @inferred mul!(copy(V), transpose(M), W) ≈ transpose(A) * W
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W

    B = LinearMap(Symmetric(rand(10, 10)))
    @test transpose(B) == B
    @test B == transpose(B)

    B = LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
    @test adjoint(B) == B
    @test B == B'
end

@testset "function maps" begin
    N = 100
    function myft(v::AbstractVector)
        # not so fast fourier transform
        N = length(v)
        w = zeros(complex(eltype(v)), N)
        for k = 1:N
            kappa = (2*(k-1)/N)*pi
            for n = 1:N
                w[k] += v[n]*exp(kappa*(n-1)*im)
            end
        end
        return w
    end
    MyFT = LinearMap{ComplexF64}(myft, N) / sqrt(N)
    U = Matrix(MyFT) # will be a unitary matrix
    @test @inferred U'U ≈ Matrix{eltype(U)}(I, N, N)

    CS = LinearMap(cumsum, 2)
    @test size(CS) == (2, 2)
    @test @inferred !issymmetric(CS)
    @test @inferred !ishermitian(CS)
    @test @inferred !isposdef(CS)
    @test @inferred !(LinearMaps.ismutating(CS))
    @test @inferred Matrix(CS) == [1. 0.; 1. 1.]
    @test @inferred Array(CS) == [1. 0.; 1. 1.]
    CS = LinearMap(cumsum, 10; ismutating=false)
    v = rand(10)
    cv = cumsum(v)
    @test CS * v == cv
    @test *(CS, v) == cv
    @test_throws ErrorException CS' * v
    CS = LinearMap(cumsum, x -> cumsum(reverse(x)), 10; ismutating=false)
    cv = cumsum(v)
    @test @inferred CS * v == cv
    @test @inferred *(CS, v) == cv
    @test @inferred CS' * v == cumsum(reverse(v))
    @test @inferred mul!(similar(v), transpose(CS), v) == cumsum(reverse(v))

    CS! = LinearMap(cumsum!, 10; ismutating=true)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test_throws ErrorException CS!'v
    @test_throws ErrorException transpose(CS!) * v

    CS! = LinearMap{ComplexF64}(cumsum!, 10; ismutating=true)
    v = rand(ComplexF64, 10)
    cv = cumsum(v)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test_throws ErrorException CS!'v
    @test_throws ErrorException adjoint(CS!) * v
    CS! = LinearMap{ComplexF64}(cumsum!, x -> cumsum!(reverse!(x)), 10; ismutating=true)
    @test @inferred LinearMaps.ismutating(CS!)
    @test @inferred CS! * v == cv
    @test @inferred *(CS!, v) == cv
    @test @inferred mul!(similar(v), CS!, v) == cv
    @test @inferred CS' * v == cumsum(reverse(v))
    @test @inferred mul!(similar(v), transpose(CS), v) == cumsum(reverse(v))
    @test @inferred mul!(similar(v), adjoint(CS), v) == cumsum(reverse(v))

    # Test fallback methods:
    L = LinearMap(x -> x, x -> x, 10)
    v = randn(10)
    @test @inferred (2 * L)' * v ≈ 2 * v
    @test @inferred transpose(2 * L) * v ≈ 2 * v
    L = LinearMap{ComplexF64}(x -> x, x -> x, 10)
    v = rand(ComplexF64, 10)
    @test @inferred (2 * L)' * v ≈ 2 * v
    @test @inferred transpose(2 * L) * v ≈ 2 * v
end

CS! = LinearMap(cumsum!, 10; ismutating=true)
v = rand(10)
u = similar(v)
mul!(u, CS!, v)
@test ((@allocated mul!(u, CS!, v)) == 0)
n = 10
L = sum(fill(CS!, n))
@test mul!(u, L, v) ≈ n * cumsum(v)
@test ((@numalloc mul!(u, L, v)) <= 1)

A = 2 * rand(ComplexF64, (10, 10)) .- 1
B = rand(size(A)...)
M = LinearMap(A)
N = LinearMap(B)
LC = M + N
v = rand(ComplexF64, 10)
w = similar(v)
mul!(w, M, v)
@test ((@allocated mul!(w, M, v)) == 0)
@testset "linear combinations" begin
    # @test_throws ErrorException LinearMaps.LinearCombination{ComplexF64}((M, N), (1, 2, 3))
    @test @inferred size(3M + 2.0N) == size(A)
    # addition
    @test @inferred convert(Matrix, LC) == A + B
    @test @inferred convert(Matrix, LC + LC) ≈ 2A + 2B
    @test @inferred convert(Matrix, M + LC) ≈ 2A + B
    @test @inferred convert(Matrix, M + M) ≈ 2A
    # subtraction
    @test @inferred Matrix(LC - LC) == zeros(size(LC))
    @test @inferred Matrix(LC - M) == B
    @test @inferred Matrix(N - LC) == -A
    @test @inferred Matrix(M - M) == zeros(size(M))
    # scalar multiplication
    @test @inferred Matrix(-M) == -A
    @test @inferred Matrix(-LC) == -A - B
    @test @inferred Matrix(3 * M) == 3 * A
    @test @inferred Matrix(M * 3) == 3 * A
    @test Matrix(3.0 * LC) ≈ Matrix(LC * 3) == 3A + 3B
    @test @inferred Matrix(3 \ M) ≈ A/3
    @test @inferred Matrix(M / 3) ≈ A/3
    @test @inferred Matrix(3 \ LC) ≈ (A + B) / 3
    @test @inferred Matrix(LC / 3) ≈ (A + B) / 3
    @test @inferred Array(3 * M' - CS!) == 3 * A' - Array(CS!)
    @test (3M - 1im*CS!)' == (3M + ((-1im)*CS!))' == M'*3 + CS!'*1im
    @test @inferred Array(M + A) == 2 * A
    @test @inferred Array(A + M) == 2 * A
    @test @inferred Array(M - A) == 0 * A
    @test Array(A - M) == 0 * A
end

# new type
import Base: *, size
import LinearAlgebra: mul!

struct SimpleFunctionMap <: LinearMap{Float64}
    f::Function
    N::Int
end
struct SimpleComplexFunctionMap <: LinearMap{Complex{Float64}}
    f::Function
    N::Int
end
Base.size(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}) = (A.N, A.N)
Base.:(*)(A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, v::Vector) = A.f(v)
mul!(y::Vector, A::Union{SimpleFunctionMap,SimpleComplexFunctionMap}, x::Vector) = copyto!(y, *(A, x))

@testset "composition" begin
    F = LinearMap(cumsum, 10; ismutating=false)
    A = 2 * rand(ComplexF64, (10, 10)) .- 1
    B = rand(size(A)...)
    M = 1 * LinearMap(A)
    N = LinearMap(B)
    @test @inferred (F * F) * v == @inferred F * (F * v)
    @test @inferred (F * A) * v == @inferred F * (A * v)
    @test @inferred (A * F) * v == @inferred A * (F * v)
    @test @inferred A * (F * F) * v == @inferred A * (F * (F * v))
    @test @inferred (F * F) * (F * F) * v == @inferred F * (F * (F * (F * v)))
    @test @inferred Matrix(M * transpose(M)) ≈ A * transpose(A)
    @test @inferred !isposdef(M * transpose(M))
    @test @inferred isposdef(M * M')
    @test @inferred issymmetric(N * N')
    @test @inferred ishermitian(N * N')
    @test @inferred !issymmetric(M' * M)
    @test @inferred ishermitian(M' * M)
    @test @inferred isposdef(transpose(F) * F)
    @test @inferred isposdef(transpose(F) * F * 3)
    @test @inferred isposdef(transpose(F) * 3 * F)
    @test @inferred !isposdef(-5*transpose(F) * F)
    @test @inferred isposdef((M * F)' * M * 4 * F)
    @test @inferred transpose(M * F) == @inferred transpose(F) * transpose(M)
    @test @inferred (4*((-3*M)*2)) == @inferred -12M*2
    @test @inferred (4*((3*(-M))*2)*(-5)) == @inferred -12M*(-10)
    L = 3 * F + 1im * A + F * M' * F
    LF = 3 * Matrix(F) + 1im * A + Matrix(F) * Matrix(M)' * Matrix(F)
    @test Array(L) ≈ LF
    R1 = rand(ComplexF64, 10, 10); L1 = LinearMap(R1)
    R2 = rand(ComplexF64, 10, 10); L2 = LinearMap(R2)
    R3 = rand(ComplexF64, 10, 10); L3 = LinearMap(R3)
    CompositeR = prod(R -> LinearMap(R), [R1, R2, R3])
    @test transpose(CompositeR) == transpose(L3) * transpose(L2) * transpose(L1)
    @test adjoint(CompositeR) == L3' * L2' * L1'
    @test adjoint(adjoint((CompositeR))) == CompositeR
    @test transpose(transpose((CompositeR))) == CompositeR
    Lt = transpose(LinearMap(CompositeR))
    @test Lt * v ≈ transpose(R3) * transpose(R2) * transpose(R1) * v
    Lc = adjoint(LinearMap(CompositeR))
    @test Lc * v ≈ R3' * R2' * R1' * v

    # test inplace operations
    w = similar(v)
    mul!(w, L, v)
    @test w ≈ LF * v

    # test new type
    F = SimpleFunctionMap(cumsum, 10)
    FC = SimpleComplexFunctionMap(cumsum, 10)
    @test @inferred ndims(F) == 2
    @test @inferred size(F, 1) == 10
    @test @inferred length(F) == 100
    @test @inferred !issymmetric(F)
    @test @inferred !ishermitian(F)
    @test @inferred !ishermitian(FC)
    @test @inferred !isposdef(F)
    w = similar(v)
    mul!(w, F, v)
    @test w == F * v
    @test_throws MethodError F' * v
    @test_throws MethodError transpose(F) * v
    @test_throws MethodError mul!(w, adjoint(F), v)
    @test_throws MethodError mul!(w, transpose(F), v)

    # test composition of several maps with shared data #31
    global sizes = ( (5, 2), (3, 3), (3, 2), (2, 2), (9, 2), (7, 1) )
    N = length(sizes) - 1
    global Lf = []
    global Lt = []
    global Lc = []

    # build list of operators [LN, ..., L2, L1] for each mode
    for (fi, i) in [ (Symbol("f$i"), i) for i in 1:N]
        @eval begin
            function ($fi)(source)
                dest = ones(prod(sizes[$i + 1]))
                tmp = reshape(source, sizes[$i])
                return conj.($i * dest)
            end
            insert!(Lf, 1, LinearMap($fi, prod(sizes[$i + 1]), prod(sizes[$i])))
            insert!(Lt, 1, LinearMap(x -> x, $fi, prod(sizes[$i]), prod(sizes[$i + 1])))
            insert!(Lc, 1, LinearMap{ComplexF64}(x -> x, $fi, prod(sizes[$i]), prod(sizes[$i + 1])))
        end
    end
    @test size(prod(Lf[1:N])) == (prod(sizes[end]), prod(sizes[1]))
    @test isreal(prod(Lf[1:N]))
    # multiply as composition and as recursion
    v1 = ones(prod(sizes[1]))
    u1 = ones(prod(sizes[1]))
    w1 = im.*ones(ComplexF64, prod(sizes[1]))
    for i = N:-1:1
        v2 = prod(Lf[i:N]) * ones(prod(sizes[1]))
        u2 = transpose(LinearMap(prod(Lt[N:-1:i]))) * ones(prod(sizes[1]))
        w2 = adjoint(LinearMap(prod(Lc[N:-1:i]))) * ones(prod(sizes[1]))

        v1 = Lf[i] * v1
        u1 = transpose(Lt[i]) * u1
        w1 = adjoint(Lc[i]) * w1

        @test v1 == v2
        @test u1 == u2
        @test w1 == w2
    end
end

A = rand(10, 20)
B = rand(ComplexF64, 10, 20)
SA = A'A + I
SB = B'B + I
L = LinearMap{Float64}(A)
MA = LinearMap(SA)
MB = LinearMap(SB)
@testset "wrapped maps" begin
    @test size(L) == size(A)
    @test @inferred !issymmetric(L)
    @test @inferred issymmetric(MA)
    @test @inferred !issymmetric(MB)
    @test @inferred isposdef(MA)
    @test @inferred isposdef(MB)
end

A = 2 * rand(ComplexF64, (10, 10)) .- 1
B = rand(size(A)...)
M = @inferred 1 * LinearMap(A)
N = LinearMap(B)
LC = M + N
v = rand(ComplexF64, 10)
w = similar(v)
@testset "identity/scaling map" begin
    Id = LinearMaps.UniformScalingMap(1, 10)
    @test_throws ErrorException LinearMaps.UniformScalingMap(1, 10, 20)
    @test_throws ErrorException LinearMaps.UniformScalingMap(1, (10, 20))
    @test size(Id) == (10, 10)
    @test @inferred isreal(Id)
    @test @inferred issymmetric(Id)
    @test @inferred ishermitian(Id)
    @test @inferred isposdef(Id)
    @test Id * v == v
    @test (2 * M' + 3 * I) * v == 2 * A'v + 3v
    @test (3 * I + 2 * M') * v == 2 * A'v + 3v
    @test (2 * M' - 3 * I) * v == 2 * A'v - 3v
    @test (3 * I - 2 * M') * v == -2 * A'v + 3v
    @test transpose(LinearMap(2 * M' + 3 * I)) * v ≈ transpose(2 * A' + 3 * I) * v
    @test LinearMap(2 * M' + 3 * I)' * v ≈ (2 * A' + 3 * I)' * v
end

@testset "noncommutative number type" begin
    using Quaternions
    x = Quaternion.(rand(10), rand(10), rand(10), rand(10))
    v = rand(10)
    A = Quaternion.(rand(10,10), rand(10,10), rand(10,10), rand(10,10))
    α = UniformScaling(Quaternion.(rand(4)...))
    β = UniformScaling(Quaternion.(rand(4)...))
    L = LinearMap(A)
    @test Array(L) == A
    @test Array(α * L) == α * A
    @test Array(L * α) == A * α
    @test Array(α * L) == α * A
    @test Array(L * α ) == A * α
    @test (α * L')' * x ≈ (α * A')' * x
    @test (α * L')' * v ≈ (α * A')' * v
    @test Array(@inferred adjoint(α * L * β)) ≈ conj(β) * A' * conj(α)
    @test Array(@inferred transpose(α * L * β)) ≈ β * transpose(A) * α
end


text:
include("gaussblur.jl")

if !isdir("data")
  print("Downloading data...")
  run(`./get-data.sh`)
  run(`mkdir output`)
  println("done.")
end

using Images
const noise_mean = 0.0021172377176794793
const sigmasq = let
  lambda = 723.0
  NA = 1.4
  FWHM = lambda/(2*NA)
  sigma = FWHM/(2*log(2.0))
  (sigma/(64*100.0))^2
end

const gb_sim = GaussBlur2D(sigmasq,64,500)

function readImages(imageDir,nImages,n_pixels_x)
  imA = Array(Float64,n_pixels_x,n_pixels_x,nImages)
  for q = 1:nImages
    imageId = @sprintf("%05D",q)
    img = imread("$imageDir/$imageId.tif")
    imA[:,:,q] = float(data(img))
  end
  imA
end

function runFW(sim,imageArray; n_cd_iters :: Int64 = 200)
  nImages = size(imageArray,3)
  results = Array(Array{Float64},nImages)
  for imageIdx = 1:nImages
    img = imageArray[:,:,imageIdx]
    target = vec(img).-noise_mean
    function callback(old_thetas,thetas, weights,output,old_obj_val)
      #evalute current OV
      new_obj_val,t = loss(LSLoss(), output - target)
      if old_obj_val - new_obj_val < 7E-5
        return true
      end
      return false
    end
    (points,weights) = ADCG(gb_sim, LSLoss(), target, 200000.0; max_iters=20, callback=callback, max_cd_iters=n_cd_iters)
    results[imageIdx] = points
    print(".")
  end
  results
end

function writeCSV(filename,results)
  header = "Ground-truth,frame,xnano,ynano,znano,intensity\n"
  nImages = size(results,1)
  csvfile = open(filename,"w")
  write(csvfile,header)
  for frame = 1:nImages
    localizations = results[frame]
    nSources = size(localizations,2)
    for i = 1:nSources
      xPos = localizations[1,i]*6400.0
      yPos = localizations[2,i]*6400.0
      write(csvfile,"$i, $frame, $yPos,$xPos, 0.0,0.0\n")
    end
  end
  close(csvfile)
end

nImages = 100;

imageArray = readImages("data/sequence",nImages,64)

print("Running FWA-MI")
results = runFW(gb_sim,imageArray)
println(" done.")
writeCSV("output/fwami.csv",results)
##ablations...
print("Running FWA-M")
results = runFW(gb_sim,imageArray; n_cd_iters=0)
println(" done.")
writeCSV("output/fwam.csv",results)
print("Running GF")
results = runFW(gb_sim,imageArray; n_cd_iters=1)
println(" done.")
writeCSV("output/gf.csv",results)


text:
"Tableau for explicit Euler method"
function getTableauExplicitEuler()
    a = zeros(Float64, 1, 1)
    b = [1.0]
    c = [0.0]
    o = 1

    TableauERK(:explicit_euler, o, a, b, c)
end

"Tableau for explicit midpoint method"
function getTableauExplicitMidpoint()
    a = [[0.0 0.0]
         [0.5 0.0]]
    b = [0.0, 1.0]
    c = [0.0, 0.5]
    o = 2

    TableauERK(:explicit_midpoint, o, a, b, c)
end

"Tableau for Runge's method"
function getTableauRunge()
    a = [[0.0 0.0]
         [1.0 0.0]]
    b = [0.5, 0.5]
    c = [0.0, 1.0]
    o = 2

    TableauERK(:runge, o, a, b, c)
end

"Tableau for Heun's method"
function getTableauHeun()
    a = [[0.0 0.0]
         [1.0 0.0]]
    b = [0.5, 0.5]
    c = [0.0, 1.0]
    o = 2

    TableauERK(:heun, o, a, b, c)
end

"Tableau for Kutta's method of order three"
function getTableauKutta()
    a = [[ 0.0 0.0 0.0]
         [ 0.5 0.0 0.0]
         [-1.0 2.0 0.0]]
    b = [1/6, 4/6, 1/6]
    c = [0.0, 0.5, 1.0]
    o = 3

    TableauERK(:kutta, o, a, b, c)
end

"Tableau for explicit Runge-Kutta method of order four (1/6 rule)"
function getTableauERK4()
    a = [[0.0 0.0 0.0 0.0]
         [0.5 0.0 0.0 0.0]
         [0.0 0.5 0.0 0.0]
         [0.0 0.0 1.0 0.0]]
    b = [1/6, 1/3, 1/3, 1/6]
    c = [0.0, 0.5, 0.5, 1.0]
    o = 4

    TableauERK(:erk4, o, a, b, c)
end

"Tableau for explicit Runge-Kutta method of order four (3/8 rule)"
function getTableauERK438()
    a = [[ 0.0  0.0  0.0  0.0]
         [ 1/3  0.0  0.0  0.0]
         [-1/3  1.0  0.0  0.0]
         [ 1.0 -1.0  1.0  0.0]]
    b = [1/8, 3/8, 3/8, 1/8]
    c = [0.0, 1/3, 2/3, 1.0]
    o = 4

    TableauERK(:erk438, o, a, b, c)
end

"Tableau for explicit Verner's method of order six"
function getTableauVerner()
    a = [[ 0.0          0.0     0.0           0.0       0.0          0.0  0.0         0.0]
         [ 1/6          0.0     0.0           0.0       0.0          0.0  0.0         0.0]
         [ 4/75         16/75   0.0           0.0       0.0          0.0  0.0         0.0]
         [ 5/6         -8/3     5/2           0.0       0.0          0.0  0.0         0.0]
         [-165/64       55/6   -425/64        85/96     0.0          0.0  0.0         0.0]
         [ 12/5        -8.0     4015/612     -11/36     88/255       0.0  0.0         0.0]
         [-8263/15000   124/75 -643/680      -81/250    2484/10625   0.0  0.0         0.0]
         [ 3501/1720   -300/43  297275/52632 -319/2322  24068/84065  0.0  3850/26703  0.0]]
    b = [3/40, 0.0, 875/2244, 23/72, 264/1955, 0.0, 125/11592, 43/616]
    c = [0.0, 1/6, 4/15, 2/3, 5/6, 1.0, 1/15, 1.0]
    o = 6

    TableauERK(:verner, o, a, b, c)
end