text:
module RefractiveIndex

using Pkg.Artifacts
using YAML
using Interpolations
using HTTP.URIs: unescapeuri
using Unitful: @u_str, uparse, uconvert, ustrip, AbstractQuantity
using Memoize
using DelimitedFiles: readdlm

import Base: getindex, show

export RefractiveMaterial

const RI_INFO_ROOT = Ref{String}()
const RI_LIB = Dict{Tuple{String, String, String}, NamedTuple{(:name, :path), Tuple{String, String}}}()
include("init.jl")
include("dispersionformulas.jl")

struct RefractiveMaterial{DF<:DispersionFormula}
    name::String
    reference::String
    comment::String
    dispersion::DF
    λrange::Tuple{Float64, Float64}
    specs::Dict{Symbol, Any}
end

const DISPERSIONFORMULAE = Dict(
    "formula 1" => Sellmeier,
    "formula 2" => Sellmeier2,
    "formula 3" => Polynomial,
    "formula 4" => RIInfo,
    "formula 5" => Cauchy,
    "formula 6" => Gases,
    "formula 7" => Herzberger,
    "formula 8" => Retro,
    "formula 9" => Exotic,
    "tabulated nk" => TabulatedNK,
    "tabulated n" => TabulatedN,
    "tabulated k" => TabulatedK,
)

function str2tuple(str)
    arr = parse.(Float64, split(str))
    ntuple(i -> arr[i], length(arr))
end

function DispersionFormula(data)
    DF = DISPERSIONFORMULAE[data[:type]]
    if haskey(data, :coefficients)
        λrange = str2tuple(data[:wavelength_range])
        return DF(str2tuple(data[:coefficients])), λrange
    else
        raw = readdlm(IOBuffer(data[:data]), ' ', Float64)
        λrange = extrema(@view raw[:, 1])
        return DF(raw), λrange
    end
end

function _RM_data(shelf, book, page)
    metadata = RI_LIB[(shelf, book, page)]
    path = joinpath(RI_INFO_ROOT[], "data", metadata.path)
    isfile(path) || @error "Specified material does not exist"
    yaml = YAML.load_file(path; dicttype=Dict{Symbol, Any})
    reference = get(yaml, :REFERENCES, "")
    comment = get(yaml, :COMMENTS, "")
    specs = get(yaml, :SPECS, Dict{Symbol, Any}())
    data = only(get(yaml, :DATA, Dict{Symbol, String}[]))
end
"""
    RefractiveMaterial(shelf, book, page)

Load the refractive index data for the material corresponding to the specified
shelf, book, and page within the [refractiveindex.info](https://refractiveindex.info/) database. The data
can be queried by calling the returned `RefractiveMaterial` object at a given wavelength.

# Examples
```julia-repl
julia> MgLiTaO3 = RefractiveMaterial("other", "Mg-LiTaO3", "Moutzouris-o")
"Mg-LiTaO3 (Moutzouris et al. 2011: n(o) 0.450-1.551 µm; 8 mol.% Mg)"

julia> MgLiTaO3(0.45) # default unit is microns
2.2373000025056826

julia> using Unitful

julia> MgLiTaO3(450u"nm") # auto-conversion from generic Unitful.jl length units
2.2373000025056826

julia> MgLiTaO3(450e-9, "m") # strings can be used to specify units (parsing is cached)
2.2373000025056826
```
"""
function RefractiveMaterial(shelf, book, page)
    metadata = RI_LIB[(shelf, book, page)]
    path = joinpath(RI_INFO_ROOT[], "data", metadata.path)
    isfile(path) || @error "Specified material does not exist"
    yaml = YAML.load_file(path; dicttype=Dict{Symbol, Any})
    reference = get(yaml, :REFERENCES, "")
    comment = get(yaml, :COMMENTS, "")
    specs = get(yaml, :SPECS, Dict{Symbol, Any}())
    data = only(get(yaml, :DATA, Dict{Symbol, String}[]))
    DF, λrange = DispersionFormula(data)

    RefractiveMaterial(
        string(book, " ($(metadata.name))"),
        reference,
        comment,
        DF,
        λrange,
        specs
    )
end

"""
    RefractiveMaterial(url::String)

Extracts the shelf, book, and page from a refractiveindex.info URL and loads
the corresponding data from the local database (does not require an active internet connection).

!!! warning
    The refractiveindex.info website is regularly updated and may contain materials not yet
    available in the local copy of the database, which is updated on a roughly annual basis.
    Future versions of this package may allow these new entries to be automatically downloaded
    on demand.

# Examples
```julia-repl
julia> Ar = RefractiveMaterial("https://refractiveindex.info/?shelf=main&book=Ar&page=Peck-15C")
"Ar (Peck and Fisher 1964: n 0.47-2.06 µm; 15 °C)"

julia> Ar(532, "nm")
1.0002679711455778
```
"""
function RefractiveMaterial(url::String)
    ue_url = unescapeuri(url)
    r = r"refractiveindex.info\/\?shelf=(?'shelf'\w+)&book=(?'book'.*)&page=(?'page'.*)"
    m = match(r, ue_url)
    isnothing(m) && @error "Invalid refractiveindex.info url"
    RefractiveMaterial(String(m["shelf"]),
                       String(m["book"]),
                       String(m["page"]))
end

show(io::IO, ::MIME"text/plain", m::RefractiveMaterial{DF}) where {DF} = show(io, m.name)
(m::RefractiveMaterial)(λ::Float64) = m.dispersion(λ)
(m::RefractiveMaterial)(λ::AbstractQuantity) = m(Float64(ustrip(uconvert(u"μm", λ))))

@memoize _dim_to_micron(dim) = ustrip(uconvert(u"μm", 1.0uparse(dim)))
(m::RefractiveMaterial)(λ, dim::String) = m(λ*_dim_to_micron(dim))

(m::RefractiveMaterial{T})(λ::Float64) where {T <: Tabulated}= m.dispersion.n(λ)
end # module


text:
#example of creating chains in flux
using Flux, BenchmarkTools

model = Chain(x -> 3x^2, x -> 2x, x -> x+1 )

# now take a gradient of this model
@btime df(x) = Flux.gradient(model,x)[1]

@show model(5)
@show df(5)

ext: jl


text:
include("common.jl")
include("compression.jl")
include("superchunks.jl")


text:
using Revise

using OpticalBlochEquations

using Plots
using SparseArrays
using LinearAlgebra
using StaticArrays
using BenchmarkTools
using PhysicalConstants.CODATA2018
import DifferentialEquations: ODEProblem, solve, DP5

const h     = PlanckConstant.val
const ħ     = h / 2π
const μB    = BohrMagneton.val
const c     = SpeedOfLightInVacuum.val
const ϵ0    = VacuumElectricPermittivity.val

const λ     = @with_unit 626 "nm"
const Γ     = @with_unit 2π * 8.3 "MHz"
const M     = @with_unit 50 "u"
const E₀    = c / λ

# Generate states
m1 = manifold(F=2, ω=0, μ=0)
m2 = manifold(F=3, ω=E₀/Γ, μ=0, Γ=Γ)
states = [m1.states..., m2.states...]

δ = -2.73Γ

# Generate lasers
l1 = define_laser([1,0,0], [0,0,1], (E₀ + δ)/Γ, 4*1.352464)

lasers = [l1]

d = zeros(12, 12, 3)
d[1, 8, 1] = 0.25819889
d[2, 9, 1] = 0.4472136
d[3, 10, 1] = 0.63245553
d[4, 11, 1] = 0.81649658
d[5, 12, 1] = 1

d[1, 7, 2] = -0.57735027
d[2, 8, 2] = -0.73029674
d[3, 9, 2] = -0.77459667
d[4, 10, 2] = -0.73029674
d[5, 11, 2] = -0.57735027

d[1, 6, 3] = 1
d[2, 7, 3] = 0.81649658
d[3, 8, 3] = 0.63245553
d[4, 9, 3] = 0.4472136
d[5, 10, 3]= 0.25819889

for q in 1:3, i in 1:5, j in 6:12
    d[j,i,q] = d[i,j,q]
end

(dρ, ρ, p) = obe(states, lasers, d)
@btime ρ!(dρ, ρ, p, 1.0)

tspan = (0., 600); tstops = 0:1:600
prob = ODEProblem(ρ!, ρ, tspan, p)
@time sol = solve(prob, alg=DP5(), save_everystep=true)

ρ₁₁ = [real(x[1,1]) for x in sol.u]
ρ₂₂ = [real(x[2,2]) for x in sol.u]
ρ₃₃ = [real(x[3,3]) for x in sol.u]
ρ₄₄ = [real(x[4,4]) for x in sol.u]
ρ₅₅ = [real(x[5,5]) for x in sol.u]
plot(sol.t, ρ₁₁)
plot!(sol.t, ρ₂₂)
plot!(sol.t, ρ₃₃)
plot!(sol.t, ρ₄₄)
plot!(sol.t, ρ₅₅, yticks=0.0:0.2:1.0)


text:
module DashQuill
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.0.3"

include("jl/'dre'_quill.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_quill",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "async-Quill.js",
    external_url = "https://unpkg.com/dash_quill@0.0.3/dash_quill/async-Quill.js",
    dynamic = nothing,
    async = :true,
    type = :js
),
DashBase.Resource(
    relative_package_path = "async-Quill.js.map",
    external_url = "https://unpkg.com/dash_quill@0.0.3/dash_quill/async-Quill.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_quill.min.js",
    external_url = "https://unpkg.com/dash_quill@0.0.3/dash_quill/dash_quill.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_quill.min.js.map",
    external_url = "https://unpkg.com/dash_quill@0.0.3/dash_quill/dash_quill.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

max_stars_repo_name: pandamodium/dash-quill

max_issues_repo_name: pandamodium/dash-quill

max_forks_repo_name: pandamodium/dash-quill


text:
using ZygoteRules, Test

foo(x) = 2x
@adjoint foo(x) = 3x, ȳ -> (3,)

# using Zygote

# @test gradient(foo, x) == (3,)


text:
function move_left(b::Buffer)
    b.cursor.pos[2] -= parse_n(b)
    after(b)
end
function move_right(b::Buffer)
    b.cursor.pos[2] += parse_n(b)
    after(b)
end
function move_up(b::Buffer)
    b.cursor.pos[1] -= parse_n(b)
    after(b)
end
function move_down(b::Buffer)
    b.cursor.pos[1] += parse_n(b)
    after(b)
end
function move_eol(b::Buffer)
    b.cursor.pos[2] = width(b)
    after(b)
end
function move_bol(b::Buffer)
    b.cursor.pos[2] = 1
    after(b)
end

function go(b::Buffer)
    b.mode[2] = b.mode[1]
    setmode(b, go_mode)
end
function gobottom(b::Buffer)
    b.cursor.pos[1] = height(b)
    after(b)
    escape(b)
end

function enter_findmode(b::Buffer)
    b.mode[2] = b.mode[1]
    setmode(b, find_mode)
end

function enter_searchmode(b::Buffer)
    b.mode[2] = b.mode[1]
    b.state[:search] = ""
    setmode(b, search_mode)
end

function next(b::Buffer)
    search(b)
end

#Caution: This doesn't work the same as it does in vim
#In vim, an empty line is included with w but not e, these are different operators
#Here, we treat the e operation as a combination of e and an eow call
function move_word(b::Buffer)
    n = parse_n(b.args)
    move_to_nth_word(b, n)
    after(b)
end
function move_eow(b::Buffer)
    n = parse_n(b.args)
    move_to_nth_eow(b, n)
    after(b)
end
function back_word(b::Buffer)
    n = parse_n(b.args)
    back_to_nth_word(b, n)
    after(b)
end
function back_eow(b::Buffer)
    n = parse_n(b.args)
    back_to_nth_eow(b, n)
    after(b)
end

function move_Word(b::Buffer)
    n = parse_n(b.args)
    move_to_nth_Word(b, n)
    after(b)
end
function move_eoW(b::Buffer)
    n = parse_n(b.args)
    move_to_nth_eoW(b, n)
    after(b)
end
function back_Word(b::Buffer)
    n = parse_n(b.args)
    back_to_nth_Word(b, n)
    after(b)
end
function back_eoW(b::Buffer)
    n = parse_n(b.args)
    back_to_nth_eoW(b)
    after(b)
end

movements = Dict('h' => move_left,
                 'j' => move_down,
                 'k' => move_up,
                 'l' => move_right,
                 '$' => move_eol,
                 'g' => go,
                 'G' => gobottom,
                 'f' => enter_findmode,
                 '/' => enter_searchmode,
                 'n' => next,
                 'w' => move_word,
                 'e' => move_eow,
                 'b' => back_eow,
                 'W' => move_Word,
                 'E' => move_eoW,
                 'B' => back_eoW,
                )


text:
module MLverse

using Reexport
@reexport using LearnBase
@reexport using LossFunctions
@reexport using PenaltyFunctions
@reexport using LearningStrategies
@reexport using DataGenerator

end


text:
# Julia wrapper for header: IviScope.h
# Automatically generated using Clang.jl

function IviScope_init(logicalName, idQuery, resetDevice, vi)
    ccall((:IviScope_init, libniscope), ViStatus, (ViRsrc, ViBoolean, ViBoolean, Ptr{ViSession}), logicalName, idQuery, resetDevice, vi)
end

function IviScope_close(vi)
    ccall((:IviScope_close, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_reset(vi)
    ccall((:IviScope_reset, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_self_test(vi, selfTestResult, selfTestMessage)
    ccall((:IviScope_self_test, libniscope), ViStatus, (ViSession, Ptr{ViInt16}, Ptr{ViChar}), vi, selfTestResult, selfTestMessage)
end

function IviScope_error_query(vi, errorCode, errorMessage)
    ccall((:IviScope_error_query, libniscope), ViStatus, (ViSession, Ptr{ViInt32}, Ptr{ViChar}), vi, errorCode, errorMessage)
end

function IviScope_error_message(vi, statusCode, message)
    ccall((:IviScope_error_message, libniscope), ViStatus, (ViSession, ViStatus, Ptr{ViChar}), vi, statusCode, message)
end

function IviScope_revision_query(vi, driverRev, instrRev)
    ccall((:IviScope_revision_query, libniscope), ViStatus, (ViSession, Ptr{ViChar}, Ptr{ViChar}), vi, driverRev, instrRev)
end

function IviScope_InvalidateAllAttributes(vi)
    ccall((:IviScope_InvalidateAllAttributes, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_ResetWithDefaults(vi)
    ccall((:IviScope_ResetWithDefaults, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_Disable(vi)
    ccall((:IviScope_Disable, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_InitWithOptions(logicalName, IDQuery, resetDevice, optionString, vi)
    ccall((:IviScope_InitWithOptions, libniscope), ViStatus, (ViRsrc, ViBoolean, ViBoolean, ViConstString, Ptr{ViSession}), logicalName, IDQuery, resetDevice, optionString, vi)
end

function IviScope_GetAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviScope_GetAttributeViInt32, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViInt32}), vi, channelName, attributeId, value)
end

function IviScope_SetAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviScope_SetAttributeViInt32, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32), vi, channelName, attributeId, value)
end

function IviScope_CheckAttributeViInt32(vi, channelName, attributeId, value)
    ccall((:IviScope_CheckAttributeViInt32, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32), vi, channelName, attributeId, value)
end

function IviScope_GetAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviScope_GetAttributeViInt64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViInt64}), vi, channelName, attributeId, value)
end

function IviScope_SetAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviScope_SetAttributeViInt64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt64), vi, channelName, attributeId, value)
end

function IviScope_CheckAttributeViInt64(vi, channelName, attributeId, value)
    ccall((:IviScope_CheckAttributeViInt64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt64), vi, channelName, attributeId, value)
end

function IviScope_GetAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviScope_GetAttributeViReal64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViReal64}), vi, channelName, attributeId, value)
end

function IviScope_SetAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviScope_SetAttributeViReal64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViReal64), vi, channelName, attributeId, value)
end

function IviScope_CheckAttributeViReal64(vi, channelName, attributeId, value)
    ccall((:IviScope_CheckAttributeViReal64, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViReal64), vi, channelName, attributeId, value)
end

function IviScope_GetAttributeViString(vi, channelName, attributeId, bufferSize, value)
    ccall((:IviScope_GetAttributeViString, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViInt32, Ptr{ViChar}), vi, channelName, attributeId, bufferSize, value)
end

function IviScope_SetAttributeViString(vi, channelName, attributeId, value)
    ccall((:IviScope_SetAttributeViString, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViConstString), vi, channelName, attributeId, value)
end

function IviScope_CheckAttributeViString(vi, channelName, attributeId, value)
    ccall((:IviScope_CheckAttributeViString, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViConstString), vi, channelName, attributeId, value)
end

function IviScope_GetAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviScope_GetAttributeViBoolean, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViBoolean}), vi, channelName, attributeId, value)
end

function IviScope_SetAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviScope_SetAttributeViBoolean, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViBoolean), vi, channelName, attributeId, value)
end

function IviScope_CheckAttributeViBoolean(vi, channelName, attributeId, value)
    ccall((:IviScope_CheckAttributeViBoolean, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViBoolean), vi, channelName, attributeId, value)
end

function IviScope_GetAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviScope_GetAttributeViSession, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, Ptr{ViSession}), vi, channelName, attributeId, value)
end

function IviScope_SetAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviScope_SetAttributeViSession, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViSession), vi, channelName, attributeId, value)
end

function IviScope_CheckAttributeViSession(vi, channelName, attributeId, value)
    ccall((:IviScope_CheckAttributeViSession, libniscope), ViStatus, (ViSession, ViConstString, ViAttr, ViSession), vi, channelName, attributeId, value)
end

function IviScope_LockSession(vi, callerHasLock)
    ccall((:IviScope_LockSession, libniscope), ViStatus, (ViSession, Ptr{ViBoolean}), vi, callerHasLock)
end

function IviScope_UnlockSession(vi, callerHasLock)
    ccall((:IviScope_UnlockSession, libniscope), ViStatus, (ViSession, Ptr{ViBoolean}), vi, callerHasLock)
end

function IviScope_GetError(vi, errorCode, bufferSize, description)
    ccall((:IviScope_GetError, libniscope), ViStatus, (ViSession, Ptr{ViStatus}, ViInt32, Ptr{ViChar}), vi, errorCode, bufferSize, description)
end

function IviScope_ClearError(vi)
    ccall((:IviScope_ClearError, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_GetNextInterchangeWarning(vi, bufferSize, warning)
    ccall((:IviScope_GetNextInterchangeWarning, libniscope), ViStatus, (ViSession, ViInt32, Ptr{ViChar}), vi, bufferSize, warning)
end

function IviScope_ClearInterchangeWarnings(vi)
    ccall((:IviScope_ClearInterchangeWarnings, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_ResetInterchangeCheck(vi)
    ccall((:IviScope_ResetInterchangeCheck, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_GetNextCoercionRecord(vi, bufferSize, record)
    ccall((:IviScope_GetNextCoercionRecord, libniscope), ViStatus, (ViSession, ViInt32, Ptr{ViChar}), vi, bufferSize, record)
end

function IviScope_GetSpecificDriverCHandle(vi, specificDriverCHandle)
    ccall((:IviScope_GetSpecificDriverCHandle, libniscope), ViStatus, (ViSession, Ptr{ViSession}), vi, specificDriverCHandle)
end

function IviScope_GetSpecificDriverIUnknownPtr(vi, specificDriverIUnknownPtr)
    ccall((:IviScope_GetSpecificDriverIUnknownPtr, libniscope), ViStatus, (ViSession, Ptr{Cvoid}), vi, specificDriverIUnknownPtr)
end

function IviScope_ConfigureChannel(vi, channel, vrange, offset, coupling, probeAttenuation, enabled)
    ccall((:IviScope_ConfigureChannel, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViInt32, ViReal64, ViBoolean), vi, channel, vrange, offset, coupling, probeAttenuation, enabled)
end

function IviScope_ConfigureChanCharacteristics(vi, channel, inputImpedance, maxInputFrequency)
    ccall((:IviScope_ConfigureChanCharacteristics, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64), vi, channel, inputImpedance, maxInputFrequency)
end

function IviScope_ConfigureAcquisitionType(vi, acquisitionType)
    ccall((:IviScope_ConfigureAcquisitionType, libniscope), ViStatus, (ViSession, ViInt32), vi, acquisitionType)
end

function IviScope_ConfigureAcquisitionRecord(vi, timePerRecord, minimumRecordLength, acqStartTime)
    ccall((:IviScope_ConfigureAcquisitionRecord, libniscope), ViStatus, (ViSession, ViReal64, ViInt32, ViReal64), vi, timePerRecord, minimumRecordLength, acqStartTime)
end

function IviScope_ActualRecordLength(vi, actualRecordLength)
    ccall((:IviScope_ActualRecordLength, libniscope), ViStatus, (ViSession, Ptr{ViInt32}), vi, actualRecordLength)
end

function IviScope_SampleRate(vi, sampleRate)
    ccall((:IviScope_SampleRate, libniscope), ViStatus, (ViSession, Ptr{ViReal64}), vi, sampleRate)
end

function IviScope_ConfigureTrigger(vi, triggerType, holdoff)
    ccall((:IviScope_ConfigureTrigger, libniscope), ViStatus, (ViSession, ViInt32, ViReal64), vi, triggerType, holdoff)
end

function IviScope_ConfigureTriggerCoupling(vi, coupling)
    ccall((:IviScope_ConfigureTriggerCoupling, libniscope), ViStatus, (ViSession, ViInt32), vi, coupling)
end

function IviScope_ConfigureEdgeTriggerSource(vi, source, level, slope)
    ccall((:IviScope_ConfigureEdgeTriggerSource, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViInt32), vi, source, level, slope)
end

function IviScope_ReadWaveform(vi, channel, waveformSize, maxTime, waveform, actualPoints, initialX, xIncrement)
    ccall((:IviScope_ReadWaveform, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, ViInt32, Ptr{ViReal64}, Ptr{ViInt32}, Ptr{ViReal64}, Ptr{ViReal64}), vi, channel, waveformSize, maxTime, waveform, actualPoints, initialX, xIncrement)
end

function IviScope_Abort(vi)
    ccall((:IviScope_Abort, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_InitiateAcquisition(vi)
    ccall((:IviScope_InitiateAcquisition, libniscope), ViStatus, (ViSession,), vi)
end

function IviScope_AcquisitionStatus(vi, status)
    ccall((:IviScope_AcquisitionStatus, libniscope), ViStatus, (ViSession, Ptr{ViInt32}), vi, status)
end

function IviScope_FetchWaveform(vi, channel, waveformSize, waveform, actualPoints, initialX, xIncrement)
    ccall((:IviScope_FetchWaveform, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{ViReal64}, Ptr{ViInt32}, Ptr{ViReal64}, Ptr{ViReal64}), vi, channel, waveformSize, waveform, actualPoints, initialX, xIncrement)
end

function IviScope_IsInvalidWfmElement(vi, elementValue, isInvalid)
    ccall((:IviScope_IsInvalidWfmElement, libniscope), ViStatus, (ViSession, ViReal64, Ptr{ViBoolean}), vi, elementValue, isInvalid)
end

function IviScope_GetChannelName(vi, index, bufferSize, name)
    ccall((:IviScope_GetChannelName, libniscope), ViStatus, (ViSession, ViInt32, ViInt32, Ptr{ViChar}), vi, index, bufferSize, name)
end

function IviScope_ConfigureTVTriggerSource(vi, source, TVSignalFormat, TVEvent, TVPolarity)
    ccall((:IviScope_ConfigureTVTriggerSource, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, ViInt32, ViInt32), vi, source, TVSignalFormat, TVEvent, TVPolarity)
end

function IviScope_ConfigureTVTriggerLineNumber(vi, lineNumber)
    ccall((:IviScope_ConfigureTVTriggerLineNumber, libniscope), ViStatus, (ViSession, ViInt32), vi, lineNumber)
end

function IviScope_ConfigureRuntTriggerSource(vi, source, runtLowThreshold, runtHighThreshold, runtPolarity)
    ccall((:IviScope_ConfigureRuntTriggerSource, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViInt32), vi, source, runtLowThreshold, runtHighThreshold, runtPolarity)
end

function IviScope_ConfigureGlitchTriggerSource(vi, source, level, glitchWidth, glitchPolarity, glitchCondition)
    ccall((:IviScope_ConfigureGlitchTriggerSource, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViInt32, ViInt32), vi, source, level, glitchWidth, glitchPolarity, glitchCondition)
end

function IviScope_ConfigureWidthTriggerSource(vi, source, level, widthLowThreshold, widthHighThreshold, widthPolarity, widthCondition)
    ccall((:IviScope_ConfigureWidthTriggerSource, libniscope), ViStatus, (ViSession, ViConstString, ViReal64, ViReal64, ViReal64, ViInt32, ViInt32), vi, source, level, widthLowThreshold, widthHighThreshold, widthPolarity, widthCondition)
end

function IviScope_ConfigureAcLineTriggerSlope(vi, slope)
    ccall((:IviScope_ConfigureAcLineTriggerSlope, libniscope), ViStatus, (ViSession, ViInt32), vi, slope)
end

function IviScope_ConfigureTriggerModifier(vi, triggerModifier)
    ccall((:IviScope_ConfigureTriggerModifier, libniscope), ViStatus, (ViSession, ViInt32), vi, triggerModifier)
end

function IviScope_ConfigureNumEnvelopes(vi, numberOfEnvelopes)
    ccall((:IviScope_ConfigureNumEnvelopes, libniscope), ViStatus, (ViSession, ViInt32), vi, numberOfEnvelopes)
end

function IviScope_ReadMinMaxWaveform(vi, channel, waveformSize, maxTime, minWaveform, maxWaveform, actualPoints, initialX, xIncrement)
    ccall((:IviScope_ReadMinMaxWaveform, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, ViInt32, Ptr{ViReal64}, Ptr{ViReal64}, Ptr{ViInt32}, Ptr{ViReal64}, Ptr{ViReal64}), vi, channel, waveformSize, maxTime, minWaveform, maxWaveform, actualPoints, initialX, xIncrement)
end

function IviScope_FetchMinMaxWaveform(vi, channel, waveformSize, minWaveform, maxWaveform, actualPoints, initialX, xIncrement)
    ccall((:IviScope_FetchMinMaxWaveform, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{ViReal64}, Ptr{ViReal64}, Ptr{ViInt32}, Ptr{ViReal64}, Ptr{ViReal64}), vi, channel, waveformSize, minWaveform, maxWaveform, actualPoints, initialX, xIncrement)
end

function IviScope_ConfigureRefLevels(vi, lowRef, midRef, highRef)
    ccall((:IviScope_ConfigureRefLevels, libniscope), ViStatus, (ViSession, ViReal64, ViReal64, ViReal64), vi, lowRef, midRef, highRef)
end

function IviScope_ReadWaveformMeasurement(vi, channel, measurementFunction, maxTime, measurement)
    ccall((:IviScope_ReadWaveformMeasurement, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, ViInt32, Ptr{ViReal64}), vi, channel, measurementFunction, maxTime, measurement)
end

function IviScope_FetchWaveformMeasurement(vi, channel, measurementFunction, measurement)
    ccall((:IviScope_FetchWaveformMeasurement, libniscope), ViStatus, (ViSession, ViConstString, ViInt32, Ptr{ViReal64}), vi, channel, measurementFunction, measurement)
end

function IviScope_ConfigureNumAverages(vi, numberOfAverages)
    ccall((:IviScope_ConfigureNumAverages, libniscope), ViStatus, (ViSession, ViInt32), vi, numberOfAverages)
end

function IviScope_ConfigureInitiateContinuous(vi, continuousAcquisition)
    ccall((:IviScope_ConfigureInitiateContinuous, libniscope), ViStatus, (ViSession, ViBoolean), vi, continuousAcquisition)
end

function IviScope_ConfigureInterpolation(vi, interpolation)
    ccall((:IviScope_ConfigureInterpolation, libniscope), ViStatus, (ViSession, ViInt32), vi, interpolation)
end

function IviScope_SampleMode(vi, sampleMode)
    ccall((:IviScope_SampleMode, libniscope), ViStatus, (ViSession, Ptr{ViInt32}), vi, sampleMode)
end

function IviScope_AutoProbeSenseValue(vi, channel, autoProbeSenseValue)
    ccall((:IviScope_AutoProbeSenseValue, libniscope), ViStatus, (ViSession, ViConstString, Ptr{ViReal64}), vi, channel, autoProbeSenseValue)
end

function IviScope_AutoSetup(vi)
    ccall((:IviScope_AutoSetup, libniscope), ViStatus, (ViSession,), vi)
end


text:
@symbol_func function alpha_n_star(cur_reactor::AbstractReactor)
  cur_scaling = cur_reactor.mode_scaling

  cur_alpha_star = 1.0

  cur_alpha_star += cur_scaling[:n]

  cur_alpha_star -= 2 * cur_scaling[:P]

  cur_alpha_star
end


text:
module Metadata

import Dates

using FlatBuffers

@enum(DType, BOOL = 0, INT8 = 1, INT16 = 2, INT32 = 3, INT64 = 4,
  UINT8 = 5, UINT16 = 6, UINT32 = 7, UINT64 = 8,
  FLOAT = 9, DOUBLE = 10,  UTF8 = 11,  BINARY = 12,
  CATEGORY = 13, TIMESTAMP = 14, DATE = 15, TIME = 16)

@enum(Encoding, PLAIN = 0, DICTIONARY = 1)

@enum(TimeUnit, SECOND = 0, MILLISECOND = 1, MICROSECOND = 2, NANOSECOND = 3)

# FlatBuffers.enumsizeof(::Type{TimeUnit}) = UInt8

mutable struct PrimitiveArray
    dtype::DType
    encoding::Encoding
    offset::Int64
    length::Int64
    null_count::Int64
    total_bytes::Int64
end

mutable struct CategoryMetadata
    levels::PrimitiveArray
    ordered::Bool
end

@DEFAULT CategoryMetadata ordered=false

mutable struct TimestampMetadata
    unit::TimeUnit
    timezone::String
end

mutable struct DateMetadata
end

mutable struct TimeMetadata
    unit::TimeUnit
end

@UNION TypeMetadata (Nothing,CategoryMetadata,TimestampMetadata,DateMetadata,TimeMetadata)

mutable struct Column
    name::String
    values::PrimitiveArray
    metadata_type::Int8
    metadata::TypeMetadata
    user_metadata::String
end

function Column(name::String, values::PrimitiveArray, metadata::TypeMetadata=nothing,
                user_metadata::String="")
    Column(name, values, FlatBuffers.typeorder(TypeMetadata, typeof(metadata)),
           metadata, user_metadata)
end

mutable struct CTable
    description::String
    num_rows::Int64
    columns::Vector{Column}
    version::Int32
    metadata::String
end

end # module

# wesm/feather/cpp/src/metadata_generated.h
# wesm/feather/cpp/src/types.h
const JULIA_TYPE_DICT = Dict{Metadata.DType,DataType}(
    Metadata.BOOL      => Bool,
    Metadata.INT8      => Int8,
    Metadata.INT16     => Int16,
    Metadata.INT32     => Int32,
    Metadata.INT64     => Int64,
    Metadata.UINT8     => UInt8,
    Metadata.UINT16    => UInt16,
    Metadata.UINT32    => UInt32,
    Metadata.UINT64    => UInt64,
    Metadata.FLOAT     => Float32,
    Metadata.DOUBLE    => Float64,
    Metadata.UTF8      => String,  # can also be WeakRefString{UInt8}
    Metadata.BINARY    => Vector{UInt8},
    Metadata.CATEGORY  => Int64,
    Metadata.TIMESTAMP => Int64,
    Metadata.DATE      => Int64,
    Metadata.TIME      => Int64
)

const METADATA_TYPE_DICT = Dict{DataType,Metadata.DType}(
    Bool    => Metadata.BOOL,
    Int8    => Metadata.INT8,
    Int16   => Metadata.INT16,
    Int32   => Metadata.INT32,
    Int64   => Metadata.INT64,
    UInt8   => Metadata.UINT8,
    UInt16  => Metadata.UINT16,
    UInt32  => Metadata.UINT32,
    UInt64  => Metadata.UINT64,
    Float32 => Metadata.FLOAT,
    Float64 => Metadata.DOUBLE,
    String  => Metadata.UTF8,
    Vector{UInt8} => Metadata.BINARY,
    Dates.Time => Metadata.INT64,
    Dates.DateTime => Metadata.INT64,
    Dates.Date => Metadata.INT32,
    # WeakRefString{UInt8} => Metadata.UTF8  # not currently being used
)

const NON_PRIMITIVE_TYPES = Set([Metadata.UTF8, Metadata.BINARY])

const JULIA_TIME_DICT = Dict{Metadata.TimeUnit,DataType}(
    Metadata.SECOND => Dates.Second,
    Metadata.MILLISECOND => Dates.Millisecond,
    Metadata.MICROSECOND => Dates.Microsecond,
    Metadata.NANOSECOND => Dates.Nanosecond
)
const METADATA_TIME_DICT = Dict{DataType,Metadata.TimeUnit}(v=>k for (k,v) in JULIA_TIME_DICT)

isprimitivetype(t::Metadata.DType) = t ∉ NON_PRIMITIVE_TYPES

juliatype(meta::Nothing, values_type::Metadata.DType) = JULIA_TYPE_DICT[values_type]
juliatype(values_type::Metadata.DType) = juliatype(nothing, values_type)
function juliatype(meta::Metadata.CategoryMetadata, values_type::Metadata.DType)
    JULIA_TYPE_DICT[meta.levels.dtype]
end
function juliatype(meta::Metadata.TimestampMetadata, values_type::Metadata.DType)
    Timestamp{JULIA_TIME_DICT[meta.unit]}
end
function juliatype(meta::Metadata.TimeMetadata, values_type::Metadata.DType)
    TimeOfDay{JULIA_TIME_DICT[meta.unit],JULIA_TYPE_DICT[values_type]}
end
juliatype(meta::Metadata.DateMetadata, values_type::Metadata.DType) = Datestamp

function juliatype(col::Metadata.Column)
    T = juliatype(col.metadata, col.values.dtype)
    col.values.null_count == 0 ? T : Union{T,Missing}
end

function feathertype(::Type{T}) where T
    if !haskey(METADATA_TYPE_DICT, T)
        throw(ArgumentError("Type $T is not supported by the Feather format."))
    else
        return METADATA_TYPE_DICT[T]
    end
end
feathertype(::Type{Union{T,Missing}}) where T = feathertype(T)
feathertype(::Type{<:Arrow.Datestamp}) = Metadata.INT32
feathertype(::Type{<:Arrow.Timestamp}) = Metadata.INT64
feathertype(::Type{<:Arrow.TimeOfDay{P,Int32}}) where P = Metadata.INT32
feathertype(::Type{<:Arrow.TimeOfDay{P,Int64}}) where P = Metadata.INT64

getmetadata(io::IO, ::Type{T}, A::ArrowVector) where T = nothing
getmetadata(io::IO, ::Type{Union{T,Missing}}, A::ArrowVector) where T = getmetadata(io, T, A)
getmetadata(io::IO, ::Type{Arrow.Datestamp}, A::ArrowVector) = Metadata.DateMetadata()
function getmetadata(io::IO, ::Type{Arrow.Timestamp{T}}, A::ArrowVector) where T
    Metadata.TimestampMetadata(METADATA_TIME_DICT[T], "")
end
function getmetadata(io::IO, ::Type{Arrow.TimeOfDay{P,T}}, A::ArrowVector) where {P,T}
    Metadata.TimeMetadata(METADATA_TIME_DICT[P])
end
# WARNINGMDATA_TYPE_DICT Arrow standard says nothing about specifying whether DictEncoding is ordered!
function getmetadata(io::IO, ::Type{T}, A::DictEncoding) where T
    vals = writecontents(Metadata.PrimitiveArray, io, levels(A))
    Metadata.CategoryMetadata(vals, true)
end

getmetadata(io::IO, ::Type{Union{Missing, T}}, A::DictEncoding) where T = getmetadata(io, T, A)


text:
module MimiSSPs

global g_ssp_datasets = Dict{Symbol,Any}()
global g_emissions_scenario_datasets = Dict{Symbol,Any}()

include("components/SSPs.jl")
include("components/RegionAggregatorSum.jl")

end # module


text:
function asracos_opt!(objective::Objective, parameter::Parameter)
    asracos = ASRacos(parameter.evaluation_server_num)
    rc = asracos.rc
    rc.objective = objective
    rc.parameter = parameter
    ub = isnull(parameter.uncertain_bits)? 1 : parameter.uncertain_bits

    # require calculator server
    ip = parameter.control_server_ip
    port = parameter.control_server_port
    cs_send = connect(ip, port)
    println(cs_send, "client: require servers#")
    readline(cs_send)
    println(cs_send, string(parameter.evaluation_server_num, "#"))
    msg = readline(cs_send)

    servers_msg = readline(cs_send)
    if servers_msg == "No available evaluation server"
        zoolog("Error: no available evaluation server")
        return Solution()
    end
    servers = split(servers_msg, " ")
    println("get $(length(servers)) servers")
    # close the socket
    close(cs_send)
    parameter.ip_port = RemoteChannel(()->Channel(length(servers)))
    for server in servers
        put!(parameter.ip_port, server)
    end
    asracos_init_attribute!(asracos, parameter)
    asracos_init_sample_set!(asracos, ub)
    println("Initialization succeeds")
    finish = SharedArray{Bool}(1)
    finish[1] = false
    # addprocs(1)
    @spawn asracos_updater!(asracos, parameter.budget, ub, finish)
    i = parameter.train_size
    br = false
    while true
        i += 1
        if finish[1] == true
            break
        end
        ip_port = take!(parameter.ip_port)
        sol = take!(asracos.sample_set)
        @spawn begin
            try
                br = compute_fx(sol, ip_port, parameter)
                put!(parameter.ip_port, ip_port)
                put!(asracos.result_set, sol)
            catch e
                println("Exception")
                println(e)
                cs_exception = connect(ip, port)
                println(cs_exception, "client: restart#")
                readline(cs_exception)
                println(cs_exception, string(servers_msg, "#"))
                close(cs_exception)
                return Solution()
            end
        end
    end
    # finish task
    result = take!(asracos.asyn_result)
    objective.history = take!(asracos.history)
    parameter.positive_data = take!(asracos.chan_positive_data)
    parameter.negative_data = take!(asracos.chan_negative_data)
    cs_receive = connect(ip, port)
    println(cs_receive, "client: return servers#")
    readline(cs_receive)
    println(cs_receive, string(servers_msg, "#"))
    close(cs_receive)
    return result
end

function compute_fx(sol::Solution, ip_port, parameter::Parameter)
    ip, port = get_ip_port(ip_port)
    client = connect(ip, port)

    # send calculate info
    println(client, "client: calculate#")
    msg = readline(client)
    br = false
    if check_exception(msg) == true
        br = true
    end

    println(client, "asracos#")
    msg = readline(client)
    if check_exception(msg) == true
        br = true
    end

    # send objective_file:func
    if br == false
        smsg = string(parameter.objective_file, ":", parameter.func, "#")
        println(client, smsg)
        msg = readline(client)
        if check_exception(msg) == true
            br = true
        end
    end

    # send x
    if br == false
        str = list2str(sol.x)
        println(client, str)
        receive = readline(client)
        if check_exception(receive) == true
            br = true
        end
    end
    if br == false
        value = parse(Float64, receive)
        sol.value = value
    end
    close(client)
    return br
end

function asracos_updater!(asracos::ASRacos, budget, ub, finish)
    rc = asracos.rc
    parameter = rc.parameter
    history = []
    t = parameter.train_size + 1
    strategy = parameter.replace_strategy
    time_log1 = now()
    interval = 10
    time_sum = interval
    output_file = parameter.output_file
    f = Nullable()
    if !isnull(output_file)
        f = open(output_file, "w")
    end
    br = false

    while(t <= budget)
        sol = take!(asracos.result_set)
        push!(history, sol.value)
        bad_ele = sracos_replace!(rc.positive_data, sol, "pos")
        sracos_replace!(rc.negative_data, bad_ele, "neg", strategy=strategy)
        rc.best_solution = rc.positive_data[1]
	    time_log2 = now()
        time_pass = Dates.value(time_log2 - time_log1) / 1000
        zoolog("Budget $(t): time=$(floor(time_pass))s, value=$(sol.value), best_solution_value=$(rc.best_solution.value)")
		str = "Budget $(t): time=$(floor(time_pass))s, value=$(sol.value), best_solution_value=$(rc.best_solution.value)\nbest_x=$(rc.best_solution.x)\n"
		if !isnull(f)
			write(f, str)
			flush(f)
		end
		if !isnull(parameter.time_limit) && time_pass > parameter.time_limit
			zoolog("Exceed time limit: $(parameter.time_limit)")
			br = true
		end
		 if rand(rng, Float64) < rc.parameter.probability
			 classifier = RacosClassification(rc.objective.dim, rc.positive_data,
			 rc.negative_data, ub=ub)
			 mixed_classification(classifier)
			 solution, distinct_flag = distinct_sample_classifier(rc, classifier, data_num=rc.parameter.train_size)
		 else
			 solution, distinct_flag = distinct_sample(rc, rc.objective.dim)
		 end
		 if distinct_flag == false
			 zoolog("ERROR: dimension limited")
			 break
		 end
		 if isnull(solution)
			 zoolog("ERROR: solution null")
			 break
		 end
		 put!(asracos.sample_set, solution)
		 t += 1
		 if br == true
			 break
		 end
     end
     finish[1] = true
     if !isnull(f)
         close(f)
     end
     put!(asracos.asyn_result, rc.best_solution)
     put!(asracos.history, history)
	 put!(asracos.chan_positive_data, rc.positive_data)
	 put!(asracos.chan_negative_data, rc.negative_data)
end

# function asracos_sample(rc, ub)
#     if rand(rng, Float64) < rc.parameter.probability
#         classifier = RacosClassification(rc.objective.dim, rc.positive_data,
#         rc.negative_data, ub=ub)
#         zoolog("before classification")
#         mixed_classification(classifier)
#         zoolog("after classification")
#         solution, distinct_flag = distinct_sample_classifier(rc, classifier, data_num=rc.parameter.train_size)
#     else
#         solution, distinct_flag = distinct_sample(rc, rc.objective.dim)
#     end
#     #painc stop
#     if distinct_flag == false
#         zoolog("ERROR: dimension limited")
#     end
#     if isnull(solution)
#         zoolog("ERROR: solution null")
#     end
#     return solution
# end

function asracos_init_attribute!(asracos::ASRacos, parameter::Parameter)
    f = open("init.txt", "w")
    rc = asracos.rc
    # otherwise generate random solutions
    iteration_num = rc.parameter.train_size
    data_temp = rc.parameter.init_sample
    init_num = 0
    remote_data = RemoteChannel(()->Channel(iteration_num))
    remote_result = RemoteChannel(()->Channel(iteration_num))
    if !isnull(data_temp)
        init_num = length(data_temp) < iteration_num? length(data_temp):iteration_num
        for i = 1:init_num
            # str = "initial sample: $(i), value=$(data_temp[i].value)\nx=$(sol.x)\n"
            put!(remote_data, data_temp[i])
        end
    end
    i = 1
    while i <= iteration_num - init_num
        # distinct_flag: True means sample is distinct(can be use),
        # False means sample is distinct, you should sample again.
        sol, distinct_flag = distinct_sample_from_set(rc, rc.objective.dim, rc.data,
            data_num=iteration_num)
        # panic stop
        if isnull(sol)
            break
        end
        if distinct_flag
            put!(remote_data, sol)
            i += 1
        end
    end
    fn = RemoteChannel(()->Channel(1))
    for i = 1:iteration_num
        d = take!(remote_data)
        if d.value != 0
            put!(remote_result, d)
            if i == iteration_num
                put!(fn, 1)
            end
            continue
        end
        ip_port = take!(parameter.ip_port)
        @spawn begin
            compute_fx(d, ip_port, parameter)
            put!(parameter.ip_port, ip_port)
            put!(remote_result, d)
            if i == iteration_num
                put!(fn, 1)
            end
        end

    end
    result = take!(fn)
    f = open("init.txt", "w")
    for i = 1:iteration_num
        d = take!(remote_result)
        push!(rc.data, d)
        str_print = "init sample: $(i), value=$(d.value)"
        str = "init sample: $(i), value=$(d.value)\n x=$(d.x)\n"
        zoolog(str_print)
        write(f, str)
        flush(f)
    end
    # zoolog("after taking result")
    close(f)
    selection!(rc)
    return
end

function construct_init_sample(init_file)
    f = open(init_file)
    lines = readlines(f)
    count = floor(Int, length(lines)/2)
    result = []
    for i in 1:count
        value_num = 2 * i - 1
        x_num = 2 * i
        value = eval(parse(split(lines[value_num], "=")[2]))
        x = eval(parse(split(lines[x_num], "=")[2]))
        sol = Solution(x=x, value=value)
        push!(result, sol)
    end
    return result
end

function check_exception(msg)
    if length(msg) < 9
        return false
    end
    res = msg[1:9]
    if res == "Exception"
        println(msg)
        return true
    end
    return false
end

function get_ip_port(ip_port)
    temp = split(ip_port, ":")
    ip = temp[1]
    port = parse(Int64, temp[2])
    return ip, port
end

function list2str(list)
    result = ""
    for i = 1:length(list)
        if i == 1
            result = string(list[i])
        else
            result = string(result, " ", string(list[i]))
        end
    end
    result = string(result, "#")
    result
end


text:
"""
    algorithm_cache(algorithm, problem, network)

Builds the cache for the given algorithm
"""
function algorithm_cache end

"""
    setup_algorithm!(alg_cache, ∇C, ∇logψ, par_cache)
"""
setup_algorithm!(g, data, grad, par_cache) = setup_algorithm!(g, data, par_cache)


text:
using Test

@testset "euler" begin
  include("wave_1d.jl")
  include("wave_2d.jl")
  include("wave_3d.jl")
  include("isentropicvortex.jl")
end


text:
using Revise
using Test, QuHamiltonian

@vertex i, j
@macroexpand @sum(Z[i]Z[j])

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
@testset "Test number of weights" begin
    for n = 1:5
        σ_parameters = MerweScaled(n=n, α=1.0, β=0.0, κ=3.0 - n)
        target_num = 2n + 1
        @test num_σ_points(σ_parameters) == target_num
        @test size(σ_parameters.m_w, 1) == target_num
        @test size(σ_parameters.Σ_w, 1) == target_num
    end
end

@testset "Test that weights are normalized" begin
    for n = 1:5
        σ_parameters = MerweScaled(n=n, α=1.0, β=0.0, κ=3.0 - n)
        @test sum(σ_parameters.m_w) - 1.0 < 1e-6
        @test sum(σ_parameters.Σ_w) - 1.0 < 1e-6
    end
end

max_stars_repo_name: pxl-th/KalmanFilters

max_issues_repo_name: pxl-th/KalmanFilters

max_forks_repo_name: pxl-th/KalmanFilters


text:
export resample, resample_by_FFT, resample_by_RFFT, upsample2_abs2, upsample2, upsample2_1D

"""
    resample(arr, new_size [, normalize])

Calculates the `sinc` interpolation of an `arr` on a new array size
`new_size`.
It is a re-evaluation of the Fourier series at new grid points.
`new_size` can be arbitrary. Each dimension is then independently either up or downsampled.

This method is based on FFTs and therefore implicitly assumes periodic
boundaries and a finite frequency support.

`normalize=true` by default multiplies by an appropriate factor so that
the array size is included in the scaling. This results in an array having roughly
the same mean intensity.

## Basic Principle
If `size(new_size)[i] > size(arr)[i]`, we apply zero padding in Fourier space.

If `size(new_size)[i] < size(arr)[i]`, we cut out a centered part of the
Fourier spectrum.

We apply some tricks at the boundary to increase accuracy of highest frequencies.

For real arrays we use `rfft` based operations, for complex one we use `fft` based ones.

# Examples

sinc interpolation of 2 datapoints result in an approximation of cosine.

```jldoctest
julia> resample([2.0, 0.0], (6,))
6-element Vector{Float64}:
 2.0
 1.5
 0.5
 0.0
 0.5
 1.5

julia> resample([2.0, 0.0], (6,)) ≈ 1 .+ cos.(2π .* (0:5)./6)
true
```
"""
function resample(arr::AbstractArray{T, N}, new_size; normalize=true) where {T, N}
    if new_size == size(arr)
        return copy(arr)
    end
    # for complex arrays we need a full FFT
    if T <: Complex
        arr_out = resample_by_FFT(arr, Tuple(new_size))
    else
        arr_out = resample_by_RFFT(arr, Tuple(new_size))
    end
    # normalize that values scale accordingly
    # this violates energy!
    if normalize
        arr_out .*= length(arr_out) ./ length(arr)
    end
    return arr_out
end

"""
    resample_by_RFFT(mat, new_size)

Does a resampling based on `rfft`. This function is called by `resampling`.
"""
function resample_by_RFFT(mat, new_size) where {T}
    old_size=size(mat)
    rf = rffts(mat)
    irffts(select_region_rft(rf,old_size,new_size), new_size[1])
end

"""
    resample_by_FFT(mat, new_size)

Does a resampling based on `fft`. This function is called by `resampling`.
"""
function resample_by_FFT(mat, new_size)
    old_size = size(mat)
    # for real arrays we apply an operation so that mat_fixed_before is hermitian
    mat_fixed_before = ft_fix_before(ffts(mat),old_size,new_size)
    mat_pad = ft_pad(mat_fixed_before,new_size)
    # afterwards we add the highest pos. frequency to the highest lowest one
    res_f = ft_fix_after(mat_pad, old_size,new_size)
    res = iffts(res_f)
    return res
end

function upsample2_1D(mat::AbstractArray{T, N}, dim=1, fix_center=false) where {T,N}
    newsize = Tuple((d==dim) ? 2*size(mat,d) : size(mat,d) for d in 1:N)
    res = zeros(eltype(mat), newsize)
    if fix_center && isodd(size(mat,dim))
        selectdim(res,dim,2:2:size(res,dim)) .= mat
        shifts = Tuple((d==dim) ? 0.5 : 0.0 for d in 1:N)
        selectdim(res,dim,1:2:size(res,dim)) .= shift(mat, shifts, take_real=true) # this is highly optimized and all fft of zero-shift directions are automatically avoided
    else
        selectdim(res,dim,1:2:size(res,dim)) .= mat
        shifts = Tuple((d==dim) ? -0.5 : 0.0 for d in 1:N)
        selectdim(res,dim,2:2:size(res,dim)) .= shift(mat, shifts, take_real=true) # this is highly optimized and all fft of zero-shift directions are automatically avoided
    end
    return res
end

"""
    upsample2(mat; dims=1:N)

Upsamples by a factor of two in all dimensions.
The code is optimized for speed by using subpixelshifts rather than Fourier resizing.
"""
function upsample2(mat::AbstractArray{T, N}; dims=1:N, fix_center=false) where {T,N}
    res = mat
    for d in dims
        res = upsample2_1D(res,d, fix_center)
    end
    return res
end

"""
    upsample2_abs2(mat::AbstractArray{T, N}; dims=1:N)

Upsamples by a factor of two and applies the abs2 operation. The code is optimized for speed.
"""
function upsample2_abs2(mat::AbstractArray{T, N}; dims=1:N) where {T,N}
    return abs2.(upsample2(mat, dims=dims))
end


text:
using NearestNeighbors
using StaticArrays

using Test
using LinearAlgebra

using Distances: Distances, Metric, evaluate, PeriodicEuclidean
struct CustomMetric1 <: Metric end
Distances.evaluate(::CustomMetric1, a::AbstractVector, b::AbstractVector) = maximum(abs.(a .- b))
function NearestNeighbors.interpolate(::CustomMetric1,
                                      a::V,
                                      b::V,
                                      x,
                                      d,
                                      ab) where {V <: AbstractVector}
    idx = (abs.(b .- a) .>= d - x)
    c = copy(Array(a))
    c[idx] = (1 - x / d) * a[idx] + (x / d) * b[idx]
    return c, true
end
struct CustomMetric2 <: Metric end
Distances.evaluate(::CustomMetric2, a::AbstractVector, b::AbstractVector) = norm(a - b) / (norm(a) + norm(b))

# TODO: Cityblock()
const metrics = [Chebyshev(), Euclidean(), Minkowski(3.5)]
const fullmetrics = [metrics; Hamming(); CustomMetric1(); CustomMetric2()]
const trees = [KDTree, BallTree]
const trees_with_brute = [BruteTree; trees]

include("test_knn.jl")
include("test_inrange.jl")
include("test_monkey.jl")
include("test_datafreetree.jl")

@testset "periodic euclidean" begin
    pred = PeriodicEuclidean([Inf, 2.5])
    l = [0.0 0.0; 0.0 2.5]
    S = BallTree(l, pred)
    @test inrange(S,[0.0,0.0], 1e-2, true) == [1, 2]
end


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

using Markdown
using InteractiveUtils

# ╔═╡ a0736b60-3c7d-11ec-1b71-27fd38cbf6e4
begin

# Import the necessary packages

using Distributions
using DelimitedFiles
using DataFrames
using Optim
using CSV
using Random
using Distributions
using Statistics
using LinearAlgebra
using ForwardDiff

end

# ╔═╡ 7b7526d0-3c9d-11ec-3575-37a66c520386
# Replication of Bruhin, Fehr, and Schunk, 2019, "Many Faces of Human Sociality: Uncovering the Distribution and Stability of Social Preferences"

# ╔═╡ 9d2d5b40-6023-11ec-0e3b-e5fe5c8d527b
# The code in this notebook performs the finite mixture estimates as in table 2

# ╔═╡ 7b36490e-3c9d-11ec-1ea6-dbf0082362dc
# Authors:
# Massimiliano Pozzi (Bocconi University, pozzi.massimiliano@studbocconi.it)
# Salvatore Nunnari (Bocconi University, https://snunnari.github.io, salvatore.nunnari@unibocconi.it)

# ╔═╡ 600cd6e0-4ea4-11ec-2071-13c721a5b90e
# This notebook works with both Neptune and Pluto; the "begin" and "end" tags in each cell are necessary for Pluto notebooks but do not play any role in Neptune

# ╔═╡ 61e5b590-4ea4-11ec-1bf8-a900f11ecb18
# This notebook was tested with the following packages versions:

# Pozzi: julia 1.6.1, Pluto 0.14.7, Neptune 0.14.0, Distributions 0.25.17, DataFrames 1.2.2, Optim 1.4.1, CSV 0.9.5, FiniteDiff 2.8.1, ForwardDiff 0.10.20

# ╔═╡ e0965590-3c7d-11ec-13b6-a7d15db3fd9f
begin

# 1. Data Cleaning and Preparation

# We import the relevant datasets containing the data on the 39 dictator games and 78 reciprocity games in Session 1 and Session 2. As the authors, we remove from these datasets those subjects that behaved very inconsistenly throughout the games. These subjects are identified from the individual estimates that we do not run in this notebook. The file "dropped_subjects_section4paragraph2.csv" contains the ID of these subjects.

# Load the data for Session 1

data, header = readdlm(raw"../input/choices_exp1.csv",',',header=true)
dt1 = DataFrame(data, vec(header))

# Import data with ID of subjects to drop

data_drop = readdlm(raw"../input/dropped_subjects_section4paragraph2.csv",',',header=false)

# Drop the session 1 subjects whose IDs are listed in the data_drop dataframe (14 individuals)

dt1 = dt1[[i for i=1:length(dt1.sid) if !(dt1.sid[i] in(data_drop))],:]

sort!(dt1,:sid) # individuals were not in order. We sort the dataframe so that the first 117 obs are for individual 1 and so on.

# Import session 2 data, remove subjects to drop and sort the dataframe as we did for session 1

data2, header2 = readdlm(raw"../input/choices_exp2.csv",',',header=true)
dt2 = DataFrame(data2, vec(header2))
dt2 = dt2[[i for i=1:length(dt2.sid) if !(dt2.sid[i] in(data_drop))],:]
sort!(dt2,:sid)

# import the aggregate estimates for session 1 and session 2 saved in the table1_julia.csv in the output folder. These will be used as a starting point for the EM-algorithm

res_s1=readdlm(raw"../output/table1_julia.csv",',',header=true)[1][:,2]
res_s1=convert(Array{Float64,1}, [res_s1[1:4];log(res_s1[5])])

res_s2=readdlm(raw"../output/table1_julia.csv",',',header=true)[1][:,6]
res_s2=convert(Array{Float64,1}, [res_s2[1:4];log(res_s2[5])])

end

# ╔═╡ d8928282-3c85-11ec-00ca-c5215c731b33
begin

# Create indicators_x and indicators_y. These are the columns s, r, q, v

indicators_x = hcat(dt1.s_x,dt1.r_x,dt1.q,dt1.v)
indicators_y = hcat(dt1.s_y,dt1.r_y,dt1.q,dt1.v)

indicators_x2 = hcat(dt2.s_x,dt2.r_x,dt2.q,dt2.v)
indicators_y2 = hcat(dt2.s_y,dt2.r_y,dt2.q,dt2.v)

end;

# ╔═╡ 47094410-4ea6-11ec-2865-47e35e532cf1
# For a more detailed explanation on the model we refer the reader to the python notebook.

# Below, we explain each step of the algoritm, but for a more rigorous discussion please refer to Dempster, A.P, N. M. Laird and D. B. Rubin, 1977, "Maximum Likelihood from Incomplete Data via the EM Algorithm", Journal of the Royal Statistical Society, Series B (Methodological), 39(1): 1–38 and to McLachlan, Geoffrey, Sharon Lee, and Suren Rathnayake, 2019, "Finite Mixture Models", Annual Review of Statistics and Its Application, 6: 355–378.

# Define the component density function. This would be f(θ_k, σ_k; X, Y, C_i) above for a specific set of parameters (θ_k, σ_k).
# v is the vector of parameters (θ_k, σ_k)
# y is the choice of the player
# self_x is the payoff when choosing x (left)
# self_y is the payoff when choosing y (right)
# indicators are the ones explained above
# This function returns a 1x160 vector with the likelihood contribution for each individual (nr. of unique individual is 160)

function singlecomp(v,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	beta  = v[1:(length(v)-1)] # vector of parameters alpha,beta,gamma,delta
    sigma = exp(v[length(v)])  # sigma. exp to force it positive during search

    lli = (indicators_x * beta)  # We get the matrix with rows [αs,βr,γq,δv]
    rli = (indicators_y * beta)

	uleft  = sigma.*( (1 .-lli) .*self_x .+lli .*other_x)
	uright = sigma.*( (1 .-rli) .*self_y .+rli .*other_y)

	probs  = (exp.(uleft)./(exp.(uleft).+exp.(uright))).^y .* (exp.(uright)./(exp.(uleft).+exp.(uright))).^(1 .-y)

	indprob = reshape(probs, 117, 160) # 117x160 matrix. each column is the likelihood of the observations for a single individual

	indlike = prod(indprob, dims=1) # 1x160 vector containing the likelihood contribution for a single individual. We took the product of a column

	# in theory it should not happen that a probability is equal to zero, but for some combinations of parameters it does. If prob = 0 we set it equal to a small value

	index_p0 = [i for i=1:length(indlike) if indlike[i]==0.0]

	for i in index_p0
        indlike[i] = 1E-8
	end

	return indlike

end

# ╔═╡ ca59c720-47a5-11ec-3cf6-f1b17e96c3e9
# This function returns a 160x3 matrix. Column one contains the individual likelihood contribution using the parameters for type one. Column two contains the individual likelihood contribution for type two etc. Using the notation from above a generic element i in the first column is what we called f(θ_1, σ_1; X, Y, C_i). Remember that we assumed there are three types (K=3).
# vall is a 1x15 vector that contains the parameters (θ,σ) for component 1,2 and 3, what we called (θ_1,σ_1,θ_2,σ_2,θ_3,σ_3).

function compdens(vall,y,self_x, other_x,self_y,other_y,indicators_x,indicators_y)

	theta1 = vall[1:5]    # parameters for first component  α1, β1, γ1, δ1, σ1
	theta2 = vall[6:10]   # parameters for second component α2, β2, γ2, δ2, σ2
	theta3 = vall[11:15]  # parameters for third component  α3, β3, γ3, δ3, σ3

	comp1 = singlecomp(theta1,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	comp2 = singlecomp(theta2,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	comp3 = singlecomp(theta3,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	comp = hcat(comp1',comp2',comp3') # 160x3 matrix explained above

	return comp

end

# ╔═╡ b77e0c80-4780-11ec-2a25-cbc0a9d7ca2f
# This function returns a 160x3 matrix containing the posterior probability of an individual to be of type one (first column), two (second column) and three (third column). Element in the first row first column contains the ex-post probability of individual 1 being of type 1, element in the first row second column contains the ex-post probability of individual 1 being of type 2 etc.
# mix is a 1x3 vector containing the shares of type one, two and three in the population, what we called (π_1, π_2, π_3)

function estep(vall,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y,mix)

	comp = compdens(vall,y,self_x, other_x,self_y,other_y,indicators_x,indicators_y)

# m is a 160x3 matrix. first colunm is the weight for the first component, second column the weight for the second component etc

	m = ones(160,3) .* mix'

	numerator = m .* comp # 160x3 matrix. each element is π_k f(θ_k,σ_k;X,Y,C_i)

	denominator = ones(160,3) .* (sum(numerator, dims=2)) # 160x3 matrix. elements in the same rows are equal. They are the denominator for individual 1 up to 160

# tau is our 160x3 matrix of ex-post probabilities for individual i belonging to type k

	tau = ones(160,3) .* (numerator ./ denominator)

	return tau

end

# ╔═╡ 930a4d20-477e-11ec-38a8-d57693e5c84e
# This function returns the negative of the sum of the individual observation log likelihood weighted by the ex-post probabilities for individuak i to belong to type k. This is the function to minimize in the M-step

function maxloglike(vall,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y, tau)

	comp = compdens(vall,y,self_x, other_x,self_y,other_y,indicators_x,indicators_y)

	mll = -sum(tau .* log.(comp))

	return mll

end

# ╔═╡ 93006210-477e-11ec-09d7-192bb2e62b61
# Define the M-step function. This function minimizes maxloglike and returns two vectors. mix is a 1x3 vector containing the shares of the types in the population, namely (π_1, π_2, π_3) that will be used in the E-step. vall_res is a 1x15 vector containing the new sets of parameters (θ_1,σ_1,θ_2,σ_2,θ_3,σ_3) that will be used in the E-step

function mstep(vall,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y, tau)

	ret = optimize(x->maxloglike(x,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y, tau),
			vall,method=NelderMead(),iterations = 10000)

	mix = mean(tau, dims=1)
	vall_res = Optim.minimizer(ret)

	return mix, vall_res

end

# ╔═╡ 0cf5a4f0-47c4-11ec-2f4f-b7efee6740c6
# This function returns the finite mixture negative log likelihood we want to minimize. This is the sum over individuals of the negative log of what we called L(Ψ;X,Y,C_i). psi is a 1x17 vector containing the 5 parameters for each type (5x3) and the two mixtures. It's what we defined as Ψ.

function totallike(psi,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	vall = psi[1:end-2]
	mix  = psi[end-1:end]
	mix  = [mix; 1-sum(mix)]

	comp = compdens(vall,y,self_x, other_x,self_y,other_y,indicators_x,indicators_y)
	m = ones(160,3) .* mix'

	tll = -sum(log.(sum(m .* comp, dims=2)))

	return tll

end

# ╔═╡ 92c13630-477e-11ec-13b8-3b6951af272c
# The function estimates runs the EM-algoritm and returns a series of objects.
# theta_v is the 1x15 vector containing the estimates on the parameters for each type (θ_1,σ_1, θ_2,σ_2, θ_3,σ_3)
# mix is a 1x3 vector containing the types' shares (π_1, π_2, π_3)
# the third object returned is the value of the log likelihood obtained
# tau is a 160x3 matrix of ex-post probabilities
# psi is a 1x17 vector containing theta_v and the first two mixtures

# Note that the prints show up in the console, not in the notebook

function estimate(y,self_x, other_x, self_y,other_y,indicators_x,indicators_y,maxiter,diffcrit)

# Starting conditions. random starting guesses close to the pooled ones. dimension is 1x15. 5 parameters for each type preferences (3 in our case)

	theta_v = repeat(res_s1,3) .* rand(Uniform(0.8,1.2), 15)

# Starting guesses for the mixtures

	m0 = rand(Uniform(0.5,2.0), 3)
	mix = m0/sum(m0)

# Initial EM control flow parameters

	iter  = 1 # nr iterations
	diff  = 1 # difference
	llold = 0 # old log like
	llnew = 1 # new log like

	while iter < maxiter+1 && abs(diff) > diffcrit

		println("EM type maximization: ", iter, "/", maxiter)
		println("\nWorking on E-step. \n")

		tau_em = estep(theta_v,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y,mix)

		println("Working on M-step. \n")

		m_step = mstep(theta_v,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y, tau_em)

		mix = m_step[1][:]
		theta_v = m_step[2]
		psi = [m_step[2];mix[1:2]]

		println(" sum of mixtures:", sum(mix), "\n")

# Adjust EM control flow parameters

		iter =  iter + 1
		llold = llnew
		llnew = totallike(psi,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)
		diff = llnew - llold

		println("Achieved Log Likelihood:", round(llnew,digits=4),"\n \n")
	end

	println("End of EM-steps \n")

# Maximize the likelihood directly.

	res = optimize(x->totallike(x,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y),
			psi,method=NelderMead(),iterations = 1000)

	println("Achieved final Log Likelihood:", round(Optim.minimum(res),digits=4),"\n")

	psi = Optim.minimizer(res)
	theta_v = psi[1:end-2]
	mix = psi[end-1:end]
	mix = [mix;1-sum(mix)]

# Retrieve the taus

	tau = estep(theta_v,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y,mix)

	return theta_v, mix, Optim.minimum(res), tau, psi
end

# ╔═╡ 1d4d9680-47cb-11ec-021f-5511d2dd1dc0
begin

# This cell computes the estimates. Rarely the last minimization of the objective function returns a domain error when too close to the minimum. If that happens just restart the cell.
# Note that prints will be shown in the terminal and not in the notebook

solmixture = estimate(dt1.choice_x,dt1.self_x,dt1.other_x,dt1.self_y,dt1.other_y,indicators_x,indicators_y,15,5e-5)

end;

# ╔═╡ 77c1e520-4c3b-11ec-3db8-a58c85ca0e8f
begin

# the results we want

theta_v = solmixture[1]           # parameters for the three types
mixture = solmixture[2]           # the mixture probabilities
likelihood = solmixture[3]        # value of the likelihood
posterior_prob = solmixture[4]    # individual posterior prob. of being of a type
psi = solmixture[5]
results_sm1 = round.([psi[1:4];exp(psi[5]);psi[6:9];exp(psi[10]);psi[11:14];exp(psi[15]);psi[16:17];1-sum(psi[16:17])],digits=3)

end;

# ╔═╡ 2c418020-4c44-11ec-24f8-13b7ff1ffcdf
begin

# For a more detailed explanation on how to compute cluster robust standard errors we refer the reader to the python notebook or to David A. Freedman, 2006, "On The So-Called 'Huber Sandwich Estimator' and 'Robust Standard Errors'", The American Statistician, 60:4, 299-302).

# Compute the individual cluster robust standard errors.
# The formula is the following:

# adj * (inv_hessian * grad_con * inv_hessian).

# adj is a correction for the degree of freedoms and the number of clusters. inv_hessian is the inverse of the hessian of the negative log-likelihood evaluated at the minimum. grad_con is is a 5x5 matrix of gradient contributions. It is computed as the sum over all clusters of the outer product of the sum of the single observation gradient evaluated in the minimum for all observations of an individual cluster.

# Function to compute mixture negative log likelihood. Same as totallike but here mix is the full vector, so that the hessian is a 18x18 matrix

function totallike_mix(psi,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	# psi contains the 15 parameters, 5 for each type and the three mixture prob.

	vall = psi[1:15]   # 1x15 vector containing parameters for each type
	mix  = psi[16:end] # three mixture prob

	comp = compdens(vall,y,self_x, other_x,self_y,other_y,indicators_x,indicators_y)
	m = ones(160,3) .* mix'

	tll = -sum(log.(sum(m .* comp, dims=2)))

	return tll

end

# compute the inverse of the hessian

hessian_ms1 = ForwardDiff.hessian(x->totallike_mix(x,dt1.choice_x,dt1.self_x,dt1.other_x,dt1.self_y,dt1.other_y,indicators_x,indicators_y),[theta_v;mixture])

inv_hess_ms1 = inv(hessian_ms1)

end;

# ╔═╡ f2b8a210-4df7-11ec-1676-6dd04ccc9a65
begin

# Define an auxiliary function.
# This function takes as input a dataset and a scalar and returns a dataframe containing only observations for individual whose ID is equal to that scalar

function getind(dataset,sid)
    dataset_ind = dataset[dataset.sid .== sid,:]
    return dataset_ind
end

function ind_singlecomp(v,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	beta  = v[1:(length(v)-1)] # vector of parameters alpha,beta,gamma,delta
    sigma = exp(v[length(v)])  # sigma. exp to force it positive during search

    lli = (indicators_x * beta)  # We get the matrix with rows [αs,βr,γq,δv]
    rli = (indicators_y * beta)

	uleft  = sigma .* ( (1 .- lli) .* self_x .+ lli .* other_x)
	uright = sigma .* ( (1 .- rli) .* self_y .+ rli .* other_y)

	probs  = (exp.(uleft)./(exp.(uleft).+exp.(uright))).^y .* (exp.(uright)./(exp.(uleft).+exp.(uright))).^(1 .-y)

	indlike = prod(probs)

	if indlike==0.0 indlike=1e-8 end

	return indlike

end

function ind_compdens(vall,y,self_x, other_x,self_y,other_y,indicators_x,indicators_y)

	# vall is a 1x15 vector that contains the parameters θ for component 1,2 and 3

	theta1 = vall[1:5]    # parameters for first component  α1, β1, γ1, δ1, σ1
	theta2 = vall[6:10]   # parameters for second component α2, β2, γ2, δ2, σ2
	theta3 = vall[11:15]  # parameters for third component  α3, β3, γ3, δ3, σ3

	comp1 = ind_singlecomp(theta1,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	comp2 = ind_singlecomp(theta2,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	comp3 = ind_singlecomp(theta3,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

# comp is a 1x3 vector. The first column is component1 (where each element in column is the likelihood of individual i when we use as parameters the ones for the first component). Same for column 2 and 3.

	comp = [comp1,comp2,comp3]

end

function ind_totallike(psi,y,self_x, other_x, self_y,other_y,indicators_x,indicators_y)

	# psi contains the 15 parameters, 5 for each type and the two mixture prob.

	vall = psi[1:15] # 1x15 vector containing parameters for each type
	mix  = psi[16:end]
	comp = ind_compdens(vall,y,self_x, other_x,self_y,other_y,indicators_x,indicators_y)

	tll = -log.(sum(mix .* comp))

	return tll

end

end

# ╔═╡ 985c0c70-4c40-11ec-0c1a-a395646d2ef7
# Define the function that computes the matrix of individual gradient contribution

function grad_contr(dataset, parameters)
	G_sm1 = zeros(18,18)
	for sid in unique(dataset.sid)
		dt1_ind = getind(dataset,sid)
	 ind_grad = ForwardDiff.gradient(x->ind_totallike(x,dt1_ind.choice_x,dt1_ind.self_x,dt1_ind.other_x,dt1_ind.self_y,dt1_ind.other_y,hcat(dt1_ind.s_x,dt1_ind.r_x,dt1_ind.q,dt1_ind.v),hcat(dt1_ind.s_y,dt1_ind.r_y,dt1_ind.q,dt1_ind.v)), parameters)
	G_sm1 .+= (ind_grad*ind_grad')
	end
	return G_sm1
end

# ╔═╡ aeef235e-4c5d-11ec-0a38-fbe40474af45
begin

# Do the same for session 2

solmixture2 = estimate(dt2.choice_x,dt2.self_x,dt2.other_x,dt2.self_y,dt2.other_y,indicators_x2,indicators_y2,15,5e-3);

# the results we want

theta_v2 = solmixture2[1]           # parameters for the three types
mixture2 = solmixture2[2]           # the mixture probabilities
likelihood2 = solmixture2[3]        # value of the likelihood
posterior_prob2 = solmixture2[4]    # individual posterior prob. of being of a type
psi2 = solmixture2[5]
results_sm2 = round.([psi2[1:4];exp(psi2[5]);psi2[6:9];exp(psi2[10]);psi2[11:14];exp(psi2[15]);psi2[16:17];1-sum(psi2[16:17])],digits=3)

end;

# ╔═╡ e85bb2fe-4c6e-11ec-3e51-3929725157f8
begin

# compute the inverse of the hessian

hessian_ms2 = ForwardDiff.hessian(x->totallike_mix(x,dt2.choice_x,dt2.self_x,dt2.other_x,dt2.self_y,dt2.other_y,indicators_x2,indicators_y2),[theta_v2;mixture2])

inv_hess_ms2 = inv(hessian_ms2)

# compute gradient contribution G

grad_contr_ms2 = grad_contr(dt2, [theta_v2;mixture2])

# degree of freedoms and number of clusters adjustment

adj2 = (length(dt2.self_x)-1) / (length(dt2.self_x)-length([theta_v;mixture])) * (160/159)

# Sandwich formula for the clustered standard errors

varcov_ms2 = adj2 .* (inv_hess_ms2 * grad_contr_ms2 * inv_hess_ms2)
se_sm2 = sqrt.(diag(varcov_ms2))
se_sm2 = [se_sm2[1:4]; sqrt(exp(psi2[5])^2*se_sm2[5]^2);se_sm2[6:9]; sqrt(exp(psi2[10])^2*se_sm2[10]^2);se_sm2[11:14]; sqrt(exp(psi2[15])^2*se_sm2[15]^2);se_sm2[16:18]]

# sqrt(exp(res_s1[5])^2*se_s1[5]^2) is the formula to go from standard error exp(sigma) to standard error of sigma

se_sm2 = round.(se_sm2, digits=3)

end;

# ╔═╡ 5de4a510-4eb7-11ec-0a34-41d9b6a79e59
begin

# Hypothesis Testing

# We now do some hypothesis testing on the parameters we obtained in session 1 and session 2. We compute the z-test statistics and the corresponding p-values to check if each parameter we obtained is statistically different from zero.

# Compute the z-test statistics and the corresponding p-values to check if the parameters are statistically different from zero

# Session 1

zvalues_s1 = abs.(results_sm1 ./ se_sm1)
pvalues_s1 = 2*(1 .-cdf.(Normal(0,1),zvalues_s1))

# Session 2

zvalues_s2 = abs.(results_sm2 ./ se_sm2)
pvalues_s2 = 2*(1 .-cdf.(Normal(0,1),zvalues_s2))

end;

# ╔═╡ b5b26ef0-4dfc-11ec-3539-b70db077382f
begin

# print the results for session 1. Note that the parameters will not be necessarily in this order: strongly altruistic, moderately altruistic, behindness averse. It depends on the starting conditions of the mixtures and on how the algorithm moves.

# There are some discrepancies between the standard errors we obtain and those in Table 2 of the paper. The differences in the SEs for the point estimates of Session 1 are due to small differences in the point estimates themselves which are likely due to the use of a different minimization algorithm (we use Nelder-Mead, which is gradient-free, while the authors write down derivatives in closed form and use BFGS). The differences in the SEs for the mixture probabilities of Sessions 1 and 2 are probably due to the fact that our code takes a shortcut and computes the gradient of a subject's log likelihood instead of computing the sum of the gradients of each single choice's log likelihood for the same subject (which is the theoretically correct way of computing the SEs; see, for example, the discussion in the paper cited above). The reason for taking this shortcut is that the theoretically correct object is complex to compute, given the nature of the log likelihood function, and that the shortcut actually returns the same SEs as the authors' for the other parameters.

param_names = ["π: Type's shares in the population",
			   "α: Weight on other's payoff when behind",
	           "β: Weight on other's payoff when ahead",
	           "γ: Measure of positive reciprocity",
	           "δ: Measure of negative reciprocity",
			   "σ: Choice sensitivity"]

table21 = DataFrame()
table21.param_name = param_names
table21.str_altr_est = [results_sm1[end];results_sm1[11:15]]
table21.str_altr_se = [se_sm1[end];se_sm1[11:15]]
table21.pval_str = round.([pvalues_s1[end];pvalues_s1[11:15]],digits=3)
table21.mod_altr_est = [results_sm1[end-2];results_sm1[1:5]]
table21.mod_altr_se = [se_sm1[end-2];se_sm1[1:5]]
table21.pval_mod = round.([pvalues_s1[end-2];pvalues_s1[1:5]],digits=3)
table21.beh_aver_est = [results_sm1[end-1];results_sm1[6:10]]
table21.beh_aver_se = [se_sm1[end-1];se_sm1[6:10]]
table21.pval_beh = round.([pvalues_s1[end-1];pvalues_s1[6:10]],digits=3)
table21

end

# ╔═╡ b7f90390-4dfc-11ec-2c4d-b3dedc1ed52d
begin

# print the results for session 2

table22 = DataFrame()
table22.param_name = param_names
table22.str_altr_est = [results_sm2[end];results_sm2[11:15]]
table22.str_altr_se = [se_sm2[end];se_sm2[11:15]]
table22.pval_str = round.([pvalues_s2[end];pvalues_s2[11:15]],digits=3)
table22.mod_altr_est = [results_sm2[end-2];results_sm2[1:5]]
table22.mod_altr_se = [se_sm2[end-2];se_sm2[1:5]]
table22.pval_mod = round.([pvalues_s2[end-2];pvalues_s2[1:5]],digits=3)
table22.beh_aver_est = [results_sm2[end-1];results_sm2[6:10]]
table22.beh_aver_se = [se_sm2[end-1];se_sm2[6:10]]
table22.pval_beh = round.([pvalues_s2[end-1];pvalues_s2[6:10]],digits=3)
table22

end

# ╔═╡ 06309850-4dff-11ec-04ba-e3b11d017c68
begin

# save the results as a csv file

table2 = DataFrame()
table2.param_name = param_names
table2.str_altr_est_s1 = [results_sm1[end];results_sm1[11:15]]
table2.str_altr_se_s1 = [se_sm1[end];se_sm1[11:15]]
table2.pval_str_s1 = round.([pvalues_s1[end];pvalues_s1[11:15]],digits=3)
table2.mod_altr_est_s1 = [results_sm1[end-2];results_sm1[1:5]]
table2.mod_altr_se_s1 = [se_sm1[end-2];se_sm1[1:5]]
table2.pval_mod_s1 = round.([pvalues_s1[end-2];pvalues_s1[1:5]],digits=3)
table2.beh_aver_est_s1 = [results_sm1[end-1];results_sm1[6:10]]
table2.beh_aver_se_s1 = [se_sm1[end-1];se_sm1[6:10]]
table2.pval_beh_s1 = round.([pvalues_s1[end-1];pvalues_s1[6:10]],digits=3)
table2.str_altr_est_s2 = [results_sm2[end-1];results_sm2[6:10]]
table2.str_altr_se_s2 = [se_sm2[end-1];se_sm2[6:10]]
table2.pval_str_s2 = round.([pvalues_s2[end-1];pvalues_s2[6:10]],digits=3)
table2.mod_altr_est_s2 = [results_sm2[end];results_sm2[11:15]]
table2.mod_altr_se_s2 = [se_sm2[end];se_sm2[11:15]]
table2.pval_mod_s2 = round.([pvalues_s2[end];pvalues_s2[11:15]],digits=3)
table2.beh_aver_est_s2 = [results_sm2[end-2];results_sm2[1:5]]
table2.beh_aver_se_s2 = [se_sm2[end-2];se_sm2[1:5]]
table2.pval_beh_s2 = round.([pvalues_s2[end-2];pvalues_s2[1:5]],digits=3)

CSV.write(raw"../output/table2_julia.csv", table2)

end

# ╔═╡ Cell order:
# ╠═7b7526d0-3c9d-11ec-3575-37a66c520386
# ╠═9d2d5b40-6023-11ec-0e3b-e5fe5c8d527b
# ╠═7b36490e-3c9d-11ec-1ea6-dbf0082362dc
# ╠═600cd6e0-4ea4-11ec-2071-13c721a5b90e
# ╠═61e5b590-4ea4-11ec-1bf8-a900f11ecb18
# ╠═a0736b60-3c7d-11ec-1b71-27fd38cbf6e4
# ╠═e0965590-3c7d-11ec-13b6-a7d15db3fd9f
# ╠═d8928282-3c85-11ec-00ca-c5215c731b33
# ╠═47094410-4ea6-11ec-2865-47e35e532cf1
# ╠═ca59c720-47a5-11ec-3cf6-f1b17e96c3e9
# ╠═b77e0c80-4780-11ec-2a25-cbc0a9d7ca2f
# ╠═930a4d20-477e-11ec-38a8-d57693e5c84e
# ╠═93006210-477e-11ec-09d7-192bb2e62b61
# ╠═0cf5a4f0-47c4-11ec-2f4f-b7efee6740c6
# ╠═92c13630-477e-11ec-13b8-3b6951af272c
# ╠═1d4d9680-47cb-11ec-021f-5511d2dd1dc0
# ╠═77c1e520-4c3b-11ec-3db8-a58c85ca0e8f
# ╠═2c418020-4c44-11ec-24f8-13b7ff1ffcdf
# ╠═f2b8a210-4df7-11ec-1676-6dd04ccc9a65
# ╠═985c0c70-4c40-11ec-0c1a-a395646d2ef7
# ╠═aeef235e-4c5d-11ec-0a38-fbe40474af45
# ╠═e85bb2fe-4c6e-11ec-3e51-3929725157f8
# ╠═5de4a510-4eb7-11ec-0a34-41d9b6a79e59
# ╠═b5b26ef0-4dfc-11ec-3539-b70db077382f
# ╠═b7f90390-4dfc-11ec-2c4d-b3dedc1ed52d
# ╠═06309850-4dff-11ec-04ba-e3b11d017c68


text:
include("../src/utils/readfiles.jl")
include("../src/setting/constants.jl")
using ProgressMeter
using Base.Threads

"""
generate perturbation file according to the reference bin file
"""
function generate_perturbation(target_basedir::String, reference_basedir::String, output_basedir::String, tags::String, nproc::Int64, nspec::Int64)
    # p = Progress(nproc)
    @threads for iproc in 0:nproc - 1
        kernel_generate_perturbation(target_basedir, reference_basedir, output_basedir, tags, iproc, nspec)
        # next!(p)
    end
end

function kernel_generate_perturbation(target_basedir::String, reference_basedir::String, output_basedir::String, tags::String, iproc::Int64, nspec::Int64)
    tags_splitted = split(tags, ",")
    model_gll_reference = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_target = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_output = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    for tag in tags_splitted
        # convert tag to String
        tag = String(tag)
        sem_io_read_gll_file_1!(reference_basedir, iproc, tag, model_gll_reference)
        sem_io_read_gll_file_1!(target_basedir, iproc, tag, model_gll_target)

        model_gll_output = (model_gll_target .- model_gll_reference) ./ model_gll_reference
        sem_io_write_gll_file_1(output_basedir, iproc, tag, model_gll_output)
    end
end

"""
generate the real vaue bin file based on reference model and perturbation bin file
"""
function generate_real(target_basedir::String, reference_basedir::String, output_basedir::String, tags::String, nproc::Int64, nspec::Int64)
    tags_splitted = split(tags, ",")
    p = Progress(nproc)
    @threads for iproc in 0:nproc - 1
        model_gll_reference = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
        model_gll_target = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
        model_gll_output = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
        for tag in tags_splitted
            # convert tag to String
            tag = String(tag)
            sem_io_read_gll_file_1!(reference_basedir, iproc, tag, model_gll_reference)
            sem_io_read_gll_file_1!(target_basedir, iproc, tag, model_gll_target)

            model_gll_output = model_gll_reference .* (1 .+ model_gll_target)
            sem_io_write_gll_file_1(output_basedir, iproc, tag, model_gll_output)
        end
        next!(p)
    end
end

"""
calculate the difference between two sets of gll files.
"""
function generate_difference(target_basedir::String, reference_basedir::String, output_basedir::String, tags::String, nproc::Int64, nspec::Int64)
    tags_splitted = split(tags, ",")
    p = Progress(nproc)
    @threads for iproc in 0:nproc - 1
        model_gll_reference = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
        model_gll_target = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
        model_gll_output = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
        for tag in tags_splitted
            # convert tag to String
            tag = String(tag)
            sem_io_read_gll_file_1!(reference_basedir, iproc, tag, model_gll_reference)
            sem_io_read_gll_file_1!(target_basedir, iproc, tag, model_gll_target)

            model_gll_output = model_gll_target .- model_gll_reference
            sem_io_write_gll_file_1(output_basedir, iproc, tag, model_gll_output)
        end
        next!(p)
    end
end

"""
calculate the model with the minus sign
"""
function minus_sign(target_basedir::String, output_basedir::String, tags::String, nproc::Int64, nspec::Int64)
    tags_splitted = split(tags, ",")
    p = Progress(nproc)
    @threads for iproc in 0:nproc - 1
        model_gll_target = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
        model_gll_output = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
        for tag in tags_splitted
            # convert tag to String
            tag = String(tag)
            sem_io_read_gll_file_1!(target_basedir, iproc, tag, model_gll_target)

            model_gll_output = -model_gll_target
            sem_io_write_gll_file_1(output_basedir, iproc, tag, model_gll_output)
        end
        next!(p)
    end
end

"""
get new model from the perturbation.
"""
function generate_new(old_basedir::String, per_basedir::String, output_basedir::String, tags::String, nproc::Int64, nspec::Int64)
    p = Progress(nproc)
    @threads for iproc in 0:nproc - 1
        kernel_generate_new(old_basedir, per_basedir, output_basedir, tags, iproc, nspec)
        next!(p)
    end
end

function kernel_generate_new(old_basedir::String, per_basedir::String, output_basedir::String, tags::String, iproc::Int64, nspec::Int64)
    tags_splitted = split(tags, ",")
    model_gll_per = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_old = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_output = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    for tag in tags_splitted
        try
            # convert tag to String
            tag = String(tag)
            tag_old = tag * "_new"
            sem_io_read_gll_file_1!(old_basedir, iproc, tag_old, model_gll_old)
            # modify tag here
            tag_per = "d" * tag * tag
            sem_io_read_gll_file_1!(per_basedir, iproc, tag_per, model_gll_per)

            model_gll_output = model_gll_old .* (exp.(model_gll_per))
            sem_io_write_gll_file_1(output_basedir, iproc, tag_old, model_gll_output)
        catch e
        end
    end
end

"""
Compare two models.
"""
function compare_model(target_basedir::String, reference_basedir::String, tag::String,  nspec::Int64, iproc::Int64)
    model_gll_target = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_reference = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    sem_io_read_gll_file_1!(target_basedir, iproc, tag, model_gll_target)
    sem_io_read_gll_file_1!(reference_basedir, iproc, tag, model_gll_reference)
    model_diff = model_gll_target .- model_gll_reference
    @info "max difference" maximum(model_diff)
end

"""
convert vpv,vph,vsv,vsh to bulk_cv,bulk_ch,vsv,vsh.
"""
function vpvs2bulk_c!(target_basedir::String, output_basedir::String, nspec::Int64, iproc::Int64)
    model_gll_vpv = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_vph = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_vsv = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_vsh = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    # * read in models
    sem_io_read_gll_file_1!(target_basedir, iproc, "vpv", model_gll_vpv)
    sem_io_read_gll_file_1!(target_basedir, iproc, "vph", model_gll_vph)
    sem_io_read_gll_file_1!(target_basedir, iproc, "vsv", model_gll_vsv)
    sem_io_read_gll_file_1!(target_basedir, iproc, "vsh", model_gll_vsh)
    # * now we override vpv,vph with the bulk_c, calculate based on:
    # * ϕ^2=α^2-4/3*β^2
    model_gll_vpv = @. sqrt(model_gll_vpv^2 - 4 / 3 * model_gll_vsv^2)
    model_gll_vph = @. sqrt(model_gll_vph^2 - 4 / 3 * model_gll_vsh^2)
    # * write out files
    sem_io_write_gll_file_1(output_basedir, iproc, "bulk_cv", model_gll_vpv)
    sem_io_write_gll_file_1(output_basedir, iproc, "bulk_ch", model_gll_vph)
    sem_io_write_gll_file_1(output_basedir, iproc, "vsv", model_gll_vsv)
    sem_io_write_gll_file_1(output_basedir, iproc, "vsh", model_gll_vsh)
end

"""
convert bulk_cv,bulk_ch,vsv,vsh to vpv,vph,vsv,vsh.
"""
function bulk_c2vpvs!(target_basedir::String, output_basedir::String, nspec::Int64, iproc::Int64)
    model_gll_vpv = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_vph = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_vsv = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    model_gll_vsh = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    # * read in models
    sem_io_read_gll_file_1!(target_basedir, iproc, "bulk_cv", model_gll_vpv)
    sem_io_read_gll_file_1!(target_basedir, iproc, "bulk_ch", model_gll_vph)
    sem_io_read_gll_file_1!(target_basedir, iproc, "vsv", model_gll_vsv)
    sem_io_read_gll_file_1!(target_basedir, iproc, "vsh", model_gll_vsh)
    # * now we override bilk_c with vpv and vph, calculate based on:
    # * ϕ^2=α^2-4/3*β^2 so α^2=ϕ^2+4/3*β^2
    model_gll_vpv = @. sqrt(model_gll_vpv^2 + 4 / 3 * model_gll_vsv^2)
    model_gll_vph = @. sqrt(model_gll_vph^2 + 4 / 3 * model_gll_vsh^2)
    # * write out files
    sem_io_write_gll_file_1(output_basedir, iproc, "vpv", model_gll_vpv)
    sem_io_write_gll_file_1(output_basedir, iproc, "vph", model_gll_vph)
    sem_io_write_gll_file_1(output_basedir, iproc, "vsv", model_gll_vsv)
    sem_io_write_gll_file_1(output_basedir, iproc, "vsh", model_gll_vsh)
end

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
using Genie.Configuration, Logging

const config = Settings(
  server_port                     = 8000,
  server_host                     = "0.0.0.0",
  log_level                       = Logging.Error,
  log_to_file                     = false,
  server_handle_static_files      = true
)

if config.server_handle_static_files
  @warn("For performance reasons Genie should not serve static files (.css, .js, .jpg, .png, etc) in production.
         It is recommended to set up Apache or Nginx as a reverse proxy and cache to serve static assets.")
end

ENV["JULIA_REVISE"] = "off"

max_stars_repo_name: essenciary/performance-tests

max_issues_repo_name: essenciary/performance-tests

max_forks_repo_name: essenciary/performance-tests


text:
#! format: off
module QuantumControl
include("reexport.jl")

using QuantumPropagators
@reexport_members(QuantumPropagators)

using QuantumControlBase
@reexport_members(QuantumControlBase)

module Shapes
    # we need `QuantumControlBase.Shapes` to be available under a name that
    # doesn't clash with `QuantumControl.Shapes` in order for the
    # `@reexport_members` macro to work correctly
    using QuantumControlBase: Shapes as QuantumControlBase_Shapes
    using QuantumControlBase.Shapes
    include("reexport.jl")
    @reexport_members(QuantumControlBase_Shapes)
end

module Functionals
    using QuantumControlBase: Functionals as QuantumControlBase_Functionals
    using QuantumControlBase.Functionals
    include("reexport.jl")
    @reexport_members(QuantumControlBase_Functionals)
end

module WeylChamber
    using QuantumControlBase: WeylChamber as QuantumControlBase_WeylChamber
    using QuantumControlBase.WeylChamber
    include("reexport.jl")
    @reexport_members(QuantumControlBase_WeylChamber)
end

using Krotov
using GRAPE

include("print_versions.jl")

end


text:
# Data topmute
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: December 2017

# Updated by Ziyi Yin, November 2020, ziyi.yin@gatech.edu
#
export marineTopmute2D, judiMarineTopmute2D
export model_topmute, judiTopmute, find_water_bottom, depth_scaling, judiDepthScaling, laplace, low_filter, judiFilter

############################################ Data space preconditioners ################################################

function judiFilter(geometry, fmin, fmax)
    nsrc = length(geometry.xloc)
    N = 0
    for j=1:nsrc
        N += geometry.nt[j]*length(geometry.xloc[j])
    end
    D = joLinearFunctionFwd_T(N, N,
                              v -> low_filter(v,geometry.dt[1];fmin=fmin, fmax=fmax),
                              w -> low_filter(w,geometry.dt[1];fmin=fmin, fmax=fmax),
                              Float32,Float32,name="Data filter")
    return D
end

function low_filter(Din::Array{Float32, 1}, dt_in; fmin=0.0, fmax=25.0)
    Dout = deepcopy(Din)
    responsetype = Bandpass(fmin, fmax; fs=1e3/dt_in)
    designmethod = Butterworth(5)
    return filt(digitalfilter(responsetype, designmethod), Float32.(Dout))
end

function low_filter(Din::Array{Float32, 2}, dt_in; fmin=0.0, fmax=25.0)
    Dout = deepcopy(Din)
    responsetype = Bandpass(fmin, fmax; fs=1e3/dt_in)
    designmethod = Butterworth(5)
    for i=1:size(Dout,2)
        Dout[:, i] = filt(digitalfilter(responsetype, designmethod), Float32.(Dout[:, i]))
    end
    return Dout
end

function low_filter(Din::judiVector, dt_in; fmin=0.0, fmax=25.0)
    Dout = deepcopy(Din)
    for j=1:Dout.nsrc
		if size(Din.data[j], 2) == 1
			Dout.data[j] = low_filter(Dout.data[j], dt_in; fmin=fmin, fmax=fmax)
		else
	        for i=1:size(Din.data[j], 2)
	            Dout.data[j][:, i] = low_filter(Dout.data[j][:, i], dt_in; fmin=fmin, fmax=fmax)
	        end
		end
    end
    return Dout
end

function marineTopmute2D(Dobs::judiVector, muteStart::Integer; mute=Array{Any}(undef, 3), flipmask=false)
    # Data topmute for end-on spread marine streamer data
    Din = deepcopy(Dobs)

    for j=1:Din.nsrc

        # Design mute window
        x0 = 1f0
        xend = length(Din[j].geometry.xloc[1])
        nt = Din[j].geometry.nt[1]
        nrec = length(Din[j].geometry.xloc[1])
        drec = abs(Din[j].geometry.xloc[1][1] - Din[j].geometry.xloc[1][2])
        offsetDirectWave = 1.5f0*Din[j].geometry.t[1]
        idxOffset = Int(round(offsetDirectWave/drec))
        dx = round(idxOffset - idxOffset/10f0)

        if j==1 && ~isassigned(mute)
            z0 = muteStart - Int(round(muteStart/10))
            slope = 1.05f0*(nt - z0)/dx
            mute[1] = x0
            mute[2] = z0
            mute[3] = slope
        else#if j==1 && isassigned(mute)
            x0 = Int64(mute[1])
            z0 = Int64(mute[2])
            slope = mute[3]
        end

        mask = ones(Float32,nt,nrec)
        mask[1:z0,:] .= 0f0

        # Linear mute
        if (nrec-x0 < dx)
            x = nrec
            zIntercept = Int(round(z0+slope*(x-x0)))
            zax = z0+1:1:zIntercept
        else
            x = x0+dx
            zax = z0+1:1:nt
        end
        if length(zax) > 1
            xax = Array{Int}(round.(range(x0,stop=x,length=length(zax))))
        else
            xax = Int(round(x0))
        end
        for k=1:length(zax)
            mask[min(zax[k], nt),xax[k]:end] .= 0f0
        end
        flipmask == true && (mask = reverse(mask, dims=2))
        Din.data[j] = Din.data[j].*mask
    end
    return Din
end

function judiMarineTopmute2D(muteStart,geometry;params=Array{Any}(undef, 3),flipmask=false)
# JOLI wrapper for the linear depth scaling function
    nsrc = length(geometry.xloc)
    N = 0
    for j=1:nsrc
        N += geometry.nt[j]*length(geometry.xloc[j])
    end
    D = joLinearFunctionFwd_T(N,N,
                             v -> marineTopmute2D(v,muteStart;mute=params, flipmask=flipmask),
                             w -> marineTopmute2D(w,muteStart;mute=params, flipmask=flipmask),
                             Float32,Float32,name="Data mute")
    return D
end

############################################ Model space preconditioners ###############################################

function model_topmute(n::Tuple{Int64,Int64}, mute_end::Array{Integer,1}, length::Int64, x_orig)
# Model domain topmute for a velocity model of dimensions n = [nx, nz].
    x = deepcopy(reshape(x_orig,n))
    for j=1:n[1]
        mute_start = mute_end[j] - length
        filter = zeros(Float32, n[2])
        filter[1:mute_start-1] .= 0f0
        filter[mute_end[j]+1:end] .= 1f0
        taper_length = mute_end[j] - mute_start + 1
        taper = (1f0 .+ sin.((Float32(pi)*(0:taper_length-1))/(taper_length - 1) .- Float32(pi)/2f0))/2f0
        filter[mute_start:mute_end[j]] = taper
        global x[j,:] = x[j,:].*filter
    end
    return vec(x)
end

function model_topmute(n::Tuple{Int64,Int64}, mute_end::Int64, length::Int64, x_orig)
# Model domain topmute for a velocity model of dimensions n = [nx, nz].
    x = deepcopy(reshape(x_orig,n))
    mute_start = mute_end - length
    filter = zeros(Float32, n[2])
    filter[1:mute_start-1] .= 0f0
    filter[mute_end+1:end] .= 1f0
    taper_length = mute_end - mute_start + 1
    taper = (1f0 .+ sin.((Float32(pi)*(0:taper_length-1))/(taper_length - 1) .- Float32(pi)/2f0))/2f0
    filter[mute_start:mute_end] = taper
    for j=1:n[1]
        global x[j,:] = x[j,:].*filter
    end
    return vec(x)
end

model_topmute(n::Tuple{Int64,Int64}, mute_end::Array{Float32, 2}, length, x) = vec(mute_end) .* vec(x)

model_topmute(n::Tuple{Int64,Int64,Int64}, mute_end::Int64, length::Int64, x_orig) = model_topmute((n[1]*n[2], n[3]), mute_end, length, x_orig)

model_topmute(n::Tuple{Int64,Int64,Int64}, mute_end::Array{Integer,2}, length::Int64, x_orig) = model_topmute((n[1]*n[2], n[3]), vec(mute_end), length, x_orig)

model_topmute(n::Tuple{Int64,Int64,Int64}, mute_end::Array{Float32, 3}, length, x) = vec(mute_end) .* vec(x)

function judiTopmute(n, mute_end, length)
    # JOLI wrapper for model domain topmute
    N = prod(n)
    T = joLinearFunctionFwd_T(N,N,
                             v -> model_topmute(n, mute_end, length, v),
                             w -> model_topmute(n, mute_end, length, w),
                             Float32,Float32,name="Model topmute")
    return T
end

function find_water_bottom(m::AbstractArray{avDT,2};eps = 1e-4) where {avDT}
    #return the indices of the water bottom of a seismic image
    n = size(m)
    idx = zeros(Integer, n[1])
    for j=1:n[1]
        k=1
        while true
            if abs(m[j,k]) > eps
                idx[j] = k
                break
            end
            k += 1
        end
    end
    return idx
end

find_water_bottom(m::AbstractArray{avDT,3};eps = 1e-4) where {avDT} = reshape(find_water_bottom(reshape(m, :, size(m)[end]);eps=eps), size(m,1), size(m,2))

find_water_bottom(m::PhysicalParameter;eps = 1e-4) = find_water_bottom(m.data)

function depth_scaling(m, model)
# Linear depth scaling function for seismic images
    m_out = deepcopy(reshape(m,Int(prod(model.n)/model.n[end]),model.n[end]))
    filter = sqrt.(0f0:Float32(model.d[end]):Float32(model.d[end])*(model.n[end]-1))
    return vec(m_out.*filter')
end

function judiDepthScaling(model)
# JOLI wrapper for the linear depth scaling function
    N = prod(model.n)
    D = joLinearFunctionFwd_T(N, N,
                             v -> depth_scaling(v,model),
                             w -> depth_scaling(w,model),
                             Float32,Float32,name="Depth scaling")
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
"""
    generate_opinion()

Generate an opinion value for an agent.

See also: [`Agent`](@ref)
"""
function generate_opinion()
    return 2 * rand() - 1
end

"""
    generate_inclin_interact(lambda=log(25))

Generate a value for inclination to interact for an agent.

See also: [`Agent`](@ref)
"""
function generate_inclin_interact(
    lambda=log(25)
)
    # this function was adapted from:
    # https://www.johndcook.com/julia_rng.html
    if lambda <= 0.0
        error("mean must be positive")
    end
    return -(1 / lambda) * log(rand())
end

"""
    generate_check_regularity()

Generate a check regularity value for an agent.

See also: [`Agent`](@ref)
"""
function generate_check_regularity()
    return 1 - (rand() / 4)^2
end

"""
    generate_desired_input_count(mean)

Generate a desired input count value for an agent from the overall mean.

See also: [`Agent`](@ref)
"""
function generate_desired_input_count(
    mean
)
    # this function was adapted from:
    # https://www.johndcook.com/julia_rng.html
    if mean <= 0.0
        error("Mean must be positive")
    end
    stdev = mean / 4
    u1 = rand()
    u2 = rand()
    r = sqrt(-2.0 * log(u1))
    theta = 2.0 * pi * u2
    return trunc(Int, mean + stdev * r * sin(theta))
end

# suppress output of include()
;


text:
points2D = [SVector(1.0, 1.0), #3
            SVector(2.0, 2.0), #1
            SVector(1.2, 1.7), #2
            SVector(1.7, 1.2), #4
            SVector(1.21, 1.7)] #2

root = create_tree(points2D, BoxTreeOptions())

@test length(root.children) == 4
@test root.children[1].data[1] == 2
@test root.children[2].data[1] == 3
@test root.children[2].data[2] == 5
@test root.children[3].data[1] == 1
@test root.children[4].data[1] == 4

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
using Flux, Zygote, StochasticDiffEq
using LinearAlgebra, Statistics
println("heat neural SDE")
using Test, NeuralPDE

using Random
Random.seed!(100)

println("high-dimensional heat equation")
d = 100 # number of dimensions
x0 = fill(8.0f0,d)
tspan = (0.0f0,2.0f0)
dt = 0.5
m = 50 # number of trajectories (batch size)

g(X) = sum(X.^2)
f(X,u,σᵀ∇u,p,t) = Float32(0.0)
μ_f(X,p,t) = zero(X) #Vector d x 1
σ_f(X,p,t) = Diagonal(ones(Float32,d)) #Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls = 10 + d #hidden layer size
#sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))
# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d))
opt = Flux.ADAM(0.005)  #optimizer
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)

ans = solve(prob, pdealg, verbose=true, maxiters=500, trajectories=m,
                            alg=EM(), dt=dt, pabstol = 1f-6)

u_analytical(x,t) = sum(x.^2) .+ d*t
analytical_ans = u_analytical(x0, tspan[end])