text:
#===============================================================================

    amplification.jl

    parse data for amplification analysis from request body

    Author: Tom Price
    Date:   July 2019

===============================================================================#

import DataStructures.OrderedDict
import StaticArrays: SVector
import Memento: debug

#===============================================================================
    function definitions >>
===============================================================================#

## function called by dispatch()
## parse request body into AmpInput struct and call amp_analysis()
"Generic function called by `dispatch`."
function act(
    ::Type{Val{amplification}},
    req             ::Associative;
    out_format      ::OutputFormat = DEFAULT_AMP_OUTPUT_FORMAT
)
    debug(logger, "at act(::Type{Val{amplification}})")
    #
    ## required fields
    @get_calibration_data_from_req(amplification)
    @parse_raw_data_from_req(amplification)
    #
    ## keyword arguments
    kwargs = AMP_FIELD_DEFS |>
        sift(req_key ∘ field(:key)) |>
        mold() do x
			if isa(req[x.key], Int64) && Int_T != Int64
				x.name => Int_T(req[x.key])
			else
				x.name => req[x.key]
			end
        end
    req_key(CQ_METHOD_KEY) &&
        push!(kwargs,
            :cq_method => try
                CqMethod(req[CQ_METHOD_KEY])
            catch
                return ArgumentError("Unrecognized cq method")
            end) ## try
    req_key(CATEG_WELL_VEC_KEY) &&
        push!(kwargs,
            :categ_well_vec =>
                map(req[CATEG_WELL_VEC_KEY]) do x
                    element = str2sym.(x)
                    (length(element[2]) == 0) ?
                        element :
                        Colon()
                end) ## do x
    req_key(CTRL_WELL_DICT_KEY) &&
        push!(kwargs, :ctrl_well_dict => str2sym.(req[CTRL_WELL_DICT_KEY]))
    #
    ## baseline model parameters
    kw_bl =
        begin
            baseline_method =
                req_key(BASELINE_METHOD_KEY) &&
               req[BASELINE_METHOD_KEY]
            if      (baseline_method == SIGMOID_KEY)
                        Dict{Symbol,Any}(
                            :bl_method          =>  l4_enl,
                            :bl_fallback_func   =>  mean)
            elseif  (baseline_method == LINEAR_KEY)
                        Dict{Symbol,Any}(
                            :bl_method          =>  lin_1ft,
                            :bl_fallback_func   =>  mean)
            elseif  (baseline_method == MEDIAN_KEY)
                        Dict{Symbol,Any}(
                            :bl_method          =>  take_the_median)
            else
                Dict{Symbol,Any}()
            end
        end
    #
    ## create container for data and parameters
    ## to pass to amp_analysis()
    interface = AmpInput(
        parsed_raw_data...,
        calibration_data;
        out_format = out_format,
        amp_output = AmpOutputOption(out_format),
        kw_bl...,
        kwargs...,)
    result =
        ## issues:
        ## 1.
        ## the new code currently assumes only 1 step/ramp
        ## because as the request body is currently structured
        ## we cannot subset the fluorescence data by step_id/ramp_id
        ## 2.
        ## need to verify that the fluorescence data complies
        ## with the constraints imposed by max_cycle and well_constraint
        #
        # const sr_key =
        #     if      req_key(STEP_ID_KEY) STEP_ID_KEY
        #     elseif  req_key(RAMP_ID_KEY) RAMP_ID_KEY
        #     else throw(ArgumentError("no step/ramp information found"))
        #     end
        # const asrp_vec = [AmpStepRampProperties(:ramp, req[sr_key], DEFAULT_cyc_nums)]
        # const sr_dict =
        #     OrderedDict(
        #         map([ asrp_vec[1] ]) do asrp
        #             join([asrp.step_or_ramp, asrp.id], "_") =>
        #                 amp_analysis(
        #                     # remove MySql dependency
        #                     # db_conn, exp_id, asrp, calib_info,
        #                     # fluo_well_nums, well_nums,
        #                     amp,
        #                     asrp,
        #                     out_format == json_output ? pre_json_output : out_format))
        #         end) ## do asrp
        # ## output
        # if (out_sr_dict)
        #     final_out = sr_dict
        # else
        #     const first_sr_out = first(values(sr_dict))
        #     final_out =
        #         OrderedDict(
        #             map(fieldnames(first_sr_out)) do key
        #                 key => getfield(first_sr_out, key)
        #             end)
        # end
        first_sr_out = amp_analysis(interface)
        OrderedDict([
            map(fieldnames(first_sr_out)) do key
                key => getfield(first_sr_out, key)
            end...,
            :valid => true])
    return result |> out(out_format)
end ## act(::Type{Val{amplification}})

## called by parse_req >>

@inline str2sym(x) = isa(x, String) ? Symbol(x) : x

#==============================================================================#

"Extract dimensions of raw amplification data, then format the raw data into a
3D array as required by `calibrate`."
function parse_raw_data(::Type{Val{amplification}}, raw_dict ::Associative)
    (cycles, wells, channels) =
        map([CYCLE_NUM_KEY, WELL_NUM_KEY, CHANNEL_KEY]) do key
            raw_dict[key] |> unique ## in order of appearance
        end
    (num_cycles, num_wells, num_channels) =
        map(length, (cycles, wells, channels))
    #
    ## check that data are sorted and conformable to 3D array
    try
        assert(raw_dict[CYCLE_NUM_KEY] ==
            repeat(
                cycles,
                outer = num_wells * num_channels))
        assert(raw_dict[WELL_NUM_KEY ] ==
            repeat(
                wells,
                inner = num_cycles,
                outer = num_channels))
        assert(raw_dict[CHANNEL_KEY  ] ==
            repeat(
                channels,
                inner = num_cycles * num_wells))
    catch err
		warn(logger, "raw dict exception: " * sprint(showerror, err))
        throw(ArgumentError("The format of the fluorescence data does not " *
            "lend itself to transformation into a 3-dimensional array. " *
            "Please make sure that the data are sorted by " *
            "channel, well, and cycle number."))
    end ## try
    #
    ## reshape to 3D array
    F = eltype(first(raw_dict[FLUORESCENCE_VALUE_KEY]))
	if (F == Int64 && Int_T != Int64)
		F = Int_T
	end
    raw_data = ## formerly `fr_ary3`
        reshape(
            raw_dict[FLUORESCENCE_VALUE_KEY],
            num_cycles, num_wells, num_channels)
    #
    #
    ## kludge to index well numbers starting at 0
    return (
        RawData{F}(raw_data[cycles, wells, channels]),
        Int_T(num_cycles),
        Int_T(num_wells),
        Int_T(num_channels),
        cycles[cycles] |> SVector{num_cycles,Int_T},
        wells |> mold(Symbol ∘ Int_T) |> SVector{num_wells,Symbol},
        channels[channels] |> SVector{num_channels,Int_T})
end ## parse_raw_data(Type{Val{amplification}})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
module semiclassical

import Base: ==
import ..bases, ..operators, ..operators_dense
import ..timeevolution: integrate, recast!

using ..bases, ..states, ..operators, ..operators_dense, ..timeevolution

const QuantumState = Union{Ket, DenseOperator}
const DecayRates = Union{Void, Vector{Float64}, Matrix{Float64}}

"""
Semi-classical state.

It consists of a quantum part, which is either a `Ket` or a `DenseOperator` and
a classical part that is specified as a complex vector of arbitrary length.
"""
mutable struct State{T<:QuantumState}
    quantum::T
    classical::Vector{Complex128}
end

Base.length(state::State) = length(state.quantum) + length(state.classical)
Base.copy(state::State) = State(copy(state.quantum), copy(state.classical))

function ==(a::State{T}, b::State{T}) where T<:QuantumState
    samebases(a.quantum, b.quantum) &&
    length(a.classical)==length(b.classical) &&
    (a.classical==b.classical) &&
    (a.quantum==b.quantum)
end

operators.expect(op, state::State) = expect(op, state.quantum)
operators.variance(op, state::State) = variance(op, state.quantum)
operators.ptrace(state::State, indices::Vector{Int}) = State{DenseOperator}(ptrace(state.quantum, indices), state.classical)

operators_dense.dm(x::State{Ket}) = State{DenseOperator}(dm(x.quantum), x.classical)

"""
    semiclassical.schroedinger_dynamic(tspan, state0, fquantum, fclassical[; fout, ...])

Integrate time-dependent Schrödinger equation coupled to a classical system.

# Arguments
* `tspan`: Vector specifying the points of time for which the output should
        be displayed.
* `psi0`: Initial semi-classical state [`semiclassical.State`](@ref).
* `fquantum`: Function `f(t, psi, u) -> H` returning the time and or state
        dependent Hamiltonian.
* `fclassical`: Function `f(t, psi, u, du)` calculating the possibly time and
        state dependent derivative of the classical equations and storing it
        in the vector `du`.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is neither
        normalized nor permanent!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function schroedinger_dynamic(tspan, state0::State{Ket}, fquantum::Function, fclassical::Function;
                fout::Union{Function,Void}=nothing,
                kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)
    dschroedinger_(t, state::State{Ket}, dstate::State{Ket}) = dschroedinger_dynamic(t, state, fquantum, fclassical, dstate)
    x0 = Vector{Complex128}(length(state0))
    recast!(state0, x0)
    state = copy(state0)
    dstate = copy(state0)
    integrate(tspan_, dschroedinger_, x0, state, dstate, fout; kwargs...)
end

"""
    semiclassical.master_dynamic(tspan, state0, fquantum, fclassical; <keyword arguments>)

Integrate time-dependent master equation coupled to a classical system.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial semi-classical state [`semiclassical.State`](@ref).
* `fquantum`: Function `f(t, rho, u) -> (H, J, Jdagger)` returning the time
        and/or state dependent Hamiltonian and Jump operators.
* `fclassical`: Function `f(t, rho, u, du)` calculating the possibly time and
        state dependent derivative of the classical equations and storing it
        in the complex vector `du`.
* `fout=nothing`: If given, this function `fout(t, state)` is called every time
        an output should be displayed. ATTENTION: The given state is not
        permanent!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_dynamic(tspan, state0::State{DenseOperator}, fquantum, fclassical;
                rates::Union{Vector{Float64}, Matrix{Float64}, Void}=nothing,
                fout::Union{Function,Void}=nothing,
                tmp::DenseOperator=copy(state0.quantum),
                kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)
    function dmaster_(t, state::State{DenseOperator}, dstate::State{DenseOperator})
        dmaster_h_dynamic(t, state, fquantum, fclassical, rates, dstate, tmp)
    end
    x0 = Vector{Complex128}(length(state0))
    recast!(state0, x0)
    state = copy(state0)
    dstate = copy(state0)
    integrate(tspan_, dmaster_, x0, state, dstate, fout; kwargs...)
end

function master_dynamic(tspan, state0::State{Ket}, fquantum, fclassical; kwargs...)
    master_dynamic(tspan, dm(state0), fquantum, fclassical; kwargs...)
end

function recast!(state::State, x::Vector{Complex128})
    N = length(state.quantum)
    copy!(x, 1, state.quantum.data, 1, N)
    copy!(x, N+1, state.classical, 1, length(state.classical))
    x
end

function recast!(x::Vector{Complex128}, state::State)
    N = length(state.quantum)
    copy!(state.quantum.data, 1, x, 1, N)
    copy!(state.classical, 1, x, N+1, length(state.classical))
end

function dschroedinger_dynamic(t::Float64, state::State{Ket}, fquantum::Function,
            fclassical::Function, dstate::State{Ket})
    fquantum_(t, psi) = fquantum(t, state.quantum, state.classical)
    timeevolution.timeevolution_schroedinger.dschroedinger_dynamic(t, state.quantum, fquantum_, dstate.quantum)
    fclassical(t, state.quantum, state.classical, dstate.classical)
end

function dmaster_h_dynamic(t::Float64, state::State{DenseOperator}, fquantum::Function,
            fclassical::Function, rates::DecayRates, dstate::State{DenseOperator}, tmp::DenseOperator)
    fquantum_(t, rho) = fquantum(t, state.quantum, state.classical)
    timeevolution.timeevolution_master.dmaster_h_dynamic(t, state.quantum, fquantum_, rates, dstate.quantum, tmp)
    fclassical(t, state.quantum, state.classical, dstate.classical)
end

end # module

ext: jl


text:
membersuri(cli::Client) = "http://$(cli.host):$(cli.port)/$(cli.version)/members"

"""
    members(cli) -> Dict

Returns a dict (id, member) of members in the etcd cluster.
See the [etcd API](https://github.com/coreos/etcd/blob/master/Documentation/v2/members_api.md#list-members)
for more details.
"""
function members(cli::Client)
    resp = request(Requests.get, membersuri(cli), Dict())
    return map(m -> m["id"] => m, resp["members"]) |> Dict
end


text:
"""

# pcor_test

$(SIGNATURES)

Computes the partial correlation between two variables given a set of other
variables.

### Method
```julia
pcor_test(;
* `u::Vector{Symbol}`                  : Variables used to compute correlation
* `q::Int`                             : Number of variables in conditioning set
* `n::Int`                             : Number of observations
* `S::Matrix`                          : Sample covariance matrix
)
```
where:

  u[1], u[2]: Variables used to compute correlation between, remaining variables
  are the conditioning set

### Returns
```julia
* `res::Float64`                       : Correlation between u[1] and u[2]
```
# Extended help

### Acknowledgements

Original author:                       Giovanni M. Marchetti

Translated to Julia:                   Rob J Goedman

### License

The R package ggm is licensed under License: GPL-2.

The Julia translation is licenced under: MIT.

Part of the api, not exported.
"""
function pcor_test(d::DAG, u::SymbolList, q, n)

  r = pcor(d, u)
  df = n - 2 - q
  tval = r * sqrt(df) / sqrt(1 - r*r)
  pv = 2 * cdf(TDist(df), -abs(tval))
  (pv=tval, df=df, pvalue=pv)

end

export
  pcor_test


text:
export
    Grid,
    build,
    boundary,

    Light,

    Permittivity,
    implant!,

    Permeability,

    Simulator,
    next!,
    simulate!

const C = 299792458

struct Grid{T<:Real}
    nx::Int
    ny::Int
    nt::Int

    Δx::T
    Δy::T
    Δt::T

    max_x::T
    max_y::T
    max_t::T
end

"""
    Grid(max_x, max_y, max_t, nx, ny)

Construct a discretization grid for computational domain.

## Arguments

* `max_x`: Linear sizes of first computational domain in meters.
* `max_y`: Linear sizes of second computational domain in meters.
* `max_t`: Maximum computational time in seconds.
* `nx`: Number of discretization grid for first computational domain.
* `ny`: Number of discretization grid for second computational domain.

## Example

```jldoctest
julia> Grid(3e-6, 10e-6, 1e-12, 60, 200);
```
"""
function Grid(max_x, max_y, max_t, nx, ny)
    Δx = max_x / nx
    Δy = max_y / ny

    Δt = 1 / C / √(1/Δx^2 + 1/Δy^2)
    nt = round(Int, max_t/Δt)

    return Grid(nx, ny, nt, Δx, Δy, Δt, max_x, max_y, max_t)
end

Base.size(grid::Grid) = (grid.nx, grid.ny)
Base.size(grid::Grid, d) = d::Integer <= 2 ? size(grid)[d] : 1
function Base.axes(grid::Grid)
    axes_x = grid.Δx * (Base.OneTo(grid.nx) .- grid.nx/2)
    axes_y = grid.Δy * (Base.OneTo(grid.ny) .- grid.Δy/2)

    return (axes_x, axes_y)
end
Base.axes(grid::Grid, d) = d::Integer <= 2 ? axes(grid)[d] : 1
boundary(grid::Grid) = (grid.max_x, grid.max_y)
boundary(grid::Grid, d) = d::Integer <= 2 ? boundary(grid)[d] : 1

function build(grid::Grid)
    return cat(
        repeat(axes(grid, 1), 1, size(grid, 2), 1),
        repeat(axes(grid, 2)', size(grid, 1), 1, 1),
        dims=3
    )
end

struct Light{T<:Real}
    λ::T
    k::T
end

"""
    Light(λ)

Construct a light to propagate.

## Arguments

* `λ`: Wavelength of the light in meters.

## Example

```jldoctest
julia> Light(2.04e-6);
```
"""
function Light(λ)
    return Light(λ, 2π/λ)
end

struct Permittivity{T<:AbstractMatrix}
    ϵ::T
    ϵx::T
    ϵy::T
end

"""
    Permittivity(ϵ_const::Real, grid::Grid)

Construct a uniform permittivity for the medium.

## Arguments

* `ϵ_const`: Permittivity of the medium in F/m.
* `grid`: Discretization grid.

## Example

```jldoctest
julia> grid = Grid(3e-6, 10e-6, 1e-12, 60, 200);

julia> Permittivity(9., grid);
```
"""
function Permittivity(ϵ_const::Real, grid::Grid)
    ϵ = ϵ_const * ones(size(grid))
    ϵx = C * grid.Δt/grid.Δx ./ ϵ
    ϵy = C * grid.Δt/grid.Δy ./ ϵ

    return Permittivity(ϵ, ϵx, ϵy)
end

"""
    implant!(
        permittivity::Permittivity, ϵ_const::Real,
        xs::AbstractVector, ys::AbstractVector, rs::AbstractVector,
        grid::Grid
    )

Implant some bubble defect into the medium.

## Arguments

* `permittivity`: Permittivity object.
* `ϵ_const`: Permittivity of the defect in F/m.
* `xs`: Position of defect of first computational domain in meters.
* `ys`: Position of defect of second computational domain in meters.
* `rs`: Radius of defect in meters.
* `grid`: Discretization grid.

## Example

```jldoctest
julia> grid = Grid(3e-6, 10e-6, 1e-12, 60, 200);

julia> permittivity = Permittivity(9., grid);

julia> ϵ_defect = 1.;

julia> xs_defect = [0, 1e-6, -1e-6];

julia> ys_defect = [1e-6, 2e-6, 3e-6];

julia> rs_defect = [0.5e-6, 0.1e-6, 0.2e-6];

julia> implant!(
           permittivity, ϵ_defect,
           xs_defect, ys_defect, rs_defect,
           grid
       );
```
"""
function implant!(
    permittivity::Permittivity, ϵ_const::Real,
    xs::AbstractVector, ys::AbstractVector, rs::AbstractVector,
    grid::Grid
)
    length(xs) == length(ys) == length(rs) || throw(DimensionMismatch("xs, ys, rs must have same length"))

    in_circle(i, j) = true in [
        √(
            (axes(grid, 1)[i] - xs[c])^2 +
            (axes(grid, 2)[j] - ys[c])^2
        ) < rs[c]
        for c in 1:length(xs)
    ]

    for i in 1:size(grid, 1), j in 1:size(grid, 2)
        in_circle(i, j) && (permittivity.ϵ[i, j] = ϵ_const)
    end

    permittivity.ϵx .= C * grid.Δt/grid.Δx ./ permittivity.ϵ
    permittivity.ϵy .= C * grid.Δt/grid.Δy ./ permittivity.ϵ

    return permittivity
end

struct Permeability{T<:AbstractMatrix}
    μ::T
    μx::T
    μy::T
end

"""
    Permeability(μ_const::Real, grid::Grid)

Construct a uniform permeability for the medium.

## Arguments

* `μ_const`: Permeability of the medium in N/A².
* `grid`: Discretization grid.

## Example

```jldoctest
julia> grid = Grid(3e-6, 10e-6, 1e-12, 60, 200);

julia> Permeability(1., grid);
```
"""
function Permeability(μ_const::Real, grid::Grid)
    μ = μ_const * ones(size(grid))
    μx = C * grid.Δt/grid.Δx ./ μ
    μy = C * grid.Δt/grid.Δy ./ μ

    return Permeability(μ, μx, μy)
end

mutable struct Simulator{T<:AbstractArray}
    grid::Grid
    light::Light
    permittivity::Permittivity
    permeability::Permeability

    ez::T
    hx::T
    hy::T

    t::Int
end

function get_default_e_field(light::Light, grid::Grid; t=1)
    Δt = grid.Δt
    k = light.k

    return 0.1 * exp.(
        -axes(grid, 1)[2:end].^2 ./
        (boundary(grid, 1)/4)^2
    ) * sin(k * C*Δt*t)
end

"""
    Simulator(grid::Grid, light::Light, permittivity::Permittivity, permeability::Permeability)

## Arguments

* `grid`: Discretization grid for computational domain.
* `light`: Light to propagate.
* `permittivity`: Permittivity for the medium.
* `permeability`: Permeability for the medium.

## Example

```jldoctest
julia> grid = Grid(3e-6, 10e-6, 1e-12, 60, 200);

julia> light = Light(2.04e-6);

julia> permittivity = Permittivity(9., grid);

julia> permeability = Permeability(1., grid);

julia> Simulator(grid, light, permittivity, permeability);
"""
function Simulator(grid::Grid, light::Light, permittivity::Permittivity, permeability::Permeability)
    ez = zeros(Float64, size(grid))
    ez[2:end, 1] .= get_default_e_field(light, grid)

    return Simulator(
        grid,
        light,
        permittivity,
        permeability,

        ez,
        zeros(Float64, size(grid)),
        zeros(Float64, size(grid)),

        1
    )
end

function Simulator(;
    max_x=3e-6, max_y=10e-6, max_t=0.1e-12,
    nx=300, ny=1000,
    λ=2.04e-6,
    ϵ = 9., μ = 1.
)
    grid = Grid(max_x, max_y, max_t, nx, ny)
    light = Light(λ)
    permittivity = Permittivity(ϵ, grid)
    permeability = Permeability(μ, grid)

    return Simulator(grid, light, permittivity, permeability)
end

function next!(s::Simulator)
    s.t += 1

    ϵx, ϵy = s.permittivity.ϵx, s.permittivity.ϵy
    μx, μy = s.permeability.μx, s.permeability.μy

    s.ez[2:end, 1] .+= get_default_e_field(s.light, s.grid, t=s.t)

    s.hx[2:end-1, 2:end-1] .+= -μx[2:end-1, 2:end-1].*(s.ez[2:end-1, 2:end-1] - s.ez[2:end-1, 1:end-2])
    s.hy[2:end-1, 2:end-1] .+= +μy[2:end-1, 2:end-1].*(s.ez[2:end-1, 2:end-1] - s.ez[1:end-2, 2:end-1])

    s.ez[2:end-1, 2:end-1] .+=
        ϵx[2:end-1, 2:end-1].*(s.hy[3:end, 2:end-1] - s.hy[2:end-1, 2:end-1]) -
        ϵy[2:end-1, 2:end-1].*(s.hx[2:end-1, 3:end] - s.hx[2:end-1, 2:end-1])

    return s
end

"""
    simulate!(s::Simulator)

Run simulation from current `t` to `max_t`

## Arguments

* `s`: Simulator.

## Example

```julia
julia> grid = Grid(3e-6, 10e-6, 1e-12, 60, 200);

julia> light = Light(2.04e-6);

julia> permittivity = Permittivity(9., grid);

julia> permeability = Permeability(1., grid);

julia> s = Simulator(grid, light, permittivity, permeability);

julia> simulate!(s);
```
"""
function simulate!(s::Simulator)
    nt = s.grid.nt

    p = Progress(nt)
    for _ in (s.t+1):nt
        next!(s)
        ProgressMeter.next!(p)
    end

    return s
end


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------

@testset ExtendedTestSet "density matrix" begin

    pauli = [
        Matrix{Float64}(I, 2, 2),
        Qaintessent.matrix(X),
        Qaintessent.matrix(Y),
        Qaintessent.matrix(Z),
    ]
    @testset "density matrix pauli group matrix" begin
        @test Qaintessent.matrix(pauli_group_matrix("XIZY")) ≈
            kron(pauli[3], pauli[4], pauli[1], pauli[2])
    end

    @testset "density matrix constructor" begin
        N = 3

        v = randn(Float64, 4^N)
        d1 = DensityMatrix(v)
        d2 = DensityMatrix(v, N)
        @test d1.v ≈ d2.v
        @test d1.N == d2.N
    end

    @testset "density matrix constructor helper function" begin
        N = 4
        ψ = randn(ComplexF64, 2^N)
        # matrix representation of density matrix |ψ⟩⟨ψ|
        ρmat = reshape(kron(conj(ψ), ψ), 2^N, 2^N)
        # density matrix represented with respect to Pauli basis
        ρ = density_from_statevector(ψ)
        @test Qaintessent.matrix(ρ) ≈ ρmat
        # compute Pauli representation from matrix form
        ρalt = density_from_matrix(ρmat)
        @test ρalt.v ≈ ρ.v
    end

end


text:
using Revise
using VolterraGP
using Flux
D = 3 # number of outputs
C = 2 # order of volterra series
P = 1 # number of paramters in base kernel

train, test = generate_toy_data(20) # genrates toy data from paper
gp = GaussianProcess(threeEQs, D, C, P, train)

t = collect(0.0:0.1:1.)
Gs = collect(0.0:0.1:1.)

g = gradient(x -> sum(VolterraGP.get_phi_cov(0.1, 0.3, 2, 2, 1, 1, x)), gp)
# g[1]

b = gradient(x -> sum(VolterraGP.fill_phi(threeEQs, t, x, [0.1])), Gs)

println(gradient(VolterraGP.full_E, 0.1, 2, gp))
println(gradient(VolterraGP.full_cov, 0.1, 0.2, 2, 2, gp))
gradient(VolterraGP.kernel, 0.1, 0.1, 3, 2, gp)


text:
using Test
@time @testset "Linear deformation 1" begin include("test_linear_defor_1.jl") end
@time @testset "Linear deformation 2" begin include("test_linear_defor_2.jl") end
@time @testset "Linear deformation 3" begin include("test_linear_defor_3.jl") end
@time @testset "Linear deformation 4" begin include("test_linear_defor_4.jl") end
@time @testset "Linear deformation 5" begin include("test_linear_defor_5.jl") end
@time @testset "Linear deformation 6" begin include("test_linear_defor_6.jl") end
@time @testset "Linear deformation 7" begin include("test_linear_defor_7.jl") end
@time @testset "Linear deformation 8" begin include("test_linear_defor_8.jl") end
@time @testset "Linear deformation 9" begin include("test_linear_defor_9.jl") end
@time @testset "Linear deformation 10" begin include("test_linear_defor_10.jl") end
true


text:
include("k-means.jl")

function gaussian_weight{T<:Number}(x_i::Vector{T}, x_j::Vector{T}; sigma=1)
    exp(-norm(x_i - x_j)^2 / (2sigma^2))
end

function gaussian_weight(normval::Float64; sigma=1)
    exp(-normval^2 / (2sigma^2))
end

"""
Do spectral clustering based on Ng's algorithm (Ng, Jordan, Weiss, & others, 2002).

    @inputs the input data matrix, each column is a vector in feature space.
    @cluster_num the desired number of clusters

    returns an array, whose elements are pairs of array index and the cluster id
"""
function spectral_clustering(inputs::Matrix{Float64}, cluster_num::Int;
    neighbor_num::Int=10, sigma=1)

    W, D = build_matrix_normal_knn(inputs, k=neighbor_num, sigma=sigma)
    L_sym = Symmetric(D^(-1/2) * W * D^(-1/2))

    eigenvals, eigenvecs = eig(L_sym)

    # normalize the eigen vectors
    N = length(eigenvals)
    U = eigenvecs[:, (N - cluster_num + 1):N]
    for row = 1:size(U)[1]
        U[row, :] /= norm(U[row, :])
    end

    centroids, nearest_centroids = kmeans(cluster_num, transpose(U))

    return nearest_centroids
end

using NearestNeighbors
"""
Build a Adjacency Matrix using Normal Function as similarity function
    and only k-nearest-neighborhood have edges connected.

    @inputs input data as a matrix, each column is an input sample
    @k only the nearest k samples have an edge connected with the given point
    @sigma the parameter of the gaussian function

    return the pair of the adjacency matrix and the diagnoal matrix
"""
function build_matrix_normal_knn(inputs::Matrix{Float64}; k::Int=10, sigma=1)
    kdtree = KDTree(inputs)
    N = size(inputs)[2]
    adj = zeros(N, N)

    for i = 1:N
        idxs, dists = knn(kdtree, inputs[:, i], k, true)
        for j in 1:length(idxs)
            weight = idxs[j] == i ? 0 : gaussian_weight(dists[j], sigma=sigma)
            adj[i, idxs[j]] = weight
        end
    end

    adj = (adj + transpose(adj)) / 2.0

    degree = diagm(collect(sum(adj, 1)))

    return adj, degree
end

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libzmq"], :libzmq),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/kdheepak/ZMQBuilder/releases/download/v4.3.1"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.aarch64-linux-gnu-gcc4.tar.gz", "32110d26af10998d084e95d2e664aa1071fe916aff9144f5cd36f668e3d991d8"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.aarch64-linux-gnu-gcc7.tar.gz", "0ade64492d64466a9fb83c7beccf3c494734e74ed2a6421c9b755249b317402f"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.aarch64-linux-gnu-gcc8.tar.gz", "b0fc3eb8fdd9a359490cddc2bdeee89af9c162582f58028dd31ff8ee4c5abd9e"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.aarch64-linux-musl-gcc4.tar.gz", "b97edf71aa975493c791e031b8092a12bdcbc5a668a615eba885365e6c799677"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.aarch64-linux-musl-gcc7.tar.gz", "8abb6a0c7f3d6dea76cc19b5dbb49b20b5dd6f2410c201eecab2a8eb7ee09582"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.aarch64-linux-musl-gcc8.tar.gz", "7da883df0a898eb6e3f9d65b305de7f44341de8a3f2aa86721019c09659e7a4f"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.arm-linux-gnueabihf-gcc4.tar.gz", "6a764a07ea2a56fb4783cc46330c07f9caefa060965d7294707aca19e0a25df7"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.arm-linux-gnueabihf-gcc7.tar.gz", "edbcc82e6f9f1ebc4016e8e352c886195c45c97c54c7312bd366aa98bdb72078"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.arm-linux-gnueabihf-gcc8.tar.gz", "22cf03d9e9601c1269e5d9ae1f0b40836f2225a088d346af233ee390999e711a"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.arm-linux-musleabihf-gcc4.tar.gz", "a01549d9364fa4715b13c10b273e8eb3a4e3b4565d5aaed6f787358373c71f90"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.arm-linux-musleabihf-gcc7.tar.gz", "118c614e190c258e2bff4e4203485d82493b3722836fe57d648446557a1ed6a8"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.arm-linux-musleabihf-gcc8.tar.gz", "8526c8f3bf548333c6d31b7b520a653ff0661f0726f6fba8d18f48a68458b164"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.i686-linux-gnu-gcc4.tar.gz", "8bb5bf83793908aec33a2ab89b91b9e9f77822466b7dcb18a12215a0350a16ab"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.i686-linux-gnu-gcc7.tar.gz", "4b93e61c593c994bad3ba6ea686a737ac2d4794ac69ba057363cd61ae4944ffd"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.i686-linux-gnu-gcc8.tar.gz", "59d162a4327882f9f0225cbf994d6bfc70e359ec732a1aec9a9df0c2560f1242"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.i686-linux-musl-gcc4.tar.gz", "87f7069be7eb1965c2b2cb6abec3383cfdc9ca3d2a41ded63cb46890fdb32bd0"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.i686-linux-musl-gcc7.tar.gz", "0e8638dcbe602d2c05846507578479bac208bedea0f9c0507cddbe672e401914"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.i686-linux-musl-gcc8.tar.gz", "7e6a523d0be196f87905ea574b23ef826d1d3db51ec7f37876a4fe9b04e1fc60"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.i686-w64-mingw32-gcc4.tar.gz", "cd1996aaf802d4fc460e1e49477d5ffd323fd608af298305f498062c48a2e444"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.i686-w64-mingw32-gcc7.tar.gz", "f80568f2be31b2c13a2741364b5ea4335269567a34e6b71db6af3db524e4d4be"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.i686-w64-mingw32-gcc8.tar.gz", "9d4ab8fa1aabfadf8664be93bdf8d026e7a3cd05f241e0bd059a0afa620d867f"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.powerpc64le-linux-gnu-gcc4.tar.gz", "633ac5116e73378059656733401c898e42f6304e66a468f95a3433e2e0a8ad7f"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.powerpc64le-linux-gnu-gcc7.tar.gz", "a567ee80e0ea3bb98785d31d8f3a1c8868577c97433b1bbb4101921c21400b43"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.powerpc64le-linux-gnu-gcc8.tar.gz", "5387034d606384e38e093e5897792554e872a013510cecf8755566a2fea5b60a"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-apple-darwin14-gcc4.tar.gz", "3283783620ba8fc3cc67b6e820ed0684f33932a7b580f3f92db9a184a3122b5d"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-apple-darwin14-gcc7.tar.gz", "7d3130be3082d0d2d48f5693134b37ce635dc5a3cc8686e3b09e3678ae29b4d1"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-apple-darwin14-gcc8.tar.gz", "5e080d00520737d77dff8d3b5cce9cc301449bde017fdf1eb1497db7877b3723"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-linux-gnu-gcc4.tar.gz", "a9442cb019947cfc49cec132591fb4e49b9347cfda25919babf48b1657e91f80"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-linux-gnu-gcc7.tar.gz", "98d55b47be4d1dfea853d61ad5d32e6b85d601d7c0c94692462a6a2e459bda67"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-linux-gnu-gcc8.tar.gz", "c111e97855ae5f21719bcda8162ce4c49dea062d8ca7ab79c59f104ecb3c4c6d"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-linux-musl-gcc4.tar.gz", "b73c41fbd13f246b2b661837d61eeedfdf4665f25984f6f3842f4b1075d887ac"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-linux-musl-gcc7.tar.gz", "90984e3275674235add896fb710fa1babbf083af10a2c31cc371303acfc0c548"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-linux-musl-gcc8.tar.gz", "7f34dae4e7dc65b086c5e77de74f81317d6f31bbfdf8dde45b5477b90e2c584a"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-unknown-freebsd11.1-gcc4.tar.gz", "84516762c1dac8005727e552db0952b20a332df8ab2d3a335b3db662c64a0f68"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-unknown-freebsd11.1-gcc7.tar.gz", "6dce8c5dfaeecf4e8af2deda837578e48c4da1ac55f87efac4aac3e192a4c592"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-unknown-freebsd11.1-gcc8.tar.gz", "53a74d406e607275e1494c64277554032946fe2739bc554c35615c19346b1cb3"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-w64-mingw32-gcc4.tar.gz", "9313d63adb3f2e7559c4572e29417007bf26b57b3be1b26399a7747b25f7a5fb"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-w64-mingw32-gcc7.tar.gz", "8b1bfaf639b1e0e9dc283ac85f73b27003cc87ae986d6696b4a300a1257cf3c6"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/ZMQ.v4.3.1.x86_64-w64-mingw32-gcc8.tar.gz", "15be74887936068c0bae5d75d4f146c04d3bb47e7d03c0b82a272d679b8351fc"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
mutable struct PipelineLayout
    mVkDevice::vk.VkDevice
    mDescriptors::Vector{DescriptorSetLayout}
    mDescriptorHandles::Vector{vk.VkDescriptorSetLayout}

    mPushConstants::Vector{vk.VkPushConstantRange}
    mCreateInfo::PipelineLayoutCreateInfo

    mHandleRef::Ref{vk.VkPipelineLayout}

    function PipelineLayout(device::vk.VkDevice,
                       descriptors::Vector{DescriptorSetLayout},
                     pushConstants::Vector{vk.VkPushConstantRange})
        this = new()
        this.mVkDevice = device
        this.mDescriptors = descriptors
        this.mPushConstants = pushConstants
        this.mCreateInfo = PipelineLayoutCreateInfo()
        this.mDescriptorHandles = Vector{vk.VkDescriptorSetLayout}()

        for d in descriptors
            handle = handleRef(d)[]
            push!(this.mDescriptorHandles, handle)
            addSetLayout(this.mCreateInfo, handle)
        end
        for p in pushConstants
            addPushConstantRange(this.mCreateInfo, p)
        end

        this.mHandleRef = Ref{vk.VkPipelineLayout}()
        if vk.vkCreatePipelineLayout(this.mVkDevice, handleRef(this.mCreateInfo), C_NULL, this.mHandleRef) != vk.VK_SUCCESS
            error("Failed to create pipeline layout!")
        end
        return this
    end
end

@class PipelineLayout

function destroy(this::PipelineLayout)
    vk.vkDestroyPipelineLayout(this.mVkDevice, this.mHandleRef[], C_NULL)
end

function handleRef(this::PipelineLayout)::Ref{vk.VkPipelineLayout}
    return this.mHandleRef
end

function getLogicalDevice(this::PipelineLayout)::vk.VkDevice
    return this.mVkDevice
end

function getCreateInfo(this::PipelineLayout)::PipelineLayoutCreateInfo
    return this.mCreateInfo
end


text:
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/cvode/cvode.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __CVodeCreate(lmm::Cint, iter::Cint)
    ccall((:CVodeCreate, libsundials_cvode), CVODEMemPtr, (Cint, Cint), lmm, iter)
end

function CVodeCreate(lmm, iter)
    __CVodeCreate(convert(Cint, lmm), convert(Cint, iter))
end

function CVodeSetErrHandlerFn(cvode_mem, ehfun::CVErrHandlerFn, eh_data)
    ccall((:CVodeSetErrHandlerFn, libsundials_cvode), Cint, (CVODEMemPtr, CVErrHandlerFn, Ptr{Void}), cvode_mem, ehfun, eh_data)
end

function CVodeSetErrFile(cvode_mem, errfp)
    ccall((:CVodeSetErrFile, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{FILE}), cvode_mem, errfp)
end

function CVodeSetUserData(cvode_mem, user_data)
    ccall((:CVodeSetUserData, libsundials_cvode), Cint, (CVODEMemPtr, Any), cvode_mem, user_data)
end

function __CVodeSetMaxOrd(cvode_mem, maxord::Cint)
    ccall((:CVodeSetMaxOrd, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, maxord)
end

function CVodeSetMaxOrd(cvode_mem, maxord)
    __CVodeSetMaxOrd(cvode_mem, convert(Cint, maxord))
end

function __CVodeSetMaxNumSteps(cvode_mem, mxsteps::Clong)
    ccall((:CVodeSetMaxNumSteps, libsundials_cvode), Cint, (CVODEMemPtr, Clong), cvode_mem, mxsteps)
end

function CVodeSetMaxNumSteps(cvode_mem, mxsteps)
    __CVodeSetMaxNumSteps(cvode_mem, convert(Clong, mxsteps))
end

function __CVodeSetMaxHnilWarns(cvode_mem, mxhnil::Cint)
    ccall((:CVodeSetMaxHnilWarns, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, mxhnil)
end

function CVodeSetMaxHnilWarns(cvode_mem, mxhnil)
    __CVodeSetMaxHnilWarns(cvode_mem, convert(Cint, mxhnil))
end

function __CVodeSetStabLimDet(cvode_mem, stldet::Cint)
    ccall((:CVodeSetStabLimDet, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, stldet)
end

function CVodeSetStabLimDet(cvode_mem, stldet)
    __CVodeSetStabLimDet(cvode_mem, convert(Cint, stldet))
end

function CVodeSetInitStep(cvode_mem, hin::realtype)
    ccall((:CVodeSetInitStep, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, hin)
end

function CVodeSetMinStep(cvode_mem, hmin::realtype)
    ccall((:CVodeSetMinStep, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, hmin)
end

function CVodeSetMaxStep(cvode_mem, hmax::realtype)
    ccall((:CVodeSetMaxStep, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, hmax)
end

function CVodeSetStopTime(cvode_mem, tstop::realtype)
    ccall((:CVodeSetStopTime, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, tstop)
end

function __CVodeSetMaxErrTestFails(cvode_mem, maxnef::Cint)
    ccall((:CVodeSetMaxErrTestFails, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, maxnef)
end

function CVodeSetMaxErrTestFails(cvode_mem, maxnef)
    __CVodeSetMaxErrTestFails(cvode_mem, convert(Cint, maxnef))
end

function __CVodeSetMaxNonlinIters(cvode_mem, maxcor::Cint)
    ccall((:CVodeSetMaxNonlinIters, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, maxcor)
end

function CVodeSetMaxNonlinIters(cvode_mem, maxcor)
    __CVodeSetMaxNonlinIters(cvode_mem, convert(Cint, maxcor))
end

function __CVodeSetMaxConvFails(cvode_mem, maxncf::Cint)
    ccall((:CVodeSetMaxConvFails, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, maxncf)
end

function CVodeSetMaxConvFails(cvode_mem, maxncf)
    __CVodeSetMaxConvFails(cvode_mem, convert(Cint, maxncf))
end

function CVodeSetNonlinConvCoef(cvode_mem, nlscoef::realtype)
    ccall((:CVodeSetNonlinConvCoef, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, nlscoef)
end

function __CVodeSetIterType(cvode_mem, iter::Cint)
    ccall((:CVodeSetIterType, libsundials_cvode), Cint, (CVODEMemPtr, Cint), cvode_mem, iter)
end

function CVodeSetIterType(cvode_mem, iter)
    __CVodeSetIterType(cvode_mem, convert(Cint, iter))
end

function CVodeSetRootDirection(cvode_mem, rootdir)
    ccall((:CVodeSetRootDirection, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Cint}), cvode_mem, rootdir)
end

function CVodeSetNoInactiveRootWarn(cvode_mem)
    ccall((:CVodeSetNoInactiveRootWarn, libsundials_cvode), Cint, (CVODEMemPtr,), cvode_mem)
end

function __CVodeInit(cvode_mem, f::CVRhsFn, t0::realtype, y0::N_Vector)
    ccall((:CVodeInit, libsundials_cvode), Cint, (CVODEMemPtr, CVRhsFn, realtype, N_Vector), cvode_mem, f, t0, y0)
end

function CVodeInit(cvode_mem, f, t0, y0)
    __y0 = convert(NVector, y0)
    __CVodeInit(cvode_mem, CVRhsFn_wrapper(f), t0, convert(N_Vector, __y0))
end

function __CVodeReInit(cvode_mem, t0::realtype, y0::N_Vector)
    ccall((:CVodeReInit, libsundials_cvode), Cint, (CVODEMemPtr, realtype, N_Vector), cvode_mem, t0, y0)
end

function CVodeReInit(cvode_mem, t0, y0)
    __y0 = convert(NVector, y0)
    __CVodeReInit(cvode_mem, t0, convert(N_Vector, __y0))
end

function CVodeSStolerances(cvode_mem, reltol::realtype, abstol::realtype)
    ccall((:CVodeSStolerances, libsundials_cvode), Cint, (CVODEMemPtr, realtype, realtype), cvode_mem, reltol, abstol)
end

function __CVodeSVtolerances(cvode_mem, reltol::realtype, abstol::N_Vector)
    ccall((:CVodeSVtolerances, libsundials_cvode), Cint, (CVODEMemPtr, realtype, N_Vector), cvode_mem, reltol, abstol)
end

function CVodeSVtolerances(cvode_mem, reltol, abstol)
    __abstol = convert(NVector, abstol)
    __CVodeSVtolerances(cvode_mem, reltol, convert(N_Vector, __abstol))
end

function CVodeWFtolerances(cvode_mem, efun::CVEwtFn)
    ccall((:CVodeWFtolerances, libsundials_cvode), Cint, (CVODEMemPtr, CVEwtFn), cvode_mem, efun)
end

function __CVodeRootInit(cvode_mem, nrtfn::Cint, g::CVRootFn)
    ccall((:CVodeRootInit, libsundials_cvode), Cint, (CVODEMemPtr, Cint, CVRootFn), cvode_mem, nrtfn, g)
end

function CVodeRootInit(cvode_mem, nrtfn, g)
    __CVodeRootInit(cvode_mem, convert(Cint, nrtfn), CVRootFn_wrapper(g))
end

function __CVode(cvode_mem, tout::realtype, yout::N_Vector, tret, itask::Cint)
    ccall((:CVode, libsundials_cvode), Cint, (CVODEMemPtr, realtype, N_Vector, Ptr{realtype}, Cint), cvode_mem, tout, yout, tret, itask)
end

function CVode(cvode_mem, tout, yout, tret, itask)
    __yout = convert(NVector, yout)
    __CVode(cvode_mem, tout, convert(N_Vector, __yout), tret, convert(Cint, itask))
end

function __CVodeGetDky(cvode_mem, t::realtype, k::Cint, dky::N_Vector)
    ccall((:CVodeGetDky, libsundials_cvode), Cint, (CVODEMemPtr, realtype, Cint, N_Vector), cvode_mem, t, k, dky)
end

function CVodeGetDky(cvode_mem, t, k, dky)
    __dky = convert(NVector, dky)
    __CVodeGetDky(cvode_mem, t, convert(Cint, k), convert(N_Vector, __dky))
end

function CVodeGetWorkSpace(cvode_mem, lenrw, leniw)
    ccall((:CVodeGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrw, leniw)
end

function CVodeGetNumSteps(cvode_mem, nsteps)
    ccall((:CVodeGetNumSteps, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nsteps)
end

function CVodeGetNumRhsEvals(cvode_mem, nfevals)
    ccall((:CVodeGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevals)
end

function CVodeGetNumLinSolvSetups(cvode_mem, nlinsetups)
    ccall((:CVodeGetNumLinSolvSetups, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nlinsetups)
end

function CVodeGetNumErrTestFails(cvode_mem, netfails)
    ccall((:CVodeGetNumErrTestFails, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, netfails)
end

function CVodeGetLastOrder(cvode_mem, qlast)
    ccall((:CVodeGetLastOrder, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Cint}), cvode_mem, qlast)
end

function CVodeGetCurrentOrder(cvode_mem, qcur)
    ccall((:CVodeGetCurrentOrder, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Cint}), cvode_mem, qcur)
end

function CVodeGetNumStabLimOrderReds(cvode_mem, nslred)
    ccall((:CVodeGetNumStabLimOrderReds, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nslred)
end

function CVodeGetActualInitStep(cvode_mem, hinused)
    ccall((:CVodeGetActualInitStep, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, hinused)
end

function CVodeGetLastStep(cvode_mem, hlast)
    ccall((:CVodeGetLastStep, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, hlast)
end

function CVodeGetCurrentStep(cvode_mem, hcur)
    ccall((:CVodeGetCurrentStep, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, hcur)
end

function CVodeGetCurrentTime(cvode_mem, tcur)
    ccall((:CVodeGetCurrentTime, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, tcur)
end

function CVodeGetTolScaleFactor(cvode_mem, tolsfac)
    ccall((:CVodeGetTolScaleFactor, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{realtype}), cvode_mem, tolsfac)
end

function __CVodeGetErrWeights(cvode_mem, eweight::N_Vector)
    ccall((:CVodeGetErrWeights, libsundials_cvode), Cint, (CVODEMemPtr, N_Vector), cvode_mem, eweight)
end

function CVodeGetErrWeights(cvode_mem, eweight)
    __eweight = convert(NVector, eweight)
    __CVodeGetErrWeights(cvode_mem, convert(N_Vector, __eweight))
end

function __CVodeGetEstLocalErrors(cvode_mem, ele::N_Vector)
    ccall((:CVodeGetEstLocalErrors, libsundials_cvode), Cint, (CVODEMemPtr, N_Vector), cvode_mem, ele)
end

function CVodeGetEstLocalErrors(cvode_mem, ele)
    __ele = convert(NVector, ele)
    __CVodeGetEstLocalErrors(cvode_mem, convert(N_Vector, __ele))
end

function CVodeGetNumGEvals(cvode_mem, ngevals)
    ccall((:CVodeGetNumGEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, ngevals)
end

function CVodeGetRootInfo(cvode_mem, rootsfound)
    ccall((:CVodeGetRootInfo, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Cint}), cvode_mem, rootsfound)
end

function CVodeGetIntegratorStats(cvode_mem, nsteps, nfevals, nlinsetups, netfails, qlast, qcur, hinused, hlast, hcur, tcur)
    ccall((:CVodeGetIntegratorStats, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Cint}, Ptr{Cint}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), cvode_mem, nsteps, nfevals, nlinsetups, netfails, qlast, qcur, hinused, hlast, hcur, tcur)
end

function CVodeGetNumNonlinSolvIters(cvode_mem, nniters)
    ccall((:CVodeGetNumNonlinSolvIters, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nniters)
end

function CVodeGetNumNonlinSolvConvFails(cvode_mem, nncfails)
    ccall((:CVodeGetNumNonlinSolvConvFails, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nncfails)
end

function CVodeGetNonlinSolvStats(cvode_mem, nniters, nncfails)
    ccall((:CVodeGetNonlinSolvStats, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, nniters, nncfails)
end

function __CVodeGetReturnFlagName(flag::Clong)
    ccall((:CVodeGetReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVodeGetReturnFlagName(flag)
    __CVodeGetReturnFlagName(convert(Clong, flag))
end

function CVodeFree(cvode_mem)
    ccall((:CVodeFree, libsundials_cvode), Void, (Ref{CVODEMemPtr},), cvode_mem)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/cvode/cvode_bandpre.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function CVBandPrecInit(cvode_mem, N::sunindextype, mu::sunindextype, ml::sunindextype)
    ccall((:CVBandPrecInit, libsundials_cvode), Cint, (CVODEMemPtr, sunindextype, sunindextype, sunindextype), cvode_mem, N, mu, ml)
end

function CVBandPrecGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVBandPrecGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVBandPrecGetNumRhsEvals(cvode_mem, nfevalsBP)
    ccall((:CVBandPrecGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsBP)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/cvode/cvode_bbdpre.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function CVBBDPrecInit(cvode_mem, Nlocal::sunindextype, mudq::sunindextype, mldq::sunindextype, mukeep::sunindextype, mlkeep::sunindextype, dqrely::realtype, gloc::CVLocalFn, cfn::CVCommFn)
    ccall((:CVBBDPrecInit, libsundials_cvode), Cint, (CVODEMemPtr, sunindextype, sunindextype, sunindextype, sunindextype, sunindextype, realtype, CVLocalFn, CVCommFn), cvode_mem, Nlocal, mudq, mldq, mukeep, mlkeep, dqrely, gloc, cfn)
end

function CVBBDPrecReInit(cvode_mem, mudq::sunindextype, mldq::sunindextype, dqrely::realtype)
    ccall((:CVBBDPrecReInit, libsundials_cvode), Cint, (CVODEMemPtr, sunindextype, sunindextype, realtype), cvode_mem, mudq, mldq, dqrely)
end

function CVBBDPrecGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVBBDPrecGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVBBDPrecGetNumGfnEvals(cvode_mem, ngevalsBBDP)
    ccall((:CVBBDPrecGetNumGfnEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, ngevalsBBDP)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/cvode/cvode_diag.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function CVDiag(cvode_mem)
    ccall((:CVDiag, libsundials_cvode), Cint, (CVODEMemPtr,), cvode_mem)
end

function CVDiagGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVDiagGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVDiagGetNumRhsEvals(cvode_mem, nfevalsLS)
    ccall((:CVDiagGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsLS)
end

function CVDiagGetLastFlag(cvode_mem, flag)
    ccall((:CVDiagGetLastFlag, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, flag)
end

function __CVDiagGetReturnFlagName(flag::Clong)
    ccall((:CVDiagGetReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVDiagGetReturnFlagName(flag)
    __CVDiagGetReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/cvode/cvode_direct.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function CVDlsSetLinearSolver(cvode_mem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:CVDlsSetLinearSolver, libsundials_cvode), Cint, (CVODEMemPtr, SUNLinearSolver, SUNMatrix), cvode_mem, LS, A)
end

function CVDlsSetJacFn(cvode_mem, jac::CVDlsJacFn)
    ccall((:CVDlsSetJacFn, libsundials_cvode), Cint, (CVODEMemPtr, CVDlsJacFn), cvode_mem, jac)
end

function CVDlsGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVDlsGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVDlsGetNumJacEvals(cvode_mem, njevals)
    ccall((:CVDlsGetNumJacEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njevals)
end

function CVDlsGetNumRhsEvals(cvode_mem, nfevalsLS)
    ccall((:CVDlsGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsLS)
end

function CVDlsGetLastFlag(cvode_mem, flag)
    ccall((:CVDlsGetLastFlag, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, flag)
end

function __CVDlsGetReturnFlagName(flag::Clong)
    ccall((:CVDlsGetReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVDlsGetReturnFlagName(flag)
    __CVDlsGetReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/cvode/cvode_impl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/cvode/cvode_spils.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function CVSpilsSetLinearSolver(cvode_mem, LS::SUNLinearSolver)
    ccall((:CVSpilsSetLinearSolver, libsundials_cvode), Cint, (CVODEMemPtr, SUNLinearSolver), cvode_mem, LS)
end

function CVSpilsSetEpsLin(cvode_mem, eplifac::realtype)
    ccall((:CVSpilsSetEpsLin, libsundials_cvode), Cint, (CVODEMemPtr, realtype), cvode_mem, eplifac)
end

function CVSpilsSetPreconditioner(cvode_mem, pset::CVSpilsPrecSetupFn, psolve::CVSpilsPrecSolveFn)
    ccall((:CVSpilsSetPreconditioner, libsundials_cvode), Cint, (CVODEMemPtr, CVSpilsPrecSetupFn, CVSpilsPrecSolveFn), cvode_mem, pset, psolve)
end

function CVSpilsSetJacTimes(cvode_mem, jtsetup::CVSpilsJacTimesSetupFn, jtimes::CVSpilsJacTimesVecFn)
    ccall((:CVSpilsSetJacTimes, libsundials_cvode), Cint, (CVODEMemPtr, CVSpilsJacTimesSetupFn, CVSpilsJacTimesVecFn), cvode_mem, jtsetup, jtimes)
end

function CVSpilsGetWorkSpace(cvode_mem, lenrwLS, leniwLS)
    ccall((:CVSpilsGetWorkSpace, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}, Ptr{Clong}), cvode_mem, lenrwLS, leniwLS)
end

function CVSpilsGetNumPrecEvals(cvode_mem, npevals)
    ccall((:CVSpilsGetNumPrecEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, npevals)
end

function CVSpilsGetNumPrecSolves(cvode_mem, npsolves)
    ccall((:CVSpilsGetNumPrecSolves, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, npsolves)
end

function CVSpilsGetNumLinIters(cvode_mem, nliters)
    ccall((:CVSpilsGetNumLinIters, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nliters)
end

function CVSpilsGetNumConvFails(cvode_mem, nlcfails)
    ccall((:CVSpilsGetNumConvFails, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nlcfails)
end

function CVSpilsGetNumJTSetupEvals(cvode_mem, njtsetups)
    ccall((:CVSpilsGetNumJTSetupEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njtsetups)
end

function CVSpilsGetNumJtimesEvals(cvode_mem, njvevals)
    ccall((:CVSpilsGetNumJtimesEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, njvevals)
end

function CVSpilsGetNumRhsEvals(cvode_mem, nfevalsLS)
    ccall((:CVSpilsGetNumRhsEvals, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, nfevalsLS)
end

function CVSpilsGetLastFlag(cvode_mem, flag)
    ccall((:CVSpilsGetLastFlag, libsundials_cvode), Cint, (CVODEMemPtr, Ptr{Clong}), cvode_mem, flag)
end

function __CVSpilsGetReturnFlagName(flag::Clong)
    ccall((:CVSpilsGetReturnFlagName, libsundials_cvode), Cstring, (Clong,), flag)
end

function CVSpilsGetReturnFlagName(flag)
    __CVSpilsGetReturnFlagName(convert(Clong, flag))
end


text:
using Pkg
Pkg.activate(@__DIR__)

# parse command line arguments if given
if length(ARGS) > 0
    subdir = ARGS[1]
# otherwise prompt user to specify
else
    print("Please enter sub directory to run RDE in: ")
    subdir = readline()
end
# input validation
while !isdir(joinpath(@__DIR__, subdir))
    print("Invalid directory $subdir. Please enter sub directory to run RDE in: ")
    global subdir = readline()
end

using PyCall
pushfirst!(PyVector(pyimport("sys")["path"]), joinpath(@__DIR__, subdir))

import FrankWolfe
include("custom_oralces.jl")
include(joinpath(@__DIR__, subdir, "config.jl"))  # load indices, rates, max_iter
cd(subdir)

# Get the Python side of RDE
rde = pyimport("rde")

for idx in indices

    # Load data sample and distortion functional
    x, fname = rde.get_data_sample(idx)
    f, df, node, pred = rde.get_distortion(x)

    # Wrap objective and gradiet functions
    function func(s)
        if !(s isa Vector{eltype(x)})
            s = convert(Vector{eltype(x)}, s)
        end
        return f(s)
    end

    function grad!(storage, s)
        if !(s isa Vector{eltype(x)})
            s = convert(Vector{eltype(x)}, s)
        end
        g = df(s)
        return @. storage = g
    end

    all_s = zeros(eltype(x), (length(rates), length(x)))
    for rate in rates
        # Run FrankWolfe
        println("Running sample $idx with rate $rate")
        s0 = similar(x[:])
        s0 .= 0.0
        lmo = NonNegKSparseLMO(rate, 1.0)

        @time s, v, primal, dual_gap = FrankWolfe.frank_wolfe(
        #@time s, v, primal, dual_gap = FrankWolfe.away_frank_wolfe(
        #@time s, v, primal, dual_gap = FrankWolfe.blended_conditional_gradient(
        #@time s, v, primal, dual_gap = FrankWolfe.lazified_conditional_gradient(
            s -> func(s),
            (storage, s) -> grad!(storage, s),
            lmo,
            s0,
            ;fw_arguments...
        )
        # reset adaptive step size if necessary
        if fw_arguments.line_search isa FrankWolfe.MonotonousNonConvexStepSize
            fw_arguments.line_search.factor = 0
        end

        # Store single rate result
        all_s[indexin(rate, rates)[1], :] = s
        rde.store_single_result(s, idx, fname, rate)

    end

    # Store multiple rate results
    rde.store_collected_results(all_s, idx, node, pred, fname, rates)

end

max_stars_repo_name: jmaces/fw-rde

max_issues_repo_name: ZIB-IOL/fw-rde

max_forks_repo_name: ZIB-IOL/fw-rde


text:
module PriceIndex

greet() = print("Hello World!")

using Distributions, Parameters, LinearAlgebra, Plots, Compat

function FixedLaspeyres(p,q)
N=size(p,1)
T=size(p,2)
Pl=zeros(N+1,T)
s=zeros(N,T)
   for t in 1:T
        for n in 2:N+1
        s[n-1,t] = (p[n-1,t]*q[n-1,t])/dot(p[:,t],q[:,t])
        Pl[n,t]=Pl[n-1,t]+p[n-1,t]/p[n-1,1]*s[n-1,1]
        end
    end
return Pl[N+1,:]
end

function foo(μ = 1., σ = 2.)
    d = Normal(μ, σ)
    E = expectation(d)
    return E(x -> sin(x))
end

export foo, FixedLaspeyres

end # module


text:
# threaded heap knn graph

"""
LockHeapKNNGraph - uses locks to synchronize the heaps that store the underlying
graph edge data. The heaps themselves are *not* thread-safe.
"""
struct LockHeapKNNGraph{V<:Integer,
                        U<:Real,
                        D<:AbstractVector,
                        M<:PreMetric} <: ApproximateKNNGraph{V, U, D, M}
    heaps::Vector{BinaryMaxHeap{HeapKNNGraphEdge{V, U}}}
    locks::Vector{ReentrantLock}
    data::D
    n_neighbors::V
    metric::M
end
function LockHeapKNNGraph(data::D,
                          metric::M,
                          indices::AbstractMatrix{V},
                          distances::AbstractMatrix{U}) where {D <: AbstractVector,
                                                               M <: PreMetric,
                                                               V <: Integer,
                                                               U <: Real}
    n_neighbors = size(indices, 1)
    n_points = size(indices, 2)
    length(data) == size(indices, 2) || error("`indices` must have num columns equal to `length(data)`")
    size(indices) == size(distances) || error("`indices` and `distances` must have same shape")
    n_neighbors < n_points || error("`Must have more columns than rows`")
    # ...
    HeapType = BinaryMaxHeap{HeapKNNGraphEdge{V, U}}
    knn_heaps = HeapType[HeapType() for _ in 1:n_points]
    heap_locks = ReentrantLock[ReentrantLock() for _ in 1:n_points]
    for v in 1:n_points, i in 1:n_neighbors
        _heappush!(knn_heaps[v], HeapKNNGraphEdge(v, indices[i, v], distances[i, v]), n_neighbors)
    end
    return LockHeapKNNGraph{V, U, D, M}(knn_heaps, heap_locks, data, n_neighbors, metric)
end

function LockHeapKNNGraph(data::D, n_neighbors::Integer, metric::M) where {D <: AbstractVector,
                                                                           M <: PreMetric}

    # assert some invariants
    n_neighbors > 0 || error("LockHeapKNNGraph needs positive `n_neighbors`")
    length(data) > n_neighbors || error("LockHeapKNNGraph needs `length(data) > n_neighbors`")

    n_points = length(data)
    U = result_type(metric, first(data), first(data))
    HeapType = BinaryMaxHeap{HeapKNNGraphEdge{typeof(n_neighbors), U}}
    knn_heaps = HeapType[HeapType() for _ in 1:n_points]
    heap_locks = ReentrantLock[ReentrantLock() for _ in 1:n_points]
    # initialize approx knn heaps randomly
    Threads.@threads for i in eachindex(data)
        k_idxs = sample_neighbors(n_points, n_neighbors, exclude=[i])
        for j in k_idxs
            weight = evaluate(metric, data[i], data[j])
            lock(heap_locks[i]) do
                _heappush!(knn_heaps[i], HeapKNNGraphEdge(i, j, weight), n_neighbors)
            end
            if !(metric isa SemiMetric)
                weight = evaluate(metric, data[j], data[i])
            end
            lock(heap_locks[j]) do
                _heappush!(knn_heaps[j], HeapKNNGraphEdge(j, i, weight), n_neighbors)
            end
        end
    end
    return LockHeapKNNGraph{typeof(n_neighbors), U, D, M}(knn_heaps, heap_locks, data, n_neighbors, metric)
end

# convert matrix data input to vector of vectors (views)
function LockHeapKNNGraph(data::D, args...; kwargs...) where D <: AbstractMatrix
    return LockHeapKNNGraph(collect(eachcol(data)), args...; kwargs...)
end

"""
    inneighbors(g::LockHeapKNNGraph, v)

Similar to `inneighbors(g::HeapKNNGraph, v)`, except it acquires locks for the neighbor heaps
as it iterates in order for this to be thread-safe.

"""
@inline function LightGraphs.inneighbors(g::LockHeapKNNGraph{V}, v::V) where V
    neighbs = V[]
    for i in nv(g)
        lock(g.locks[i]) do
            for e in g.heaps[i].valtree
                if dst(e) == v
                    push!(neighbs, src(e))
                end
            end
        end
    end
    return neighbs
end

"""
    outneighbors(g::LockHeapKNNGraph, v)

Similar to `outneighbors(g::HeapKNNGraph, v)`, except locks the neighbor heap before collecting
to make this thread-safe.
"""
@inline function LightGraphs.outneighbors(g::LockHeapKNNGraph{V}, v::V) where V
    neighbs = lock(g.locks[v]) do
        dst.(g.heaps[v].valtree)
    end
    return neighbs
end

"""
    add_edge!(g::LockHeapKNNGraph, e::HeapKNNGraphEdge)

Similar to `add_edge!(g::HeapKNNGraph, e)`, but made thread-safe using locks.
"""
function LightGraphs.add_edge!(g::LockHeapKNNGraph, e::HeapKNNGraphEdge)
    # NOTE we can assume the invariants for heap knn graphs hold
    lock(g.locks[src(e)]) do
        if e < first(g.heaps[src(e)]) && !has_edge(g, src(e), dst(e))
            # we know this edge is smaller than the first, so we can start by removing that
            pop!(g.heaps[src(e)])
            push!(g.heaps[src(e)], e)
            return true
        end
        return false
    end
end

"""
    node_edges(graph::LockHeapKNNGraph, i) -> edges

Return all the outgoing edges from node i in an arbitrary order. Thread-safe.
"""
@inline function node_edges(g::LockHeapKNNGraph{V}, i::V) where V
    lock(g.locks[i]) do
        return edgetype(g)[node_edge(g, i, j) for j in one(V):g.n_neighbors]
    end
end

"""
    update_flag!(g::LockHeapKNNGraph, i, j, flag)

Update the flag of the edge at the given indices. Since the flags don't influence
the edge ordering, this can't invalidate the heap invariant. Uses locks to ensure
thread safety.
"""
function update_flag!(g::LockHeapKNNGraph{V}, i::V, j::V, flag::Bool) where V
    lock(g.locks[i]) do
        edge = node_edge(g, i, j)
        newedge = edgetype(g)(src(edge), dst(edge), weight(edge), flag)
        g.heaps[i].valtree[j] = newedge
        return newedge
    end
end


text:
include("ACQRS.jl")

include("ACQRS_sub_test.jl")


text:
## Define efficient function to check if all elements in array are equal
# method found at https://stackoverflow.com/questions/47564825/check-if-all-the-elements-of-a-julia-array-are-equal/47578613
@inline function _allequal(x::AbstractArray)::Bool
    length(x) < 2 && return true
    e1 = first(x)
    for value in x
        isequal(value, e1) || return false
    end
    return true
end


text:
module PlanningDomains

export load_domain, load_problem
export list_collections, list_domains, list_problems
export find_collections, find_domains, find_problems
export JuliaPlannersRepo, IPCInstancesRepo, PlanningDomainsRepo

import PDDL.Parser: load_domain, load_problem, parse_domain, parse_problem
import Scratch: get_scratch!, delete_scratch!
import Downloads, HTTP, JSON

## Generic types and methods ##

"Abstract type for planning domain and problem repositories."
abstract type PlanningRepository end

const REPOS_DOCSTRING = """
The following repositores are currently supported:

- [`JuliaPlannersRepo`](@ref): Repository maintained by the JuliaPlanners organization
- [`IPCInstancesRepo`](@ref): International Planning Competition domains and problems
- [`PlanningDomainsRepo`](@ref): Repository provided by http://planning.domains/

View the documentation of each repository for more information.
"""

"""
    load_domain(repository, domain)
    load_domain(repository, collection, domain)

Load a domain from a planning `repository`. In some repositories, domains
are also organized into `collection`s.

$REPOS_DOCSTRING
"""
load_domain(::Type{T}, args...) where {T <: PlanningRepository} =
    load_domain(T(), args...)

"""
    load_problem(repository, domain, problem)
    load_problem(repository, collection, domain, problem)

Load a problem (specified as a name or index) from a planning `repository`
for a given `domain`. In some repositories, domains are also organized
into `collection`s.

$REPOS_DOCSTRING
"""
load_problem(::Type{T}, args...) where {T <: PlanningRepository} =
    load_problem(T(), args...)

"""
    list_collections(repository)

List collections in a planning `repository`.

$REPOS_DOCSTRING
"""
list_collections(::Type{T}) where {T <: PlanningRepository} =
    list_collections(T())

"""
    list_domains(repository)
    list_domains(repository, collection)

List domains in a planning `repository`. In some repositories, domains
are also organized into `collection`s.

$REPOS_DOCSTRING
"""
list_domains(::Type{T}, args...) where {T <: PlanningRepository} =
    list_domains(T(), args...)

"""
    list_problems(repository, domain)
    list_problems(repository, collection, domain)

List problems from a planning `repository` for a given `domain`. In some
repositories, domains are also organized into `collection`s.

$REPOS_DOCSTRING
"""
list_problems(::Type{T}, arg, args...) where {T <: PlanningRepository} =
    list_problems(T(), arg, args...)
list_problems(::Type{T}, arg::Symbol, args...) where {T <: PlanningRepository} =
    list_problems(T(), replace(string(arg), '_' => '-'), args...)

"""
    find_collections(repository, query)

Find collections from a planning `repository` that match a `query`, which can
be provided as a regular expression or (sub)string.
"""
function find_collections(repo::PlanningRepository, args...)
    @assert !isempty(args) "Query needs to be specified."
    args, query = args[1:end-1], args[end]
    return filter!(s -> occursin(query, s), list_collections(repo, args...))
end
find_collections(::Type{T}, args...) where {T <: PlanningRepository} =
    find_collections(T(), args...)

"""
    find_domains(repository, query)
    find_domains(repository, collection, query)

Find domains from a planning `repository` that match a `query`, which can be
provided as a regular expression or (sub)string.
"""
function find_domains(repo::PlanningRepository, args...)
    @assert !isempty(args) "Query needs to be specified."
    args, query = args[1:end-1], args[end]
    return filter!(s -> occursin(query, s), list_domains(repo, args...))
end
find_domains(::Type{T}, args...) where {T <: PlanningRepository} =
    find_domains(T(), args...)

"""
    find_problems(repository, domain, query)
    find_problems(repository, collection, domain, query)

Find problems from a planning `repository` in a given `domain` that match
a `query`, which can be provided as a regular expression or (sub)string.
"""
function find_problems(repo::PlanningRepository, args...)
    @assert !isempty(args) "Query needs to be specified."
    args, query = args[1:end-1], args[end]
    return filter!(s -> occursin(query, s), list_problems(repo, args...))
end
find_problems(::Type{T}, args...) where {T <: PlanningRepository} =
    find_problems(T(), args...)

## Repository-specific implementations ##

include("cache_management.jl")
include("julia_planners_repo.jl")
include("ipc_instances_repo.jl")
include("planning_domains_repo.jl")

## Default implementations ##

const DefaultRepo = JuliaPlannersRepo
const Label = Union{AbstractString,Symbol,Int,Regex}

_sym_to_str(s::Symbol) = replace(string(s), '_' => '-')
_sym_to_str(s) = s

"""
    load_domain(domain::Symbol)
    load_domain(collection::Symbol, domain)

Load a domain from the default repository (`$DefaultRepo`). Underscores '_'
in `Symbol`s are automatically converted to dashes '-'.
"""
load_domain(arg::Symbol, args...) =
    load_domain(DefaultRepo(), _sym_to_str(arg), _sym_to_str.(args)...)

"""
    load_problem(domain::Symbol, problem)

Load a problem from the default repository (`$DefaultRepo`). Underscores '_'
in `Symbol`s are automatically converted to dashes '-'.
"""
load_problem(arg::Symbol, args...) =
    load_problem(DefaultRepo(), _sym_to_str(arg), _sym_to_str.(args)...)

"""
    list_collections()

List domains from the default repository (`$DefaultRepo`).
"""
list_collections() =
    list_collections(DefaultRepo())

"""
    list_domains()
    list_domains(collection)

List domains from the default repository (`$DefaultRepo`).
"""
list_domains() =
    list_domains(DefaultRepo())
list_domains(collection) =
    list_domains(DefaultRepo(), _sym_to_str(collection))

"""
    list_problems(domain)
    list_problems(collection, domain)

List problems from the default repository (`$DefaultRepo`).
"""
list_problems(arg, args...) =
    list_problems(DefaultRepo(), _sym_to_str(arg), _sym_to_str.(args)...)

"""
    find_collections(query::Union{AbstractString,Regex})

Find collections from the default repository (`$DefaultRepo`).
"""
find_collections(query::Union{AbstractString,Regex}) =
    find_collections(DefaultRepo(), query)

"""
    find_domains(query::Union{AbstractString,Regex})
    find_domains(collection, query::Union{AbstractString,Regex})

Find domains from the default repository (`$DefaultRepo`).
"""
find_domains(query::Union{AbstractString,Regex}) =
    find_domains(DefaultRepo(), query)
find_domains(collection::Label, query::Union{AbstractString,Regex}) =
    find_domains(DefaultRepo(), _sym_to_str(collection), query)

"""
    find_problems(domain, query::Union{AbstractString,Regex})
    find_problems(collection, domain, query::Union{AbstractString,Regex})

Find problems from the default repository (`$DefaultRepo`).
"""
find_problems(domain::Label, query::Union{AbstractString,Regex}) =
    find_problems(DefaultRepo(), _sym_to_str(domain), query)
find_problems(collection::Label, domain, query::Union{AbstractString,Regex}) =
    find_problems(DefaultRepo(), _sym_to_str(collection),
                  _sym_to_str(domain), query)

end


text:
struct NoLightHouse <: AbstractLightHouse end

struct DistanceGoalLightHouse{T <: Real} <: LossFunctionLightHouse
    startpoint::Complex{T}
    goalpoint::Complex{T}
end

function evolvetransform(lantern::DistanceGoalLightHouse{T}, current_transform::MobiusTransform{T}, playerstate::PlayerStateNoResource) where {T <: Real}
    distfunc = mobius -> complexdistsq(transform(mobius, lantern.startpoint), lantern.goalpoint)
    #@show playerstate
    mapped_orig_player_pos = transform(current_transform, playerstate.preimage_position)

    num_steps = 1
    preimage_player_pos_arr = Vector{Complex{T}}(undef, num_steps + 1)
    preimage_player_pos_arr[1] = playerstate.preimage_position
    lr = 0.005
    for l in 1:num_steps
        distparams = distfunc(current_transform)
        gradmobius = gradient(distfunc, current_transform)[1].x
        gradmobiusSV = SVector(gradmobius...)

        addequals(current_transform, -lr .* gradmobiusSV)
        preimage_player_pos_arr[l + 1] = inversetransform(current_transform, mapped_orig_player_pos)
    end
    playerstate.preimage_position = preimage_player_pos_arr[end]
    current_transform
end

function evolvetransform(lantern::DistanceGoalLightHouse{T}, current_transform::ComplexRational{T, N}, playerstate::PlayerStateNoResource) where {T <: Real, N}
    distfunc = conformal -> complexdistsq(transform(conformal, lantern.startpoint), lantern.goalpoint)
    #@show playerstate
    mapped_orig_player_pos = playerstate.mandel_position

    num_steps = 1
    preimage_player_pos_arr = Vector{Complex{T}}(undef, num_steps + 1)
    preimage_player_pos_arr[1] = playerstate.preimage_position
    lr = 0.005
    for l in 1:num_steps
        distparams = distfunc(current_transform)
        gradmobius = gradient(distfunc, current_transform)[1]
        #@show typeof(gradmobius)
        #@show gradmobius[:numerator], typeof(gradmobius[:numerator])
        #@show current_transform.numerator - (lr .* gradmobius[:numerator])
        @show norm(gradmobius[:numerator]), norm(gradmobius[:denominator])
        current_transform = ComplexRational(current_transform.numerator - (lr .* gradmobius[:numerator]),
            current_transform.denominator - (lr .* gradmobius[:denominator]))
        #gradmobiusSV = SVector(gradmobius...)

        #addequals(current_transform, -lr .* gradmobiusSV)
        preimage_player_pos_arr[l + 1] = inversetransform(current_transform, mapped_orig_player_pos, preimage_player_pos_arr[l])
    end
    playerstate.preimage_position = preimage_player_pos_arr[end]
    current_transform
end

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
# Tests boteSalvatICX, boteSalvatEdgeEnergy and boteSalvatAvailable
using Test
using BoteSalvatICX

@testset "ICX" begin

    # Comparing with values generated by xion.f
    @test isapprox(ionizationcrosssection(12,4,7.71792E+01),3.85315E-18,atol=3e-20)
    @test isapprox(ionizationcrosssection(12,4,1.29569E+04),1.24124E-18,atol=1e-20)
    @test isapprox(ionizationcrosssection(12,4,5.78762E+06),1.30393E-19,atol=1e-21)

    @test isapprox(ionizationcrosssection(23,2,2.90068E+08),5.64142E-21,atol=1e-23)
    @test isapprox(ionizationcrosssection(23,2,1.53993E+04),2.63509E-20,atol=1e-22)

    @test isapprox(ionizationcrosssection(45,7,6.68344E+08),1.62343E-20,atol=1e-22)
    @test isapprox(ionizationcrosssection(45,7,1.22321E+03),2.11571E-19,atol=2e-21)

    @test isapprox(ionizationcrosssection(67,1,1.63117E+05),1.36404E-23,atol=1e-25)
    @test isapprox(ionizationcrosssection(78,9,2.98538E+03),2.02808E-20,atol=1e-22)

    @test isapprox(ionizationcrosssection(99,3,1.00000E+09),1.39024E-22,atol=1e-24)
    @test isapprox(ionizationcrosssection(99,3,2.66073E+04),3.85171E-24,atol=1e-25)

    @test isapprox(ionizationcrosssection(1,1,6.87860E+02),1.57729E-17,atol=1e-19)

end

@testset "EdgeEnergy" begin
    # Not a very good set of values...
    @test isapprox(edgeenergy(67,9),1.36386E+3,atol=1.0)
    @test isapprox(edgeenergy(57,4),5.48366E+3,atol=1.0)
    @test isapprox(edgeenergy(42,1),1.99691E+4,atol=1.0)
    @test isapprox(edgeenergy(42,1),1.99691E+4,atol=1.0)
end

@testset "Available" begin
    @test hasedge(1,1)==true
    @test hasedge(1,2)==false
    @test hasedge(11,2)==true
    @test hasedge(12,4)==true
    @test hasedge(12,5)==false
    @test hasedge(26,6)==true
    @test all(i->hasedge(33,i),1:9)
    @test hasedge(33,10)==false
    @test all(z->hasedge(z,9),33:99)
end


text:
"""
`@layout(d, x)`

Apply the expression `x` to the widget `d`, replacing e.g. symbol `:s` with the corresponding subwidget `d[:s]`
In this context, `_` refers to the whole widget. To use actual symbols, escape them with `^`, as in `^(:a)`.
`@layout` can be combined with `@map` to have the layout update interactively as function of some widget.

## Examples

```jldoctest layout
julia> using Interact

julia> cpt = OrderedDict(:vertical => Observable(true), :b => slider(1:100), :c => button());

julia> t = Widget{:test}(cpt, output = observe(cpt[:b]));

julia> Widgets.@layout t vbox(:b, CSSUtil.vskip(1em), :c);

julia> Widgets.@layout t Widgets.@map &(:vertical) ? vbox(:b, CSSUtil.vskip(1em), :c) : hbox(:b, CSSUtil.hskip(1em), :c);
```

Use [`@layout!`](@ref) to set the widget layout in place:

```jldoctest
julia> @layout! t Widgets.@map &(:vertical) ? vbox(:b, CSSUtil.vskip(1em), :c) : hbox(:b, CSSUtil.hskip(1em), :c);
```

`@layout(x)`

Curried version of `@layout(d, x)`: anonymous function mapping `d` to `@layout(d, x)`.
"""
macro layout(args...)
    esc(layout_helper(args...))
end

function layout_helper(d, expr)
    parse_layout_call(d, expr, replace_wdg)
end

function layout_helper(expr)
    d = gensym()
    Expr(:(->), d, layout_helper(d, expr))
end

replace_wdg(d, x...) = Expr(:call, :(Widgets.component), d, x...)
replace_wdg(s) = s

"""
`@layout!(d, x)`

Set `d.layout` to match the result of `Widgets.@layout(x)`. See [`Widgets.@layout`](@ref) for more information.

## Examples

```jldoctest map
julia> using Interact

julia> t = Widget{:test}(OrderedDict(:b => slider(1:100), :c => button()));

julia> @layout! t hbox(:b, CSSUtil.hskip(1em), :c);
```
"""
macro layout!(args...)
    esc(layout!_helper(args...))
end

function layout!_helper(d, x)
    func = layout_helper(x)
    quote
        $d.layout = $func
        $d
    end
end

function layout!_helper(x)
    d = gensym()
    res = layout!_helper(d, x)
    :($d -> $res)
end

div(args...; kwargs...) = node(:div, args...; kwargs...)

node(::AbstractWidget) = nothing

"""
`node(w::Widget)`

Return primary node for widget `w`
"""
node(w::Widget) = w.scope !== nothing ? w.scope.dom : nothing

function slap_design! end

"""
`layout(f, w::Widget)`

Create a new `Widget` that is a copy of `w` and whose layout is the layout of `w` composed
with the function `f`.

## Examples

```julia
using InteractBase, CSSUtil, Widgets
w = button("OK")
Widgets.layout(w) do t
    hbox("Click here", t)
end
```
"""
function layout(f, w::Widget)
    g = layout(w)
    Widget(w, layout = x -> f(g(x)))
end

"""
`layout(w::Widget)`

Return the function that will be used to determine the layout of widget `w`.
"""
function layout(w::Widget)
    w.layout
end

layout!(f, w::Widget) = (w.layout = f; w)

(w::Widget)(args...; kwargs...) = layout(t->t(args...; kwargs...), w)


text:
using Gadfly, DataArrays, RDatasets

set_default_plot_size(6inch, 3inch)

plot(dataset("datasets", "iris"), x=:SepalLength, y=:SepalWidth,
     color=:Species, Geom.point)


text:
## port methods of series module of SymPy
## http://docs.sympy.org/0.7.2/modules/series.html

## simple methods (x, args) -> y (y coercion happens via PyCall)
series_sympy_meths = (:limit,
                      :gruntz,
                      :series, :Order, :O,
                      :residue
                      )

## fourier series
import Base: scale, truncate
export fourier_series
export shift, shiftx, scalex

fourier_series(ex::Sym, args...; kwargs...) = sympy_meth(:fourier_series, ex, args...; kwargs...)
truncate(ex::Sym; n=3) = object_meth(ex, :truncate, n=n)
shift(ex::Sym, s) = object_meth(ex, :shift, s)
shiftx(ex::Sym, x) = object_meth(ex, :shiftx, x)
scale(ex::Sym, s::AbstractArray) = scale(ex, s[1])
scale(ex::Sym, s) = object_meth(ex, :scale, s)
scalex(ex::Sym, x) = object_meth(ex, :scalex, x)


text:
# Type definitions and interface definitions

"""
    EC{T} <: AbstractArray{T,2}

Wrapper type for elastic constants held as a 6×6 Voigt matrix.  The underlying storage
is a `StaticArrays.MArray`, meaning certain operations using this structure will be
quicker than using normal `Array`s.  For instance, elastic constant transformations
such as rotations, or conversion between stiffness and compliance, are each about
twice as fast.

    EC(data::AbstractArray{T}; warn=false) where T -> ec

Create a set of elastic constants `ec` from the array `data` which must have length 36.
Symmetry is enforced when converting, taking the upper half of the matrix.  No warning
is issued for non-symmetric input, unless `warn` is `true`.
"""
struct EC{T} <: AbstractArray{T,2}
    data::MArray{Tuple{6,6},T,2,36}

    # Arbitrary input requires symmetrisation and conversion to MMatrix if not already
    function EC{T}(data; warn=false) where T
        if data isa AbstractArray
            size(data) == (6,6) || throw(ArgumentError("ECs must have dimensions 6×6"))
        elseif data isa NTuple
            length(data) == 36 || throw(ArgumentError("ECs must have length 36"))
        end
        ec = new{T}(MMatrix{6,6,T}(data))
        _make_symmetric!(ec; warn=warn)
    end
    # Construction with another EC doesn't require symmetrisation
    EC{T}(ec::EC; warn=false) where T = new{T}(ec.data)
end

"Enforce symmetry in a 6×6 array.  No checks on shape performed."
function _make_symmetric!(x; warn=false)
    if warn
        atol = √eps(eltype(x))
        @inbounds for i in 1:6, j in i+1:6
            if !isapprox(x[i,j], x[j,i], atol=atol)
                @warn("input matrix not symmetrical: taking upper half")
                break
            end
        end
    end
    # Copy upper to lower half
    @inbounds for j in 1:6, i in j+1:6
        x[i,j] = x[j,i]
    end
    x
end

EC(x::AbstractArray{T}; kwargs...) where {T<:Number} = EC{float(T)}(x; kwargs...)
EC(x::NTuple{36,T}; kwargs...) where T = EC{float(T)}(x; kwargs...)
EC(x; kwargs...) = EC{DEFAULT_FLOAT}(x; kwargs...)

Base.size(x::EC) = (6,6)
Base.length(x::EC) = 36
Base.getindex(x::EC, i...) = x.data[i...]
Base.setindex!(x::EC, val, i, j) = x.data[i,j] = x.data[j,i] = val
Base.IndexStyle(::Type{EC}) = IndexStyle(MMatrix)

Base.zero(::Type{EC{T}}) where T = EC{T}(zero(MMatrix{6,6,T}))
Base.zero(::Type{EC}) = zero(EC{DEFAULT_FLOAT})
Base.zero(::EC{T}) where T = zero(EC{T})


text:
using Laguerre
using Documenter

makedocs(;
    modules=[Laguerre],
    authors="Zhiwei Peng",
    repo="https://github.com/zpeng2/Laguerre.jl/blob/{commit}{path}#L{line}",
    sitename="Laguerre.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://zpeng2.github.io/Laguerre.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/zpeng2/Laguerre.jl",
)