text:
# taken liberally from ECOS.jl
# differences from ECOS.jl's implementation:
#    no support for exponential cone
#    final target is a dense matrix, not a sparse matrix

using MathOptInterface
using Compat.SparseArrays

const MOI = MathOptInterface
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

const MOIU = MOI.Utilities

const SF = Union{MOI.SingleVariable, MOI.ScalarAffineFunction{Float64}, MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}}
const SS = Union{MOI.EqualTo{Float64}, MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone}

# Used to build the data with allocate-load during `copy_to`.
# When `optimize!` is called, a the data is used to build `ECOSMatrix`
# and the `ModelData` struct is discarded
mutable struct ModelData
    m::Int # Number of rows/constraints
    n::Int # Number of cols/variables
    IA::Vector{Int} # List of conic rows
    JA::Vector{Int} # List of conic cols
    VA::Vector{Float64} # List of conic coefficients
    b::Vector{Float64} # List of conic coefficients
    IG::Vector{Int} # List of equality rows
    JG::Vector{Int} # List of equality cols
    VG::Vector{Float64} # List of equality coefficients
    h::Vector{Float64} # List of equality coefficients
    objconstant::Float64 # The objective is min c'x + objconstant
    c::Vector{Float64}
end

# This is tied to ECOS's internal representation
mutable struct ConeData
    f::Int # number of linear equality constraints
    l::Int # length of LP cone
    q::Int # length of SOC cone
    qa::Vector{Int} # array of second-order cone constraints
    # The following four field store model information needed to compute `ConstraintPrimal` and `ConstraintDual`
    eqsetconstant::Dict{Int, Float64}   # For the constant of EqualTo
    eqnrows::Dict{Int, Int}             # The number of rows of Zeros
    ineqsetconstant::Dict{Int, Float64} # For the constant of LessThan and GreaterThan
    ineqnrows::Dict{Int, Int}           # The number of rows of each vector sets except Zeros
    function ConeData()
        new(0, 0, 0, Int[],
            Dict{Int, Float64}(),
            Dict{Int, UnitRange{Int}}(),
            Dict{Int, Float64}(),
            Dict{Int, UnitRange{Int}}())
    end
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    cone::ConeData
    maxsense::Bool
    data::Union{Nothing, ModelData} # only non-Nothing between MOI.copy_to and MOI.optimize!
    sol::Union{Nothing, State}
    options
    function Optimizer(; kwargs...)
        new(ConeData(), false, nothing, nothing, kwargs)
    end
end

MOI.get(::Optimizer, ::MOI.SolverName) = "SOCP.jl"
MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.is_empty(optimizer::Optimizer)
    !optimizer.maxsense && optimizer.data === nothing
end
function MOI.empty!(optimizer::Optimizer)
    optimizer.maxsense = false
    optimizer.data = nothing # It should already be nothing except if an error is thrown inside copy_to
    optimizer.sol = nothing
end
function MOI.supports(::Optimizer,
                      ::Union{MOI.ObjectiveSense,
                              MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
    return true
end

function MOI.supports_constraint(::Optimizer,
                                 ::Type{MOI.VectorAffineFunction{Float64}},
                                 ::Type{<:Union{MOI.Zeros, MOI.Nonnegatives,
                                                MOI.SecondOrderCone}})
    return true
end

MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end
# Computes cone dimensions
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, MOI.Zeros}) = ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.Zeros)
    ci = cone.f
    cone.f += MOI.dimension(s)
    ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, MOI.Nonnegatives}) = ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.Nonnegatives)
    ci = cone.l
    cone.l += MOI.dimension(s)
    ci
end
constroffset(cone::ConeData, ci::CI{<:MOI.AbstractFunction, <:MOI.SecondOrderCone}) = cone.l + ci.value
function _allocate_constraint(cone::ConeData, f, s::MOI.SecondOrderCone)
    push!(cone.qa, s.dimension)
    ci = cone.q
    cone.q += MOI.dimension(s)
    ci
end
constroffset(optimizer::Optimizer, ci::CI) = constroffset(optimizer.cone, ci::CI)
function MOIU.allocate_constraint(optimizer::Optimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    CI{F, S}(_allocate_constraint(optimizer.cone, f, s))
end

# Build constraint matrix
output_index(t::MOI.VectorAffineTerm) = t.output_index
variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
coefficient(t::MOI.ScalarAffineTerm) = t.coefficient
coefficient(t::MOI.VectorAffineTerm) = coefficient(t.scalar_term)
constrrows(s::MOI.AbstractVectorSet) = 1:MOI.dimension(s)
constrrows(optimizer::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, MOI.Zeros}) = 1:optimizer.cone.eqnrows[constroffset(optimizer, ci)]
constrrows(optimizer::Optimizer, ci::CI{<:MOI.AbstractVectorFunction, <:MOI.AbstractVectorSet}) = 1:optimizer.cone.ineqnrows[constroffset(optimizer, ci)]
matrix(data::ModelData, s::MOI.Zeros) = data.b, data.IA, data.JA, data.VA
matrix(data::ModelData, s::Union{MOI.Nonnegatives, MOI.SecondOrderCone}) = data.h, data.IG, data.JG, data.VG
matrix(optimizer::Optimizer, s) = matrix(optimizer.data, s)
# ECOS orders differently than MOI the second and third dimension of the exponential cone
orderval(val, s) = val
orderidx(idx, s) = idx
expmap(i) = (1, 3, 2)[i]
function MOIU.load_constraint(optimizer::Optimizer, ci::MOI.ConstraintIndex, f::MOI.VectorAffineFunction, s::MOI.AbstractVectorSet)
    func = MOIU.canonical(f)
    I = Int[output_index(term) for term in func.terms]
    J = Int[variable_index_value(term) for term in func.terms]
    V = Float64[-coefficient(term) for term in func.terms]
    offset = constroffset(optimizer, ci)
    rows = constrrows(s)
    if s isa MOI.Zeros
        optimizer.cone.eqnrows[offset] = length(rows)
    else
        optimizer.cone.ineqnrows[offset] = length(rows)
    end
    i = offset .+ rows
    # The ECOS format is b - Ax ∈ cone
    # so minus=false for b and minus=true for A
    b, Is, Js, Vs = matrix(optimizer, s)
    b[i] .= orderval(f.constants, s)
    append!(Is, offset .+ orderidx(I, s))
    append!(Js, J)
    append!(Vs, V)
end

function MOIU.allocate_variables(optimizer::Optimizer, nvars::Integer)
    optimizer.cone = ConeData()
    VI.(1:nvars)
end

function MOIU.load_variables(optimizer::Optimizer, nvars::Integer)
    cone = optimizer.cone
    m = cone.l + cone.q
    IA = Int[]
    JA = Int[]
    VA = Float64[]
    b = zeros(cone.f)
    IG = Int[]
    JG = Int[]
    VG = Float64[]
    h = zeros(m)
    c = zeros(nvars)
    optimizer.data = ModelData(m, nvars, IA, JA, VA, b, IG, JG, VG, h, 0., c)
end

function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    optimizer.maxsense = sense == MOI.MAX_SENSE
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction,
                       ::MOI.ScalarAffineFunction{Float64})
end

function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction,
                   f::MOI.ScalarAffineFunction)
    c0 = Vector(sparsevec(variable_index_value.(f.terms), coefficient.(f.terms),
                          optimizer.data.n))
    optimizer.data.objconstant = f.constant
    optimizer.data.c = optimizer.maxsense ? -c0 : c0
    return nothing
end

function MOI.optimize!(instance::Optimizer)
    if instance.data === nothing
        # optimize! has already been called and no new model has been copied
        return
    end
    cone = instance.cone
    m = instance.data.m
    n = instance.data.n
    A = sparse(instance.data.IA, instance.data.JA, instance.data.VA, cone.f, n)
    b = instance.data.b
    G = sparse(instance.data.IG, instance.data.JG, instance.data.VG, m, n)
    h = instance.data.h
    objconstant = instance.data.objconstant
    c = instance.data.c

    cones = Socp.Cone[Socp.POC(0,cone.l)]
    offs = cone.l
    for q in cone.qa
    	push!(cones, Socp.SOC(offs, q))
    	offs += q
    end
    prob = Socp.Problem(c, A, b, G, h, cones)
    res = solve_socp(prob, Socp.SparseSolver(prob))
    instance.sol = res
end

#=
using MathOptInterface
const MOI = MathOptInterface
model = MOI.instantiate(() -> Socp.Optimizer(maxit=10000); with_bridge_type=Float64)
x = MOI.add_variable(model)
y = MOI.add_variable(model)

MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0,-5.0], [x,y]), 0.0))
MOI.add_constraint(model, MOI.SingleVariable(x), MOI.Interval(100.0, 200.0))
MOI.add_constraint(model, MOI.SingleVariable(y), MOI.Interval(80.0, 170.0))
MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], [x,y]), 0.0), MOI.GreaterThan(200.0))
=#

const reorderval = orderval
function MOI.get(instance::Optimizer, ::MOI.VariablePrimal, vi::VI)
    instance.sol.x[vi.value]
end
MOI.get(instance::Optimizer, a::MOI.VariablePrimal, vi::Vector{VI}) = MOI.get.(instance, Ref(a), vi)

# setconstant: Retrieve set constant stored in `ConeData` during `copy_to`
setconstant(instance::Optimizer, offset, ::CI{<:MOI.AbstractFunction, <:MOI.EqualTo}) = instance.cone.eqsetconstant[offset]
setconstant(instance::Optimizer, offset, ::CI) = instance.cone.ineqsetconstant[offset]
_unshift(instance::Optimizer, offset, value, ::CI) = value
_unshift(instance::Optimizer, offset, value, ci::CI{<:MOI.AbstractScalarFunction, <:MOI.AbstractScalarSet}) = value + setconstant(instance, offset, ci)
function MOI.get(instance::Optimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, MOI.Zeros})
    rows = constrrows(instance, ci)
    zeros(length(rows))
end
function MOI.get(instance::Optimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, <:MOI.EqualTo})
    offset = constroffset(instance, ci)
    setconstant(instance, offset, ci)
end
function MOI.get(instance::Optimizer, ::MOI.ConstraintPrimal, ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(instance, ci)
    rows = constrrows(instance, ci)
    _unshift(instance, offset, scalecoef(rows, reorderval(instance.sol.s[offset .+ rows], S), false, S), ci)
end

_dual(instance, ci::CI) = instance.sol.z
function MOI.get(instance::Optimizer, ::MOI.ConstraintDual, ci::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    offset = constroffset(instance, ci)
    rows = constrrows(instance, ci)
    scalecoef(rows, reorderval(_dual(instance, ci)[offset .+ rows], S), false, S)
end

MOI.get(instance::Optimizer, ::MOI.ResultCount) = 1


text:
# Automatically generated using Clang.jl

@cenum ImPlotMarker::Cint begin
    ImPlotMarker_None = 1
    ImPlotMarker_Circle = 2
    ImPlotMarker_Square = 4
    ImPlotMarker_Diamond = 8
    ImPlotMarker_Up = 16
    ImPlotMarker_Down = 32
    ImPlotMarker_Left = 64
    ImPlotMarker_Right = 128
    ImPlotMarker_Cross = 256
    ImPlotMarker_Plus = 512
    ImPlotMarker_Asterisk = 1024
end

struct ImPlotStyle
    LineWeight::Cfloat
    Marker::ImPlotMarker
    MarkerSize::Cfloat
    MarkerWeight::Cfloat
    ErrorBarSize::Cfloat
    ErrorBarWeight::Cfloat
    DigitalBitHeight::Cfloat
    DigitalBitGap::Cfloat
    Colors::NTuple{14, Cint}
end

struct ImPlotRange
    Min::Cdouble
    Max::Cdouble
end

struct ImPlotLimits
    X::ImPlotRange
    Y::ImPlotRange
end

struct ImPlotPoint
    x::Cdouble
    y::Cdouble
end

@cenum ImPlotFlags::Cint begin
    ImPlotFlags_MousePos = 1
    ImPlotFlags_Legend = 2
    ImPlotFlags_Highlight = 4
    ImPlotFlags_BoxSelect = 8
    ImPlotFlags_Query = 16
    ImPlotFlags_ContextMenu = 32
    ImPlotFlags_Crosshairs = 64
    ImPlotFlags_CullData = 128
    ImPlotFlags_AntiAliased = 256
    ImPlotFlags_NoChild = 512
    ImPlotFlags_YAxis2 = 1024
    ImPlotFlags_YAxis3 = 2048
    ImPlotFlags_Default = 175
end

@cenum ImPlotAxisFlags::Cint begin
    ImPlotAxisFlags_NULL = 0
    ImPlotAxisFlags_GridLines = 1
    ImPlotAxisFlags_TickMarks = 2
    ImPlotAxisFlags_TickLabels = 4
    ImPlotAxisFlags_Invert = 8
    ImPlotAxisFlags_LockMin = 16
    ImPlotAxisFlags_LockMax = 32
    ImPlotAxisFlags_Adaptive = 64
    ImPlotAxisFlags_LogScale = 128
    ImPlotAxisFlags_Scientific = 256
    ImPlotAxisFlags_Default = 71
    ImPlotAxisFlags_Auxiliary = 70
end

@cenum ImPlotCol::Cint begin
    ImPlotCol_Line = 0
    ImPlotCol_Fill = 1
    ImPlotCol_MarkerOutline = 2
    ImPlotCol_MarkerFill = 3
    ImPlotCol_ErrorBar = 4
    ImPlotCol_FrameBg = 5
    ImPlotCol_PlotBg = 6
    ImPlotCol_PlotBorder = 7
    ImPlotCol_XAxis = 8
    ImPlotCol_YAxis = 9
    ImPlotCol_YAxis2 = 10
    ImPlotCol_YAxis3 = 11
    ImPlotCol_Selection = 12
    ImPlotCol_Query = 13
    ImPlotCol_COUNT = 14
end

@cenum ImPlotStyleVar::Cint begin
    ImPlotStyleVar_LineWeight = 0
    ImPlotStyleVar_Marker = 1
    ImPlotStyleVar_MarkerSize = 2
    ImPlotStyleVar_MarkerWeight = 3
    ImPlotStyleVar_ErrorBarSize = 4
    ImPlotStyleVar_ErrorBarWeight = 5
    ImPlotStyleVar_DigitalBitHeight = 6
    ImPlotStyleVar_DigitalBitGap = 7
    ImPlotStyleVar_COUNT = 8
end

@cenum ImPlotColormap::Cint begin
    ImPlotColormap_Default = 0
    ImPlotColormap_Dark = 1
    ImPlotColormap_Pastel = 2
    ImPlotColormap_Paired = 3
    ImPlotColormap_Viridis = 4
    ImPlotColormap_Plasma = 5
    ImPlotColormap_Hot = 6
    ImPlotColormap_Cool = 7
    ImPlotColormap_Pink = 8
    ImPlotColormap_Jet = 9
    ImPlotColormap_COUNT = 10
end

ext: jl


text:
# ### [`ToEltype`](#)

"""
    ToEltype(T)

Converts any `AbstractArrayItem` to an `AbstractArrayItem{N, T}`.

Supports `apply!`.

## Examples

{cell=ToEltype}
```julia
using DataAugmentation

tfm = ToEltype(Float32)
item = ArrayItem(rand(Int, 10))
apply(tfm, item)
```
"""
struct ToEltype{T} <: Transform end
ToEltype(T::Type) = ToEltype{T}()

apply(::ToEltype{T}, item::AbstractArrayItem{N, <:T}; randstate = nothing) where {N, T} = item
function apply(::ToEltype{T1}, item::AbstractArrayItem{N, T2}; randstate = nothing) where {N, T1, T2}
    newdata = map(x -> convert(T1, x), itemdata(item))
    item = setdata(item, newdata)
    return item
end

function apply!(buf, ::ToEltype, item::AbstractArrayItem; randstate = nothing)
    ## copy! does type conversion under the hood
    copy!(itemdata(buf), itemdata(item))
    return buf
end

# ### [`Normalize`](#)

"""
    Normalize(means, stds)

Normalizes the last dimension of an `AbstractArrayItem{N}`.

Supports `apply!`.

## Examples

Preprocessing a 3D image with 3 color channels.

{cell=Normalize}
```julia
using DataAugmentation, Images
image = Image(rand(RGB, 20, 20, 20))
tfms = ImageToTensor() |> Normalize((0.1, -0.2, -0.1), (1,1,1.))
apply(tfms, image)
```

"""
struct Normalize{N} <: Transform
    means::SVector{N}
    stds::SVector{N}
end

function Normalize(means, stds)
    length(means) == length(stds) || error("`means` and `stds` must have same length")
    N = length(means)
    return Normalize{N}(SVector{N}(means), SVector{N}(stds))
end

function apply(tfm::Normalize, item::ArrayItem{N, T}; randstate = nothing) where {N, T}
    return ArrayItem(normalize(itemdata(item), tfm.means, tfm.stds))
end

function apply!(
        buf,
        tfm::Normalize,
        item::ArrayItem;
        randstate = nothing)
    copy!(itemdata(buf), itemdata(item))
    normalize!(itemdata(buf), tfm.means, tfm.stds)
    return buf
end

function normalize!(
        a::AbstractArray{T, N},
        means::SVector{M},
        stds::SVector{M}) where {N, T, M}
    means = reshape(convert.(T, means), (1 for _ = 2:N)..., M)
    stds = reshape(convert.(T, stds), (1 for _ = 2:N)..., M)
    a .-= means
    a ./= stds
    return a
end

normalize(a, means, stds) = normalize!(copy(a), means, stds)

function denormalize!(
        a::AbstractArray{T, N},
        means::SVector{M},
        stds::SVector{M}) where {N, T, M}
    means = reshape(convert.(T, means), (1 for _ = 2:N)..., M)
    stds = reshape(convert.(T, stds), (1 for _ = 2:N)..., M)
    a .*= stds
    a .+= means
    return a
end

denormalize(a, means, stds) = denormalize!(copy(a), means, stds)

# ### [`NormalizeIntensity`]

"""
    NormalizeIntensity()

Normalizes the pixels of an array based on calculated mean and std.
"""

struct NormalizeIntensity <: Transform end

function apply(tfm::NormalizeIntensity, item::ArrayItem; randstate = nothing)
    array = itemdata(item)
    slices = ones(Bool, size(array))
    means = mean(array[slices])
    stds = std(array[slices])
    array[slices] = (array[slices] .- means) / stds
    return ArrayItem(array)
end

# ### [`ImageToTensor`]

"""
    ImageToTensor()

Expands an `Image{N, T}` of size `sz` to an `ArrayItem{N+1}` with
size `(sz..., ch)` where `ch` is the number of color channels of `T`.

Supports `apply!`.

## Examples

{cell=ImageToTensor}
```julia
using DataAugmentation, Images

image = Image(rand(RGB, 50, 50))
tfm = ImageToTensor()
apply(tfm, image)
```

"""
struct ImageToTensor{T} <: Transform end

ImageToTensor(T::Type{<:Number} = Float32) = ImageToTensor{T}()

function apply(::ImageToTensor{T}, image::Image; randstate = nothing) where T
    return ArrayItem(imagetotensor(itemdata(image), T))
end

function apply!(buf, ::ImageToTensor, image::Image; randstate = nothing)
    imagetotensor!(buf.data, image.data)
    return buf
end

function imagetotensor(image::AbstractArray{C, N}, T = Float32) where {C<:Color, N}
    T.(permuteddimsview(_channelview(image), ((i for i in 2:N+1)..., 1)))
end

#=
function imagetotensor(image::AbstractArray{C, N}, T = Float32) where {TC, C<:Color{TC, 1}, N}
    return T.(_channelview(image))
end
=#

# TODO: relax color type constraint, implement for other colors
# single-channel colors need a `channelview` that also expands the array
function imagetotensor!(buf, image::AbstractArray{<:Color, N}) where N
    permutedims!(
        buf,
        _channelview(image),
        (2:N+1..., 1))
end

function tensortoimage(a::AbstractArray)
    nchannels = size(a)[end]
    if nchannels == 3
        return tensortoimage(RGB, a)
    elseif nchannels == 1
        return tensortoimage(Gray, a)
    else
        error("Found image tensor with $nchannels color channels. Pass in color type
                explicitly.")
    end
end

function tensortoimage(C::Type{<:Color}, a::AbstractArray{T, N}) where {T, N}
    perm = (N, 1:N-1...)
    return _colorview(C, permuteddimsview(a, perm))
end

function _channelview(img)
    chview = channelview(img)
    # for single-channel colors, expand the color dimension anyway
    if size(img) == size(chview)
        chview = reshape(chview, 1, size(chview)...)
    end
    return chview
end

function _colorview(C::Type{<:Color}, img) where T
    if size(img, 1) == 1
        img = reshape(img, size(img)[2:end])
    end
    return colorview(C, img)
end

# OneHot encoding

"""
    OneHot([T = Float32])

One-hot encodes a `MaskMulti` with `n` classes and size `sz` into
an array item of size `(sz..., n)` with element type `T`. Supports [`apply!`](#).

```julia
item = MaskMulti(rand(1:4, 100, 100), 1:4)
apply(OneHot(), item)
```
"""
struct OneHot{T} <: Transform end
OneHot() = OneHot{Float32}()

function apply(tfm::OneHot{T}, item::MaskMulti; randstate = nothing) where T
    mask = parent(itemdata(item))
    a = zeros(T, size(mask)..., length(item.classes))
    for I in CartesianIndices(mask)
        a[I, mask[I]] = one(T)
    end

    return ArrayItem(a)
end

function apply!(buf, tfm::OneHot{T}, item::MaskMulti; randstate = nothing) where T
    mask = parent(itemdata(item))
    a = itemdata(buf)
    fill!(a, zero(T))

    for I in CartesianIndices(mask)
        a[I, mask[I]] = one(T)
    end

    return buf
end

function onehot(T, x::Int, n::Int)
    v = fill(zero(T), n)
    v[x] = one(T)
    return v
end
onehot(x, n) = onehot(Float32, x, n)


text:
using Test, MicrostructureNoise

@test MicrostructureNoise.piecewise([0,1,2],[10,11]) == ([0, 1, 1, 2], [10, 10, 11, 11])


text:
module Hive

include("sql_exports.jl")
include("sql_imports.jl") # DBMS Structured SubQuery _to_sql _placeholder _placeholders

const DatabaseID = DBMS.Hive

"""
    to_sql(query::Structured)::String
"""
to_sql(query::Structured)::String = _to_sql(DatabaseID(), query)

"""
    to_sql(subquery::SubQuery)::String
"""
to_sql(subquery::SubQuery)::String = _to_sql(DatabaseID(), subquery)

placeholder(nth::Int) = _placeholder(DatabaseID(), nth)
placeholders(dims::Int) = _placeholders(DatabaseID(), dims)

end # Octo.Adapters.Hive


text:
"""
    sweep_contract_dangling(LTN::LabelledTensorNetwork, χ, τ;
        fast=false, valid=false, planar=false, connected=false, report=false)
    sweep_contract_dangling(TN::TensorNetwork, χ, τ;
        fast=false, valid=false, planar=false, connected=false, report=false)

A modified version of `sweep_contract` which can be used to evaluate open tensor networks
which possess dangling legs. This is done by halting the sweep procedure before the final
tensor is contracted and returning the MPS approximating the network so far. To do this a
dummy tensor must be present in the network which is above the rest of the network, i.e. has
a higher y coordinate.

Similar to sweep_contract underflow/overflow issues are avoided by return the MPS in a
floating point format. The output is returned as a tuple `(mps::MPS, i::Int64)` where `mps`
has a 2-norm in [1,2), representing the MPS `mps*2^i`.

For other details see the documentation of `sweep_contract`.
"""
sweep_contract_dangling(LTN::LabelledTensorNetwork, χ::Int, τ::Int;
    fast=false, valid=false, planar=false, connected=false, report=false) =
sweep_contract_dangling!(deepcopy(LTN), χ, τ;
    fast=fast, planar=planar, connected=connected, report=report)

sweep_contract_dangling(TN::TensorNetwork, χ::Int, τ::Int;
    fast=false, valid=false, planar=false, connected=false, report=false) =
sweep_contract_dangling!(deepcopy(TN), χ, τ;
    fast=fast, planar=planar, connected=connected, report=report)

"""
    sweep_contract_dangling!(LTN::LabelledTensorNetwork, χ, τ;
        fast=false, valid=false, planar=false, connected=false, report=false)
    sweep_contract_dangling!(TN::TensorNetwork, χ, τ;
        fast=false, valid=false, planar=false, connected=false, report=false)

The mutating form of `sweep_contract_dangling`.
"""
sweep_contract_dangling!(LTN::LabelledTensorNetwork, χ::Int, τ::Int;
    fast=false, valid=false, planar=false, connected=false, report=false) =
sweep_contract_dangling!(delabel(LTN), χ, τ;
        fast=fast, planar=planar, connected=connected, report=report)

function sweep_contract_dangling!(TN::TensorNetwork, χ::Int, τ::Int;
    fast=false,valid=false,planar=false,connected=false,report=false)::Tuple{MPS,Int}
    if !fast
        valid || checkvalid(TN)
        planar || planarise!(TN)
        connected || connect!(hull!(TN))
    end

    sort!(TN)

    resexp = 0
    count = 0

    MPS_t = [ones(1,1,1)]
    MPS_i = Int[]

    for (i,t) ∈ enumerate(TN)
        if i==length(TN)
            break;
        end
        ind_up = Int[]
        ind_do = Int[]
        for n ∈ t.adj
            if TN[n]>t
                push!(ind_up, n)
            elseif TN[n]<t
                push!(ind_do, n)
            else
                throw(InvalidTNError("Overlapping tensors"))
            end
        end
        sort!(ind_up, by=λ->atan(TN[λ].x-t.x,TN[λ].y-t.y))
        sort!(ind_do, by=λ->atan(TN[λ].x-t.x,t.y-TN[λ].y))
        σ = permutebetween(t.adj, [ind_do; ind_up])
        t.arr = permutedims(t.arr, σ)
        s = size(t.arr)
        t.arr = reshape(t.arr,(prod(s[1:length(ind_do)]),s[length(ind_do)+1:end]...))

        if isempty(MPS_i)
            MPS_t = splitMPStensor(MPS_t[1][1]*reshape(t.arr,(size(t.arr)...,1)))
            MPS_i = ind_up
        else
            lo = findfirst(isequal(i), MPS_i)
            hi = findlast(isequal(i), MPS_i)

            isnothing(lo) && throw(InvalidTNError("Disconnected TN"))

            X::Array{Float64} = MPS_t[lo]
            for j ∈ lo+1:hi
                finalsize = (size(X,1),size(X,2)*size(MPS_t[j],2),size(MPS_t[j],3))
                X = reshape(X,(size(X,1)*size(X,2),size(X,3)))*
                    reshape(MPS_t[j],(size(MPS_t[j],1),size(MPS_t[j],2)*size(MPS_t[j],3)))
                X = reshape(X,finalsize)
            end
            X = permutedims(X,[1,3,2])
            M = reshape(t.arr,(size(t.arr,1),prod(size(t.arr)[2:end])))
            X = reshape(
                reshape(X,(size(X,1)*size(X,2),size(X,3)))*M,
                (size(X,1),size(X,2),size(M,2))
            )
            X = permutedims(X,[1,3,2])
            X = reshape(X,(size(X,1),size(t.arr)[2:end]...,size(X,3)))

            MPS_i = [MPS_i[1:lo-1]; ind_up; MPS_i[hi+1:end]]
            if ndims(X)!=2
                MPS_t = [MPS_t[1:lo-1]; splitMPStensor(X); MPS_t[hi+1:end]]
            elseif isempty(MPS_i)
                MPS_t=[reshape([X[1]],(1,1,1))]
            elseif lo>1
                s = size(MPS_t[lo-1])
                MPS_t[lo-1] = reshape(
                    reshape(MPS_t[lo-1],(s[1]*s[2],s[3]))*X,
                    (s[1],s[2],size(X,2))
                )
                MPS_t = [MPS_t[1:lo-1]; MPS_t[hi+1:end]]
            else
                s = size(MPS_t[hi+1])
                MPS_t[hi+1] = reshape(
                    X*reshape(MPS_t[hi+1],(s[1],s[2]*s[3])),
                    (size(X,1),s[2],s[3])
                )
                MPS_t = [MPS_t[1:lo-1]; MPS_t[hi+1:end]]
            end

            if any(size.(MPS_t,3).>τ)
                count += 1
                truncMPS!(MPS_t, χ)
                h = Int(floor(log2(LinearAlgebra.norm(MPS_t[1]))))
                resexp += h
                MPS_t[1] /= exp2(h)
            end
        end
    end

    report && println("Number of truncations: $count")

    truncMPS!(MPS_t, χ)
    h = Int(floor(log2(LinearAlgebra.norm(MPS_t[1]))))
    resexp += h
    MPS_t[1] /= exp2(h)

    return (MPS_t,resexp);
end


text:
abstract type AbstractBoundaryCondition end

struct DefaultBC <: AbstractBoundaryCondition end

Base.@kwdef struct BulkFormulaTemperature{𝒯, 𝒰, 𝒱} <: AbstractBoundaryCondition
    drag_coef_temperature::𝒯
    drag_coef_moisture::𝒰
    surface_temperature::𝒱
end

Base.@kwdef struct CoupledPrimarySlabOceanBC{𝒯, 𝒰} <: AbstractBoundaryCondition
    drag_coef_temperature::𝒯
    drag_coef_moisture::𝒰
end

abstract type TemperatureBC end
struct Insulating <: TemperatureBC end
struct CoupledSecondaryAtmosModelBC <: TemperatureBC end

function numerical_boundary_flux_first_order!(
    numerical_flux::NumericalFluxFirstOrder,
    ::DefaultBC,
    balance_law::DryAtmosModel,
    fluxᵀn::Vars{S},
    n̂::SVector,
    state⁻::Vars{S},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    aux⁺::Vars{A},
    t,
    direction,
    state1⁻::Vars{S},
    aux1⁻::Vars{A},
) where {S, A}
    state⁺.ρ = state⁻.ρ
    state⁺.ρe = state⁻.ρe
    state⁺.ρq = state⁻.ρq

    # project and reflect for impenetrable condition, but
    # leave tangential component untouched
    ρu⁻ = state⁻.ρu
    state⁺.ρu = ρu⁻ - 2n̂ ⋅ ρu⁻ .* SVector(n̂)

    numerical_flux_first_order!(numerical_flux, balance_law, fluxᵀn, n̂, state⁻, aux⁻, state⁺, aux⁺, t, direction)
end

function numerical_boundary_flux_first_order!(
    numerical_flux::NumericalFluxFirstOrder,
    bctype::BulkFormulaTemperature,
    model::DryAtmosModel,
    fluxᵀn::Vars{S},
    n̂::SVector,
    state⁻::Vars{S},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    aux⁺::Vars{A},
    t,
    direction,
    state1⁻::Vars{S},
    aux1⁻::Vars{A},
) where {S, A}
    # Impenetrable free-slip condition to reflect and project momentum
    # at the boundary
    numerical_boundary_flux_first_order!(
        numerical_flux,
        DefaultBC(),
        model,
        fluxᵀn,
        n̂,
        state⁻,
        aux⁻,
        state⁺,
        aux⁺,
        t,
        direction,
        state1⁻,
        aux1⁻,
    )

    E, H = calc_ocean_sfc_fluxes(model.physics, bctype, state⁻, aux⁻)
    LH_v0 = model.physics.parameters.LH_v0

    fluxᵀn.ρ -= E / LH_v0
    fluxᵀn.ρe -= E + H
    fluxᵀn.ρq -= E / LH_v0
end

function numerical_boundary_flux_first_order!(
    numerical_flux::NumericalFluxFirstOrder,
    bctype::CoupledPrimarySlabOceanBC,
    model::DryAtmosModel,
    fluxᵀn::Vars{S},
    n̂::SVector,
    state⁻::Vars{S},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    aux⁺::Vars{A},
    t,
    direction,
    state1⁻::Vars{S},
    aux1⁻::Vars{A},
) where {S, A}
    # Impenetrable free-slip condition to reflect and project momentum
    # at the boundary

    numerical_boundary_flux_first_order!(
        numerical_flux,
        DefaultBC(),
        model,
        fluxᵀn,
        n̂,
        state⁻,
        aux⁻,
        state⁺,
        aux⁺,
        t,
        direction,
        state1⁻,
        aux1⁻,
    )

    # The following will be moved the the second-order kernel (as a Neumann BC) once available
    LH_v0 = model.physics.parameters.LH_v0
    E, H = calc_ocean_sfc_fluxes(model.physics, bctype, state⁻, aux⁻) #[W/m^2]

    fluxᵀn.ρ -= E / LH_v0
    fluxᵀn.ρe -= E + H
    fluxᵀn.ρq -= E / LH_v0

end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
mutable struct pasty
	_id::String
	language::String
	title::String
	code::String
end

struct edit
	_id::String
	editId::Number
	editType::Number
	metadata::Vector{String}
	edit::String
	editedAt::Number
end

mutable struct pasteEditInfo
	title::String
	isPrivate::Bool
	isPublic::Bool
	tags::Vector{String}
	pasties::Vector{pasty}
end

mutable struct paste
	_id::String
	ownerId::String
	title::String
	createdAt::UInt64
	expiresIn::String
	deletesAt::UInt64
	stars::UInt64
	isPrivate::Bool
	isPublic::Bool
	encrypted::Bool
	tags::Vector{String}
	pasties::Vector{pasty}
	edits::Vector{edit}
end

struct PastyCreateInfo
	title::String
	language::String
	code::String
end

struct PasteCreateInfo
	title::String # Optional
	expiresIn::String # Optional
	isPrivate::Bool # Optional
	isPublic::Bool # Optional
	tags::Vector{String} # Tags
	pasties::Vector{PastyCreateInfo} # Mandatory
end


text:
pyplot()
cd(pwd)
dir = "results/"

# Plot parameter recovery
recoveryPlots = plotrecovery(results, (mu=0, sigma=1), (:sampler, :Nd); save=true, dir=dir)

# Plot mean run time as a function of number of data points (Nd) for each sampler
meantimePlot = plotsummary(results, :Nd, :time, (:sampler,); save=true, dir=dir, yscale=:log10)

# Plot mean allocations as a function of number of data points (Nd) for each sampler
meanallocPlot = plotsummary(results, :Nd, :allocations, (:sampler,); save=false, dir=dir, yscale=:log10,
  ylabel="Allocations (log scale)"
)

# Plot mean ess per second as a function of number of data points (Nd) for each sampler
meanESSPSPlot = plotsummary(results, :Nd, :ess_ps, (:sampler,); save=true, dir=dir)

# Plot mean ess as a function of data points (Nd) for each sampler
meanESSPlot = plotsummary(results, :Nd, :ess, (:sampler,); save=true, dir=dir)

# Plot density of effective sample size as function of number of data points (Nd) for each sampler
essPlots = plotdensity(results, :ess, (:sampler, :Nd); save=true, dir=dir)

# Plot density of rhat as function of number of data points (Nd) for each sampler
rhatPlots = plotdensity(results, :r_hat, (:sampler, :Nd); save=true, dir=dir)

# Plot density of time as function of number of data points (Nd) for each sampler
timePlots = plotdensity(results, :time, (:sampler, :Nd); save=true, dir=dir)

# Plot density of gc time percent as function of number of data points (Nd) for each sampler
gcPlots = plotdensity(results, :gcpercent, (:sampler, :Nd); save=true, dir=dir)

# Plot density of memory allocations as function of number of data points (Nd) for each sampler
memPlots = plotsummary(results, :Nd, :allocations, (:sampler,); save=true, dir=dir, xscale=:log10,
  xlabel="Allocations (log scale)"
)

# Plot density of memory allocations as function of number of data points (Nd) for each sampler
memPlots = plotdensity(results, :allocations, (:sampler, :Nd); save=true, dir=dir, xscale=:log10,
  xlabel="Allocations (log scale)"
)

# Plot density of megabytes allocated as function of number of data points (Nd) for each sampler
megPlots = plotdensity(results, :megabytes, (:sampler, :Nd); save=true, dir=dir)

# Scatter plot of epsilon and effective sample size as function of number of data points (Nd) for each sampler
scatterPlots = plotscatter(results, :epsilon, :ess, (:sampler, :Nd); save=true, dir=dir)

# Density plot for tree_depth
depthPlots = plotdensity(results, :tree_depth,(:sampler, :Nd); save=true, dir=dir)
epsPlots = plotdensity(results, :epsilon,(:sampler, :Nd); save=true, dir=dir)

# Plot mean ess per second of number of data points (Nd) for each sampler
ess = plotsummary(results, :Nd, :ess, (:sampler,); save=true, dir=dir)


text:
function write_vtk(vector::Array{Float64, 1}, name::String, dofmap::DofMap, filename::String, mode::String)

    UFL_TO_VTK = Dict(ufl.quadrilateral => [[1, 2, 4, 3], 9])

    num_cells = num_entities(dofmap.mesh, 2)
    num_vertices = num_entities(dofmap.mesh, 0)
    vert_per_cell = length(dofmap.mesh.topology[(2, 0)][1, :])

    perm = UFL_TO_VTK[dofmap.element.cell()][1]
    vtk_cell_id = UFL_TO_VTK[dofmap.element.cell()][2]
    gdim = length(dofmap.mesh.vertices[1, :])

    open(filename, mode) do f
        write(f, "# vtk DataFile Version 2.0\n")
        write(f, "Minidolfin.jl output\n")
        write(f, "ASCII\n")
        write(f, "DATASET UNSTRUCTURED_GRID\n")

        write(f, "POINTS $num_vertices float\n")
        for vertices in eachrow(dofmap.mesh.vertices)
            write(f, "$(join(vcat(vertices, zeros(3-gdim)), " "))\n")
        end

        write(f, "CELLS $num_cells $((vert_per_cell + 1) * num_cells)\n")
        for vertices_id in eachrow(dofmap.mesh.topology[(2, 0)])
            # -1 for Julia index shift
            write(f, "$vert_per_cell $(join(view(vertices_id .- 1, perm), " "))\n")
        end

        write(f, "CELL_TYPES $num_cells\n")
        for cell_id in 1:num_cells
            write(f, "$vtk_cell_id\n")
        end

        write(f, "POINT_DATA $num_vertices\n")
        write(f, "SCALARS $name double 1\n")
        write(f, "LOOKUP_TABLE default\n")
        for val in vector
            write(f, "$val\n")
        end
    end

end


text:
module CompileModules

export compilecache, loadcache, ensure_compilecache

using Logging

defaultcache(input::AbstractString) = input * ".ji"

# This is a lightly modified copy of compilecache from Base/loading.jl
function compilecache(input::AbstractString,
                      output::AbstractString=defaultcache(input))
    concrete_deps = copy(Base._concrete_dependencies)
    for (key, mod) in Base.loaded_modules
        if !(mod === Main || mod === Core || mod === Base)
            push!(concrete_deps, key => Base.module_build_id(mod))
        end
    end

    verbosity = isinteractive() ? Logging.Info : Logging.Debug
    if isfile(output)
        @logmsg verbosity "Recompiling cache file $output for $input"
    else
        @logmsg verbosity "Precompiling $input"
    end
    p = Base.create_expr_cache(input, output, concrete_deps, nothing)
    if success(p)
        # append checksum to the end of the .ji file:
        open(output, "a+") do f
            write(f, Base._crc32c(seekstart(f)))
        end
    elseif p.exitcode == 125
        throw(Base.PrecompilableError())
    else
        error("Failed to precompile $input to $output.")
    end
    output
end#function

loadcache(cache::AbstractString=defaultcache(PROGRAM_FILE)) =
    Base._require_from_serialized(cache)

function ensure_compilecache(input::AbstractString,
                             output::AbstractString=defaultcache(input))
    # TODO: handle deps being stale
    if !isfile(output) || Base.stale_cachefile(input, output) === true
        compilecache(input, output)
    end#if
    output
end#function

end#module


text:
using SatelliteToolbox
using LinearAlgebra

#=
"""
shootingtime = shootingtime_decision(satplace_plans, targetplace, currenttime, limit_time, dt)

撮影時刻を決定する関数

# Argments
 - `satplace_plans`：現在時刻~撮影期限までの衛星位置ベクトル@ECEF
 - `targetplace`：撮影対象の位置ベクトル@ECEF
 - `currenttime`：現在時刻（ループ回数）
 - `limit_time`：撮影期限までの時間 （ループ回数　 現在時刻 ~ (currenttime + limit_time)*dt  の間で撮影する）
 - `dt`：シミュレータループ時間間隔
# Return
- `shootingtime`：撮影時刻
- `shooting_vec`：撮影時衛星->撮影地点方向ベクトル@ECEF

# Conditions
 - `cam_viewangle`：カメラの視野角 [deg]
 - `sat_axisval`:衛星の撮影時姿勢の可動域 [deg]
"""

function shootingtime_decision(satplace_plans, targetplace, currenttime, limit_time, cam_viewangle)
    cam_viewangle = 40
    sat_axisval = 30
    shootingtime = 0
    cθ_shooting = 0
    shootingtime = 0
    shooting_vec = zeros(1,3)

    cθ_lim = cosd(cam_viewangle/4 + sat_axisval)
    print("cθ_lim : ")
    println(cθ_lim)
    for i = 1:limit_time
        vec_stot= targetplace - satplace_plans[i, :]
        #=
        if i < 10
            println()
            print("vec_stot : ")
            println(vec_stot)
            print("targetplace : ")
            println(targetplace)
            print("satplace : ")
            println(satplace_plans[i , :])
        end
        =#

        cθ = dot(vec_stot, -satplace_plans[i, :])
        cθ = cθ / (norm(vec_stot) * norm(satplace_plans[i, :]))
        cθ_st = dot(targetplace, satplace_plans[i, :])
        cθ_st = cθ_st / (norm(targetplace) * norm(satplace_plans[i, :]))
        #=
        if i%10 == 0
            print("satplace : ")
            println(ECEFtoGeodetic(satplace_plans[i,:]))
            print("cθ : ")
            println(cθ)
            print("cθ_st : ")
            println(cθ_st)
        end
        =#
        if cθ < cθ_lim        # 撮影範囲外にいる時は除外
            continue
        elseif cθ_st < 0      # 地球の裏側にいる時は除外
            continue
        elseif cθ > cθ_shooting
            cθ_shooting = cθ
            shootingtime = i
            shooting_vec = vec_stot
        end
    end

    shooting_vec = shooting_vec / norm(shooting_vec)

    print("cθ_shooting : ")
    println(cθ_shooting)

    return shootingtime, shooting_vec
end
=#

"""
shootingtime = shootingtime_decision2(satplace_plans, targetplace, currenttime, limit_time, cam_viewangle, sat_axisval)

撮影時刻を決定する関数

# Argments
 - `satplace_plans`：現在時刻~撮影期限までの衛星位置ベクトル@ECEF
 - `targetplace`：撮影対象の位置ベクトル@ECEF
 - `currenttime`：現在時刻（ループ回数）
 - `limit_time`：撮影期限までの時間 （ループ回数　 現在時刻 ~ (currenttime + limit_time)*dt  の間で撮影する）
 - `cam_viewangle`：カメラの視野角 [deg]
 - `sat_axisval`：衛星のポインティング時回転角度限界 [deg]
# Return
- `shootingtime`：撮影時刻(呼び出し元でのループ回数)
"""

function shootingtime_decision2(satplace_plans, targetplace, currenttime, limit_time,cam_viewangle, sat_axisval)
    shootingtime = 0

    shootingtime = 1
    satplace_geod = ECEFtoGeodetic(satplace_plans[1, :])
    targetplace_geod = ECEFtoGeodetic(targetplace)
    mindistance = sqrt((satplace_geod[1]-targetplace[1])^2+(satplace_geod[2]-targetplace[2])^2)

    # 衛星が撮影地点に最も近くなる時間を検索
    for i = 2:limit_time
        satplace_geod = ECEFtoGeodetic(satplace_plans[i, :])
        distance = sqrt((satplace_geod[1]-targetplace_geod[1])^2+(satplace_geod[2]-targetplace_geod[2])^2)
        if distance < mindistance
            shootingtime = i
            mindistance = distance
        end
    end

    #撮影範囲の限界値
    cθ_lim = cosd(cam_viewangle/4 + sat_axisval)

    # 撮影範囲内にいるか判断

    shooting_vec = targetplace - satplace_plans[shootingtime, :]
    shooting_vec = shooting_vec / norm(shooting_vec)
    print("shooting_vec : ")
    println(shooting_vec)
    print("satplace : ")
    println(satplace_plans[shootingtime, :])
    cθ = dot(shooting_vec, -satplace_plans[shootingtime, :])
    cθ = cθ / (norm(shooting_vec) * norm(satplace_plans[shootingtime, :]))

    if cθ < cθ_lim
        shooting_vec = [0., 0., 0.]
        shootingtime = 1
    else
        shootingtime = shootingtime + currenttime -1
    end

    #=
    r = sqrt(shooting_vec[1]^2+shooting_vec[2]^2)
    tφ = r/shooting_vec[3]
    if tφ > tand(sat_axisval)
        shooting_vec[1] = shooting_vec[1] * shooting_vec[3] * tand(sat_axisval) / r
        shooting_vec[2] = shooting_vec[2] * shooting_vec[3] * tand(sat_axisval) / r

        shooting_vec = shooting_vec / norm(shooting_vec)
    end
    =#
    return shootingtime #, shooting_vec
end

"""
targetqua = targetqua_dicision(shooting_vec, current_attqua)

撮影時の目標姿勢クォータニオンを求める

# Argments
- `targetplace`：撮影対象の位置ベクトル@ECEF
- `satplace`：現在時刻~撮影期限までの衛星位置ベクトル@ECEF
- `satvelo` : 衛星の速度ベクトル @ECEF
- `sat_axisval`：衛星のポインティング時回転角度限界 [deg]

# return
 - `targetqua` : 撮影時の目標姿勢クォータニオン

# Conditions
 - `camera_initialdir` = 衛星が基準姿勢時の、カメラ撮影面の単位法線ベクトル @SEOF

"""
function targetqua_dicision(targetplace, satplace, satvelo, sat_axisval)
    camera_initialdir = [0., 0., 1.]
    targetqua = SatelliteToolbox.Quaternion(0., 1., 0., 0.)

    shooting_vec = targetplace - satplace
    shooting_vec = shooting_vec / norm(shooting_vec)
    println("shootingvec:", shooting_vec)
    #=
    r = sqrt(shooting_vec[1]^2 + shooting_vec[2]^2)
    tφ = r/shooting_vec[3]
    if tφ > tand(sat_axisval)
        shooting_vec[1] = shooting_vec[1] * shooting_vec[3] * tand(sat_axisval) / r
        shooting_vec[2] = shooting_vec[2] * shooting_vec[3] * tand(sat_axisval) / r

        shooting_vec = shooting_vec / norm(shooting_vec)
    end
    =#
    shooting_vec_seof = ecef_to_DCM(satplace, satvelo,true) * shooting_vec

    c2θ = dot(camera_initialdir, shooting_vec_seof) / (norm(camera_initialdir)*norm(shooting_vec_seof))
    println("c2θ", c2θ)
    axis_seof = cross(camera_initialdir, shooting_vec_seof)
    axis_seof = axis_seof / norm(axis_seof)
    println("axis_seof", axis_seof)

    cθ = sqrt((1+c2θ)/2)
    # θが回転限界を超えていたらθを回転限界まで抑える
    if cθ < cosd(sat_axisval)
        cθ = cosd(sat_axisval)
    end
    sθ = sqrt(1- cθ^2)

    targetqua = SatelliteToolbox.Quaternion(cθ, axis_seof[1]*sθ, axis_seof[2]*sθ, axis_seof[3]*sθ)

    return targetqua
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eeago_optimizer/optimize/optimize_lp_cone.jl
# Contains the optimize! routines for LP, SOCP, (and in the future MILP and
# MISOCP) type problems. This also includes functions to add variables,
# linear constraints, soc constraints, and unpack solutions.
#############################################################################

#=
LP          -> COPY TO RELAXED SOLVER AND SOLVE
MILP        -> COPY TO RELAXED SOLVER AND SOLVE
SOCP        -> COPY TO RELAXED SOLVER AND SOLVE
MISOCP      -> COPY TO RELAXED SOLVER AND SOLVE
DIFF_CVX    -> COPY TO NLP SOLVER AND SOLVE (POTENTIAL MULTISTART)
DIFF_NCVX   -> APPLY GLOBAL SOLVER (UNLESS USER REQUEST LOCAL SOLVE THEN NLP)
MINCVX      -> APPLY GLOBAL SOLVER (LOCAL SOLVE OPTION FUTURE FEATURE)
=#

function _unpack_final_subsolve!(m::Optimizer, opt::T; adjust_bnd::Bool = true) where T

    # TODO: ASSUMES NO BRIDGING NECESSARY, MAY NEED TO ADD SUPPORT LATER FOR BRIDGES
    m._termination_status_code = MOI.get(opt, MOI.TerminationStatus())
    m._primal_status_code = MOI.get(opt, MOI.PrimalStatus())

    if MOI.get(opt, MOI.ResultCount()) > 0
        variable_indices = [m._input_to_solution_map[vi] for vi in MOI.get(opt, MOI.ListOfVariableIndices())]
        m._solution = MOI.get(opt, MOI.VariablePrimal(), variable_indices)
        for (F, S) in MOI.get(m, MOI.ListOfConstraints())
            if !(F == SV)
                for ci in MOI.get(m, MOI.ListOfConstraintIndices{F,S}())
                    subsolver_ci = m._input_to_solution_map[ci]
                    primal_val = MOI.get(opt, MOI.ConstraintPrimal(), subsolver_ci)
                     _set_cons_primal!(m, ci, primal_val)
                 end
             end
        end
        m._objective_value = MOI.get(opt, MOI.ObjectiveValue())
        m._objective_bound = MOI.get(opt, MOI.ObjectiveBound())
    end

    return nothing
end

#=
ScalarCoefficientChange:
- GLPK, CPLEX, Gurobi supports
- Cbc doesn't
=#
#DONE

function _bridge_optimizer!(v::Val{MILP}, m)
    #MOIB.add_bridge(m, MOIB.Constraint.SOS1toZOBridge{Float64})
    #MOIB.add_bridge(m, MOIB.Constraint.SOS2toZOBridge{Float64})
    return
end

for S in (LP, SOCP, SDP)
    @eval _bridge_optimizer!(v::Val{$S}, m) = nothing
end
for S in (DIFF_CVX, MICVX)
    @eval function _bridge_optimizer!(v::Val{$S}, m)
        MOIB.add_bridge(m, MOIB.Constraint.SOCtoNonConvexQuadBridge{Float64})
        MOIB.add_bridge(m, MOIB.Constraint.RSOCtoNonConvexQuadBridge{Float64})
        return
    end
end

for (T, sub_optimizer) in ((LP, :_mip_optimizer),
                           (MILP, :_mip_optimizer),
                           (SOCP, :_socp_optimizer),
                           (SDP, :_semidefinite_optimizer),
                           (DIFF_CVX, :_nlp_optimizer),
                           (MICVX, :_minlp_optimizer))

    @eval function optimize!(::Val{$T}, m::Optimizer)

        opt = MOI.instantiate(($sub_optimizer)(m), with_bridge_type = Float64)
        #set_config!(m, opt)
        _bridge_optimizer!(Val{$T}(), opt)
        m._input_to_solution_map = MOI.copy_to(opt, _input_model(m); copy_names = false)

        if _verbosity(m) < 5
            MOI.set(opt, MOI.Silent(), true)
        end
        m._parse_time = time() - m._start_time

        MOI.optimize!(opt)
        _unpack_final_subsolve!(m, opt)
    end
end


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

"""
    DensityRatioWeighter(tdata; [options])

Density ratio weights based on empirical distribution of
variables in target data `tdata`.

## Optional parameters

* `variables` - Variables to consider (default to all)
* `estimator` - Density ratio estimator (default to `LSIF()`)
* `optlib`    - Optimization library (default to `default_optlib(estimator)`)

### Notes

Estimators from `DensityRatioEstimation.jl` are supported.
"""
struct DensityRatioWeighter{DΩ<:AbstractData} <: AbstractWeighter
  tdata::DΩ
  vars::Vector{Symbol}
  dre
  optlib
end

function DensityRatioWeighter(tdata::DΩ;
                              variables=nothing,
                              estimator=LSIF(),
                              optlib=default_optlib(estimator)) where {DΩ<:AbstractData}
  validvars = collect(keys(GeoStatsBase.variables(tdata)))
  wvars = variables ≠ nothing ? variables : validvars
  @assert wvars ⊆ validvars "invalid variables ($wvars) for spatial data"
  DensityRatioWeighter{DΩ}(tdata, wvars, estimator, optlib)
end

function weight(sdata::AbstractData, weighter::DensityRatioWeighter)
  # retrieve method parameters
  tdata  = weighter.tdata
  vars   = weighter.vars
  dre    = weighter.dre
  optlib = weighter.optlib

  @assert vars ⊆ keys(variables(sdata)) "invalid variables ($vars) for spatial data"

  # numerator and denominator samples
  Ω_nu = view(tdata, vars)
  Ω_de = view(sdata, vars)

  # TODO: eliminate this explicit conversion after
  # https://github.com/JuliaEarth/GeoStats.jl/projects/2
  x_nu = collect(eachrow(Ω_nu[1:npoints(Ω_nu),vars]))
  x_de = collect(eachrow(Ω_de[1:npoints(Ω_de),vars]))

  # perform denstiy ratio estimation
  ratios = densratio(x_nu, x_de, dre, optlib=optlib)

  SpatialWeights(domain(sdata), ratios)
end


text:
# %%
using Tmatrix

# inputs
rx, rz = 1e-6, 1.3e-6
n_max = 3
k1_r = 1e7;
k1_i = 0.0;
k2_r = 1.5e7;
k2_i = 1e3;

k1 = Complex(k1_r, k1_i)
k2 = Complex(k2_r, k2_i)

BigFloat_precision_list = 2 .^ [6, 7, 8, 9]

println("***************************************************************")
println("***************************************************************")
println("Calculating T-matrix using Complex numbers")
println("***************************************************************")
println("***************************************************************")

println("using Float64 =====================================")
println("No symmetry exploited:")
@time T_no_symmetry_float_64 = calculate_Tmatrix_for_spheroid(
    rx,
    rz,
    n_max,
    k1,
    k2,
    rotationally_symmetric = false,
    symmetric_about_plane_perpendicular_z = false,
);

# %%
println("exploiting symmetry:")
@time T_with_symmetry_float_64 = calculate_Tmatrix_for_spheroid(
    rx,
    rz,
    n_max,
    k1,
    k2,
    rotationally_symmetric = true,
    symmetric_about_plane_perpendicular_z = true,
);

println("using BigFloat =====================================")
println("exploiting symmetry:")
T_with_symmetry_big_float = []
for precision_here in BigFloat_precision_list
    println("using precision $precision_here")
    @time T = calculate_Tmatrix_for_spheroid(
        rx,
        rz,
        n_max,
        k1,
        k2,
        rotationally_symmetric = true,
        symmetric_about_plane_perpendicular_z = true,
        BigFloat_precision = precision_here,
    )
    append!(T_with_symmetry_big_float, T)
    println()
end

"""
n_array = collect(1:10)
elapsed_time = zeros(size(n_array))
for n = 1:length(n_array)
elapsed_time[n] = @elapsed calculate_Tmatrix_for_spheroid(rx, rz, n_max, k1, k2, rotationally_symmetric=true, symmetric_about_plane_perpendicular_z=true, BigFloat_precision=8)
end
"""
println()
println()
println()
println("***************************************************************")
println("***************************************************************")
println("Calculating T-matrix using real numbers")
println("***************************************************************")
println("***************************************************************")

println("using Float64 =====================================")
println("No symmetry exploited:")
@time T_no_symmetry_float_64_SeparateRealImag =
    calculate_Tmatrix_for_spheroid_SeparateRealImag(
        rx,
        rz,
        n_max,
        k1_r,
        k1_i,
        k2_r,
        k2_i;
        rotationally_symmetric = false,
        symmetric_about_plane_perpendicular_z = false,
    );

println("exploiting symmetry:")
@time T_with_symmetry_float_64_SeparateRealImag =
    calculate_Tmatrix_for_spheroid_SeparateRealImag(
        rx,
        rz,
        n_max,
        k1_r,
        k1_i,
        k2_r,
        k2_i;
        rotationally_symmetric = true,
        symmetric_about_plane_perpendicular_z = true,
    );

println("using BigFloat =====================================")
println("exploiting symmetry:")
T_with_symmetry_big_float_SeparateRealImag = []
for precision_here in BigFloat_precision_list
    println("using precision $precision_here")
    @time T = calculate_Tmatrix_for_spheroid_SeparateRealImag(
        rx,
        rz,
        n_max,
        k1_r,
        k1_i,
        k2_r,
        k2_i;
        rotationally_symmetric = true,
        symmetric_about_plane_perpendicular_z = true,
        BigFloat_precision = precision_here,
    )
    append!(T_with_symmetry_big_float_SeparateRealImag, T)
end


text:
# SPDX-License-Identifier: MIT
# Copyright (c) 2020-2021: Pablo Zubieta
# See LICENSE.md at https://github.com/pabloferz/EnhancedSamplingMD.jl

module Backends

using PyCall

# Signal to Revise.jl that this module should be tracked as a package
__revise_mode__ = :eval

# Exports
export ContextWrapper, supported_backends

# Types
"""
Manages access to the backend-dependent simulation context.
"""
struct ContextWrapper
    backend::Module
    context::PyObject
    sampler::PyObject

    function ContextWrapper(context, sampling_method, callback; kwargs...)
        backend = load_backend(context)
        sampler = backend.bind(context, sampling_method, callback; kwargs...)
        return new(backend, context, sampler)
    end
end

# Methods
supported_backends() = (:HOOMD, :OpenMM)

"""    load_backend(context)

Loads the appropriate backend depending on the simulation context.
"""
function load_backend(context)
    module_name = py"type"(context).__module__
    backend_name = if startswith(module_name, "hoomd")
        :HOOMD
    elseif startswith(module_name, r"(simtk\.|)openmm")
        :OpenMM
    else
        :Invalid
    end
    if backend_name === :Invalid
        backends = join(supported_backends(), ", ")
        throw(ArgumentError("Invalid backend, supported options are $backends"))
    end
    path = joinpath(@__DIR__, string(backend, ".jl"))
    return include(path)::Module
end

end  # module Backends


text:
using Test
using MPI
using ClimateMachine.Mesh.Topologies

function main()
    MPI.Init()
    comm = MPI.COMM_WORLD
    crank = MPI.Comm_rank(comm)
    csize = MPI.Comm_size(comm)

    @assert csize == 3

    topology = BrickTopology(
        comm,
        (0:4, 5:9);
        boundary = ((1, 2), (3, 4)),
        periodicity = (false, true),
    )

    elems = topology.elems
    realelems = topology.realelems
    ghostelems = topology.ghostelems
    sendelems = topology.sendelems
    elemtocoord = topology.elemtocoord
    elemtoelem = topology.elemtoelem
    elemtoface = topology.elemtoface
    elemtoordr = topology.elemtoordr
    elemtobndy = topology.elemtobndy
    nabrtorank = topology.nabrtorank
    nabrtorecv = topology.nabrtorecv
    nabrtosend = topology.nabrtosend

    globalelemtoface = [
        1 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2
        1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 2
        4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
        3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
    ]

    globalelemtoordr = ones(Int, size(globalelemtoface))

    globalelemtocoord = Array{Int}(undef, 2, 4, 16)
    globalelemtocoord[:, :, 1] = [0 1 0 1; 5 5 6 6]
    globalelemtocoord[:, :, 2] = [1 2 1 2; 5 5 6 6]
    globalelemtocoord[:, :, 3] = [1 2 1 2; 6 6 7 7]
    globalelemtocoord[:, :, 4] = [0 1 0 1; 6 6 7 7]
    globalelemtocoord[:, :, 5] = [0 1 0 1; 7 7 8 8]
    globalelemtocoord[:, :, 6] = [0 1 0 1; 8 8 9 9]
    globalelemtocoord[:, :, 7] = [1 2 1 2; 8 8 9 9]
    globalelemtocoord[:, :, 8] = [1 2 1 2; 7 7 8 8]
    globalelemtocoord[:, :, 9] = [2 3 2 3; 7 7 8 8]
    globalelemtocoord[:, :, 10] = [2 3 2 3; 8 8 9 9]
    globalelemtocoord[:, :, 11] = [3 4 3 4; 8 8 9 9]
    globalelemtocoord[:, :, 12] = [3 4 3 4; 7 7 8 8]
    globalelemtocoord[:, :, 13] = [3 4 3 4; 6 6 7 7]
    globalelemtocoord[:, :, 14] = [2 3 2 3; 6 6 7 7]
    globalelemtocoord[:, :, 15] = [2 3 2 3; 5 5 6 6]
    globalelemtocoord[:, :, 16] = [3 4 3 4; 5 5 6 6]

    globalelemtobndy = [
        1 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 2 2 2 0 0 2
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    ]

    if crank == 0
        nrealelem = 5
        globalelems = [1, 2, 3, 4, 5, 6, 7, 8, 14, 15]
        elemtoelem_expect = [
            1 1 4 4 5 6 7 8 9 10
            2 10 9 3 8 6 7 8 9 10
            6 7 2 1 4 6 7 8 9 10
            4 3 8 5 6 6 7 8 9 10
        ]
        nabrtorank_expect = [1, 2]
        nabrtorecv_expect = UnitRange{Int}[1:3, 4:5]
        nabrtosend_expect = UnitRange{Int}[1:4, 5:6]
    elseif crank == 1
        nrealelem = 5
        globalelems = [6, 7, 8, 9, 10, 1, 2, 3, 5, 11, 12, 14, 15]
        elemtoelem_expect = [
            1 1 9 3 2 6 7 8 9 10 11 12 13
            2 5 4 11 10 6 7 8 9 10 11 12 13
            9 3 8 12 4 6 7 8 9 10 11 12 13
            6 7 2 5 13 6 7 8 9 10 11 12 13
        ]
        nabrtorank_expect = [0, 2]
        nabrtorecv_expect = UnitRange{Int}[1:4, 5:8]
        nabrtosend_expect = UnitRange{Int}[1:3, 4:5]
    elseif crank == 2
        nrealelem = 6
        globalelems = [11, 12, 13, 14, 15, 16, 2, 3, 9, 10]
        elemtoelem_expect = [
            10 9 4 8 7 5 7 8 9 10
            1 2 3 3 6 6 7 8 9 10
            2 3 6 5 10 1 7 8 9 10
            6 1 2 9 4 3 7 8 9 10
        ]
        nabrtorank_expect = [0, 1]
        nabrtorecv_expect = UnitRange{Int}[1:2, 3:4]
        nabrtosend_expect = UnitRange{Int}[1:2, 3:6]
    end

    @test elems == 1:length(globalelems)
    @test realelems == 1:nrealelem
    @test ghostelems == (nrealelem + 1):length(globalelems)
    @test elemtocoord == globalelemtocoord[:, :, globalelems]
    @test elemtoface[:, realelems] ==
          globalelemtoface[:, globalelems[realelems]]
    @test elemtoelem == elemtoelem_expect
    @test elemtobndy == globalelemtobndy[:, globalelems]
    @test elemtoordr == ones(eltype(elemtoordr), size(elemtoordr))
    @test nabrtorank == nabrtorank_expect
    @test nabrtorecv == nabrtorecv_expect
    @test nabrtosend == nabrtosend_expect

    @test collect(realelems) ==
          sort(union(topology.exteriorelems, topology.interiorelems))
    @test unique(sort(sendelems)) == topology.exteriorelems
    @test length(intersect(topology.exteriorelems, topology.interiorelems)) == 0
end

main()


text:
let
  # One step of dmrg (dmrg itself
  # leads to an error during precompilation)
  # Add `precompile(Tuple{typeof(dmrg),MPO,MPS,Sweeps})`
  # to precompile.jl by hand.
  N = 4
  sites = siteinds("S=1", N)
  ampo = AutoMPO()
  for j in 1:(N - 1)
    ampo .+= ("Sz", j, "Sz", j + 1)
    ampo .+= (0.5, "S+", j, "S-", j + 1)
    ampo .+= (0.5, "S-", j, "S+", j + 1)
  end
  H = MPO(ampo, sites)
  psi0 = randomMPS(sites, 10)
  sweeps = Sweeps(1)
  maxdim!(sweeps, 10)
  cutoff!(sweeps, 1E-11)
  PH = ProjMPO(H)
  which_decomp = nothing
  obs = NoObserver()
  quiet = false
  eigsolve_tol = 1e-14
  eigsolve_krylovdim = 3
  eigsolve_maxiter = 1
  eigsolve_verbosity = 0
  ishermitian = true
  eigsolve_which_eigenvalue = :SR
  psi = copy(psi0)
  N = length(psi)
  ITensors.position!(PH, psi0, 1)
  energy = 0.0
  sw = 1
  for sw in 1:nsweep(sweeps)
    # This loop causes problems for precompilation
    #for (b, ha) in sweepnext(N)
    b, ha = 1, 1
    ITensors.position!(PH, psi, b)
    phi = psi[b] * psi[b + 1]
    vals, vecs = ITensors.KrylovKit.eigsolve(
      PH,
      phi,
      1,
      eigsolve_which_eigenvalue;
      ishermitian=ishermitian,
      tol=eigsolve_tol,
      krylovdim=eigsolve_krylovdim,
      maxiter=eigsolve_maxiter,
    )
    energy, phi = vals[1], vecs[1]
    ortho = ha == 1 ? "left" : "right"
    drho = nothing
    spec = replacebond!(
      psi,
      b,
      phi;
      maxdim=maxdim(sweeps, sw),
      mindim=mindim(sweeps, sw),
      cutoff=cutoff(sweeps, sw),
      eigen_perturbation=drho,
      ortho=ortho,
      which_decomp=which_decomp,
    )
    measure!(
      obs; energy=energy, psi=psi, bond=b, sweep=sw, half_sweep=ha, spec=spec, quiet=quiet
    )
    #end
  end
end


text:
using Symbolics
using Test
using LinearAlgebra

# Transpose and multiply

@syms A[1:2, 1:3] B[1:3, 1:1] C::Matrix{Number} x[1:2] y[1:3] z::Vector{Number}

AA = rand(2,3); BB = rand(3,2); xx = rand(2); yy = rand(3)

adjmul_leaf_nodes = [A, B, C, x, y, z, AA, BB, xx, yy]

function adjmul_rand_leaf()
    a = rand(adjmul_leaf_nodes)
    rand(Bool) ?  a' : a
end

function rand_mul_expr(a=adjmul_rand_leaf(),
                       b=adjmul_rand_leaf())

    try
        a * b
    catch err1
        try
            rand(size(a)...) * rand(size(b)...)
        catch err2
            if typeof(err1) == typeof(err2)
                @test typeof(err1) == typeof(err2)
            else
                @test_skip typeof(err1) == typeof(err2)
            end
        end
    end
    sz = try size(a * b) catch err; nothing end

    if sz !== nothing
        try
            size(a)
            size(b)
            @goto test_size
        catch err
            return # no size known
        end
        @label test_size

        if (a isa Adjoint{<:Any, <:AbstractVector} && ndims(b) == 1) || Symbolics.isdot(a, b)
            if size(a*b) != ()
                println("a * b is wrong:")
                @show a b
                @show typeof(a) typeof(b)
                return @test size(a*b) == ()
            else
                return @test true
            end
        end

        if size(a * b) == size(rand(size(a)...) * rand(size(b)...))
            @test true
        else
            println("a * b is wrong:")
            @show a b
            @show typeof(a) typeof(b)
            @test size(a * b) == size(rand(size(a)...) * rand(size(b)...))
        end
    end
end

@testset "fuzz ' and *" begin
    for i=1:1000
        rand_mul_expr()
    end
end

## Mapreduce
#

AA = rand(2,3); BB = rand(3,2); xx = rand(2); yy = rand(3)
CC = rand(2,3,4)

mapreduce_leaf_nodes = [A, B, C, x, y, z, AA, BB, CC, xx, yy]

function mapreduce_rand_leaf()
    a = rand(adjmul_leaf_nodes)
    rand(Bool) ?  a' : a
end

function rand_map_reduce(x=mapreduce_rand_leaf())
    n = ndims(x)

    reducedims = rand(1:n+1, rand(0:n))

    init=rand([nothing, 3])
    function test_run(x)
        if init==nothing
            sum(rand(size(x)...), dims=reducedims)
        else
            sum(rand(size(x)...), dims=reducedims, init=init)
        end
    end

    @test size(test_run(rand(size(x)...))) == size(test_run(x))
end


text:
s2int(s) = parse(Int, s)

function parse_board(s)
    lines = split(s, "\n")
    grid = map(s2int, hcat(map(split, lines)...))
    [map(Set, eachrow(grid)); map(Set, eachcol(grid))]
end

function slurp_bingo(filename)
    lines = readlines(filename)
    numbers = map(s2int, split(lines[1], ','))
    boards_string = join(lines[3:end], "\n")
    boards = split(boards_string, "\n\n")
    numbers, map(parse_board, boards)
end

function solve(filename)
    numbers, boards = slurp_bingo(filename)
    for number in numbers
        for board in boards
            broadcast(row -> delete!(row, number), board)
            if min(map(length, board)...) == 0
                return reduce(+, union(board...)) * number
            end
        end
    end
end

solve("4eg.txt")
solve("4.txt")

score(board, number) = reduce(+, union(board...)) * number
completed(board) = min(length.(board)...) == 0
remove(board, number) = broadcast(row -> delete!(row, number), board)

function solve2(filename)
    numbers, boards = slurp_bingo(filename)
    for number in numbers
        for board in boards
            remove(board, number)
            if all(completed, boards)
                return score(board, number)
            end
        end
    end
end

solve2("4eg.txt")
solve2("4.txt") == 2980


text:
using Base.Test
using PValueAdjust

pvalues = [0.05, 0.05, 0.05]
newPvalues = padjust(pvalues, Bonferroni)
@test_approx_eq newPvalues pvalues * 3.0

@test all(padjust([1,1,1], Bonferroni) .<= 1)

@test_throws DomainError padjust([-1], Bonferroni)
@test_throws DomainError padjust([2], Bonferroni)


text:
module Antenna

using LinearAlgebra
using Interpolations
using Chain
# pattern file
using CSV
using DataFrames
using Match
using Reexport
using Memoize
using FLoops
using Peaks

include("type.jl")
include("utils.jl")
include("plotting.jl")

export
    cal_pattern,
    anten_read,
    rotate_vec_in_sph,
    rotate_vec_in_cart,
    rotate_pattern,
    directivity,
    radiation_intensity,
    radiated_power,
    SLL_maxgain,
    gain,
    Iₛ

# array function
# array factor
AF(point, θ, ϕ, k) = exp(1im * k * (point[1] * sin(θ)cos(ϕ) + point[2] * sin(θ)sin(ϕ) + point[3] * cos(θ)))
# current source
Iₛ(point, θₜ, ϕₜ, k) = AF(point, θₜ, ϕₜ, k)^-1
# Iₛ(point, θₜ, ϕₜ, k) = 1
# apply rotation
# XXX untest
rotate_vec_in_sph = (θ, ϕ, M) -> begin
    x, y, z = sph2cart(θ, ϕ, 1)
    vec = M * [x, y, z]
    θ1, ϕ1, .. = cart2sph(vec[1], vec[2], vec[3])
    # XXX not a elegen way to aviod lossing information
    # what if there is a rotation, boudary condition is when θ = 0
    # XXXX need a mapping function
    if θ == 0
        ϕ = ϕ
    elseif θ1 == 0
        ϕ = ϕ1
    end
    (θ1, ϕ1)
end
rotate_vec_in_cart = (vec, M) -> M * vec
vec_θ(θ, ϕ) = [cos(θ)cos(ϕ), cos(θ)sin(ϕ), -sin(θ)]
vec_ϕ(θ, ϕ) = [-sin(ϕ), cos(ϕ), 0.0]

rotate_pattern = (pattern::anten_pattern, coord::Matrix{Float64}) -> begin

    # θ_grid::Matrix{Float64} = [θ for θ in θ_default, ϕ in ϕ_default]
    # ϕ_grid::Matrix{Float64} = [ϕ for θ in θ_default, ϕ in ϕ_default]

    Gθ_grid::Matrix{Vector{ComplexF64}} = [
        begin
            θ, ϕ = rotate_vec_in_sph(θ, ϕ, coord')
            (coord * vec_θ(θ, ϕ)) * pattern.θ(θ, ϕ)
        end
        for (θ, ϕ) = zip(θ_grid, ϕ_grid)
    ]

    Gϕ_grid::Matrix{Vector{ComplexF64}} = [
        begin
            θ, ϕ = rotate_vec_in_sph(θ, ϕ, coord')
            (coord * vec_ϕ(θ, ϕ)) * pattern.ϕ(θ, ϕ)
        end
        for (θ, ϕ) = zip(θ_grid, ϕ_grid)
    ]

    vec_θ₁_map_grid::Matrix{Vector{Float64}} = [vec_θ(θ, ϕ) for (θ, ϕ) = zip(θ_grid, ϕ_grid)]
    vec_ϕ₁_map_grid::Matrix{Vector{Float64}} = [vec_ϕ(θ, ϕ) for (θ, ϕ) = zip(θ_grid, ϕ_grid)]

    Gθ = (dot.(Gθ_grid, vec_θ₁_map_grid)) .+
         (dot.(Gϕ_grid, vec_θ₁_map_grid))
    Gϕ = (dot.(Gϕ_grid, vec_ϕ₁_map_grid)) .+
         (dot.(Gθ_grid, vec_ϕ₁_map_grid))

    result = anten_pattern(
        θ=LinearInterpolation((θ_default, ϕ_default), Gθ),
        ϕ=LinearInterpolation((θ_default, ϕ_default), Gϕ)
    )
end

# calculate the global pattern
# para∑ can be write like this
# para∑(Pi(point:p, point:pattern.θ , θ, ϕ, θₜ, ϕₜ, k), (point, θ,ϕ, θₜ, ϕₜ,k))
function cal_pattern(point::Vector{anten_point},  θₜ, ϕₜ; k::Float64=k , θ=θ_default, ϕ=ϕ_default)

    # apply_rotation to the pattern
   @floop for p in point
        p.pattern = rotate_pattern(p.pattern, p.local_coord)
    end

    # calculate result
    result_θ = zeros(ComplexF64, length(θ), length(ϕ))
    result_ϕ = zeros(ComplexF64, length(θ), length(ϕ))
    tmp =  zeros(ComplexF64, length(θ), length(ϕ))

  for  p_i in point
        pattern_θ = p_i.pattern.θ
        pattern_ϕ = p_i.pattern.ϕ
        co = p_i.coeffi * Iₛ(p_i.p, θₜ, ϕₜ, k)
       @inbounds for index_c in  CartesianIndices(tmp)
            θi = θ[index_c[1]]
            ϕi = ϕ[index_c[2]]
            tmp[index_c] = co * AF(p_i.p, θi, ϕi, k)
        end
        result_θ .+= tmp.* pattern_θ.(θ_grid,ϕ_grid)
        result_ϕ .+= tmp.* pattern_ϕ.(θ_grid,ϕ_grid)
    end

    anten_pattern(
        θ=LinearInterpolation((θ, ϕ), result_θ .|>abs),
        ϕ=LinearInterpolation((θ, ϕ), result_ϕ .|>abs),
    )
end

# order: the column number for θ,ϕ,Sθ,Sϕ, for HFSS is [2, 1, 4, 3]
function anten_read(filepath, type="hfss"; unit="abs", factor=1, order=[2, 1, 4, 3])
    #pick data from the dataframe
    unit_corrector = unit -> begin
        if unit == "db"
            x -> 10^(x / 10) |> sqrt
        else
            x -> x * factor
        end
    end
    convert_from_hfss = x -> begin
        (
            θ=deg2rad.(x[:, order[1]]),
            ϕ=deg2rad.(x[:, order[2]]),
            Gθ=x[:, order[3]] .|> unit_corrector(unit),
            Gϕ=x[:, order[4]] .|> unit_corrector(unit)
        )
    end

    # read file
    raw_anten = CSV.File(filepath) |> DataFrame |>
                @match type begin
                    "hfss" => convert_from_hfss
                    _ => convert_from_hfss
                end

    _3Dvec_2_1Dvec = x -> LinRange(minimum(x), maximum(x), size(unique(x), 1))

    translate2Interpolation = (θ, ϕ, G) -> @chain begin
        reshape(G, (size(ϕ, 1), size(θ, 1)))
        permutedims((order[1], order[2]))
        LinearInterpolation((θ, ϕ), _)
    end
    θ = _3Dvec_2_1Dvec(raw_anten.θ)
    ϕ = _3Dvec_2_1Dvec(raw_anten.ϕ)

    anten_pattern(
        θ=translate2Interpolation(θ, ϕ, raw_anten.Gθ),
        ϕ=translate2Interpolation(θ, ϕ, raw_anten.Gϕ)
    )
end

# function below may has performance problem that need to be check because of the usage of lambda inside the comprehension
# unit of radiation intensity is w/sr (watt/unit solid angle) refer to Antenna Theory 2-12a
function radiation_intensity(pattern::anten_pattern)
    (θ, ϕ) -> 1 / 2(120pi) * (abs(pattern.θ(θ, ϕ))^2 + abs(pattern.ϕ(θ, ϕ))^2)
end
# unit of radiated power is w (watt) refer to Antenna Theory 2-13
@memoize function radiated_power(pattern::anten_pattern)
    result = sum([sin(θ) * U * 2pi * pi * 1 / size(ϕ_default, 1)size(θ_default, 1)
                  for (U, θ) in zip(radiation_intensity(pattern).(θ_grid, ϕ_grid), θ_grid)])
end

# no unit refer to Antenna Theory 2-16
@memoize function directivity(pattern::anten_pattern)
    (θ, ϕ) -> 4pi * radiation_intensity(pattern)(θ, ϕ) / radiated_power(pattern)
end
# no unit refer to Antenna Theory 2-16
gain = (pattern, inciden_power) -> function (θ, ϕ)
    4pi * radiation_intensity(pattern)(θ, ϕ) / inciden_power
end

function SLL_maxgain(pattern::anten_pattern;θ=0.0, ϕ=0.0)
    d = directivity(pattern)
    curve = [d((mod_angle_deg(θ,ϕ) .|>deg2rad)...) for θ in -100:100, ϕ in ϕ]
    val = maxima(curve) |> sort!
    val[end]/val[end-1], val[end]
end

end


text:
module Crayons

export Crayon, CrayonStack, merge, @crayon_str

include("crayon.jl")
include("downcasts.jl")
include("crayon_stack.jl")
include("crayon_wrapper.jl")
include("test_prints.jl")
include("logo.jl")
include("consts.jl")
include("macro.jl")

end # module


text:
#=
Compute difference between values at points.

This hides knowledge of the values at points.
E.G. the values are colors.

!!! This is deep in loops, i.e. in the bottleneck, it is important to be fast.

In math terms, a metric.
Here we are not concerned if the metric does not have properties of a true metric,
We want property "non-negative".
We don't care property "symmetric" and other properties of a metric.

More generally, the metric should be a property of the tensor.
=#

#=
Apply a metric between colors at points in images,
returning a large value when points are out of bounds.
Require patchPointTarget to be in bounds,
but not require patchPointCorpus to be in bounds,

The original implemented its own metric and iterated over pixelels.
Here we don't use the word "pixel", which is fraught with ambiguity.

Here use Julia built-in colordiff()

Wolfram says CIE2000 is 3x slower than CIE94.
In Julia DE_2000 versus DE_94

TODO more generally, this should take a comparator function,
instead of specializing to Color.
=#
#=
function pointDifference(
    patchPointTarget::CartesianIndex{DimensionCount},
    patchPointCorpus::CartesianIndex{DimensionCount}
    ) where {DimensionCount}

    result = 0

    # Require both points to be in bounds

    result = colordiff(
        patchPointCorpus.maskedImage.image[ patchPointCorpus.point ],
        patchPointTarget.maskedImage.image[ patchPointTarget.point ],
        metric=DE_94() )
    # @debug "Color diff" result patchPointCorpus.point patchPointTarget.point

    # TODO original algorithm also used weight maps

    return result
end
=#

#=
Specialized: first argument is color.
More generally, call it a "value."
=#
function pointDifference(targetColor, patchPointCorpus )
    result = 0
    # Require both points to be in bounds
    corpusColor = patchPointCorpus.maskedImage.image[ patchPointCorpus.point ]
    result = colordiff(
        corpusColor,
        targetColor,
        metric )
    # WAS metric=DE_94() i.e. computing the metric deep in the loops
    return result
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# # Heights_2 problem with restricted prior on mu.

using DynamicHMCModels

ProjDir = @__DIR__
cd(ProjDir)

# Import the dataset.

data = DataFrame(CSV.read(joinpath("..", "..", "data", "Howell1.csv"), delim=';'));

# Use only adults and standardize

df = filter(row -> row[:age] >= 18, data);

# Flat `σ`, see below.

Base.@kwdef mutable struct Heights_2{Ty <: AbstractVector}
    "Observations."
    y::Ty
end;

# Write a function to return properly dimensioned transformation.

function make_transformation(model::Heights_2)
    as((σ = asℝ₊, μ  = as(Real, 100, 250)), )
end

model = Heights_2(;y = df[:, :height])

# Then make the type callable with the parameters *as a single argument*. Very constraint prior on μ. Flat σ.

function (model::Heights_2)(θ)
    @unpack y = model   # extract the data
    @unpack μ, σ = θ
    loglikelihood(Normal(μ, σ), y) + logpdf(Normal(178, 0.1), μ) +
    logpdf(Uniform(0, 50), σ)
end;

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display
println()

DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics) |> display
println()

a3d = Array{Float64, 3}(undef, 1000, 2, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1, j] = values(posterior[i].μ)
    a3d[i, 2, j] = values(posterior[i].σ)
  end
end

pnames = ["μ", "σ"]
sections =   Dict(
  :parameters => pnames,
)
chns = create_mcmcchains(a3d, pnames, sections, start=1);

# cmdstan result

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
         Mean         SD       Naive SE       MCSE      ESS
sigma  24.604616 0.946911707 0.0149719887 0.0162406632 1000
   mu 177.864069 0.102284043 0.0016172527 0.0013514459 1000

Quantiles:
         2.5%       25.0%     50.0%     75.0%     97.5%
sigma  22.826377  23.942275  24.56935  25.2294  26.528368
   mu 177.665000 177.797000 177.86400 177.9310 178.066000
";

describe(chns)

# end of m4.2d.jl


text:
using DifferentialEquations

function f(du,u,p,t)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end

function g(du,u,p,t)
  du[1] = p[3]*u[1]
  du[2] = p[4]*u[2]
end

p = [1.5,1.0,0.1,0.1];
prob = SDEProblem(f,g,[1.0,1.0],(0.0,10.0),p);

function prob_func(prob,i,repeat)
  x = 0.3rand(2)
  remake(prob,p=[p[1:2];x])
end

ensemble_prob = EnsembleProblem(prob,prob_func=prob_func);

include("../../utils.jl")

@gctime solve(ensemble_prob,SRIW1(),trajectories=1_000_000).u[end].u[end]

#= Results on Julia master with 18 threads at time of PR:
┌─────────┬────────────┬─────────┬──────────────┬───────────────────┬──────────┬────────────┐
│         │ total time │ gc time │ max GC pause │ time to safepoint │ max heap │ percent gc │
│         │         ms │      ms │           ms │                ms │       MB │          % │
├─────────┼────────────┼─────────┼──────────────┼───────────────────┼──────────┼────────────┤
│ minimum │      44143 │   18857 │         9822 │                 0 │    34017 │         43 │
│  median │      44160 │   19057 │         9919 │                 0 │    34118 │         43 │
│ maximum │      44429 │   19152 │         9948 │                 0 │    34202 │         43 │
└─────────┴────────────┴─────────┴──────────────┴───────────────────┴──────────┴────────────┘
=#

max_stars_repo_name: vtjnash/GCBenchmarks

max_issues_repo_name: vtjnash/GCBenchmarks

max_forks_repo_name: vtjnash/GCBenchmarks


text:
# generated
module ko_KR # LOLTools.DataDragon
include(normpath(@__DIR__, "champion.jl"))
include(normpath(@__DIR__, "item.jl"))
include(normpath(@__DIR__, "summoner.jl"))
end # module LOLTools.DataDragon.ko_KR


text:
function animatepopopt(dir::String)
  Plots.scalefontsizes()
  sctrpoplim = (args...;kwargs...) -> scatter(args...;kwargs...,xlims=(1, 100), ylims=(0,1), clims=(4, 9))
  sctroptlim = (args...;kwargs...) -> scatter(args...;kwargs...,xlims=(-6, 0), ylims=(0,1))

  scp = ScatterPop(sctrpoplim, dir);
  sco = ScatterOpt(sctroptlim, dir);
  l = @layout [a b]
  mpl = MultiPlot((ps...) -> plot(ps...; layout=(2,1), size=(600,800)), scp, sco; init=false)

  anim = @animate for i in eachindex(scp.data)
    pl = NaiveGAflux.plotgen(mpl, i)
    pl.subplots[1].attr[:title] = "Generation $i"
  end
  gif(anim, joinpath(dir, "evo.gif"), fps = 2)
end


text:
function build(m::Model)
    m._nlexpressions = []
    jm = Complementarity.MCPModel()

    build_variables!(m, jm)

    build_implicitconstraints!(m, jm)

    build_zeroprofit!(m, jm)

    build_marketclearance!(m, jm)

    build_incomebalance!(m, jm)

    return jm
end


text:
include("Cells.jl")
mutable struct Identifier
  id::Int64
  auth::Int64
  username::Int64
  status::Int64
end
mutable struct JSession
  ID::Identifier
  KeyBindings::Dict
  Extensions::Dict
  username::String
  Cells::Dict
end


text:
int_rules_1_1_1_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.1.1*(a+b*x)^m =#
    @apply_utils Antiderivative(1 / ~x, ~x) => log(~x)
    @apply_utils Antiderivative((~x) ^ ~(m'), ~x) => (~x) ^ (~m + 1) / (~m + 1) <-- FreeQ(~m, ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(1 / (~a + ~(b') * ~x), ~x) => log(RemoveContent(~a + ~b * ~x, ~x)) / ~b <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * ~x) ^ ~m, ~x) => (~a + ~b * ~x) ^ (~m + 1) / (~b * (~m + 1)) <-- FreeQ([~a, ~b, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * ~u) ^ ~m, ~x) => (1 / Coefficient(~u, ~x, 1)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~m, ~x), ~x, ~u) <-- FreeQ([~a, ~b, ~m], ~x) && (LinearQ(~u, ~x) && NeQ(~u, ~x))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.1 (a+b x)^m.jl


text:
module smac_util

using Rayuela
using HDF5

export run_demos_query_base, run_demos_train_query_base

# === Saving functions ===
function save_results_pq_query_base(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, train_error, recall)
  for i = 1:length(C)
    h5write(bpath, "$(trial)/C_$i", C[i])
  end
  h5write(bpath, "$(trial)/B", convert(Matrix{UInt8}, B.-1))
  h5write(bpath, "$(trial)/train_error", train_error)
  h5write(bpath, "$(trial)/recall", recall)
end

function save_results_pq(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, train_error, B_base, recall)
  h5write(bpath, "$(trial)/B_base", convert(Matrix{UInt8}, B_base.-1))
  save_results_pq_query_base(bpath, trial, C, B, train_error, recall)
end

function save_results_opq_query_base(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, R::Matrix{Float32}, train_error, recall)
  h5write(bpath, "$(trial)/R", R)
  save_results_pq_query_base(bpath, trial, C, B, train_error, recall)
end

function save_results_opq(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, R::Matrix{Float32}, train_error, B_base, recall)
  h5write(bpath, "$(trial)/B_base", convert(Matrix{UInt8}, B_base.-1))
  save_results_opq_query_base(bpath, trial, C, B, R, train_error, recall)
end

function save_results_lsq_query_base(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, R::Matrix{Float32}, train_error, opq_error, recall)
  h5write(bpath, "$(trial)/opq_base", opq_error)
  save_results_opq_query_base(bpath, trial, C, B, R, train_error, recall)
end

function save_results_lsq(
  bpath::String, trial::Integer, C::Vector{Matrix{Float32}}, B, R::Matrix{Float32}, train_error, opq_error, B_base, recall)
  h5write(bpath, "$(trial)/B_base", convert(Matrix{UInt8}, B_base.-1))
  save_results_lsq_query_base(bpath, trial, C, B, R, train_error, opq_error, recall)
end

function load_chainq(fname::String, m::Integer, trial::Integer)
  B = h5read(fname, "$trial/B"); B = convert(Matrix{Int16}, B); B.+=1
  R = h5read(fname, "$trial/R")
  chainq_error = h5read(fname, "$trial/train_error")
  C = Vector{Matrix{Float32}}(m)
  for i=1:(m); C[i] = h5read(fname, "$trial/C_$i"); end
  return C, B, R, chainq_error
end

function load_rvq(fname::String, m::Integer, trial::Integer)
  B = h5read(fname, "$trial/B"); B = convert(Matrix{Int16}, B); B.+=1
  rvq_error = h5read(fname, "$trial/train_error")
  C = Vector{Matrix{Float32}}(m)
  for i=1:(m); C[i] = h5read(fname, "$trial/C_$i"); end
  return C, B, rvq_error
end

# === Data loading function ===
function load_experiment_data(
  dataset_name::String,
  ntrain::Integer, nbase::Integer, nquery::Integer, V::Bool=false)

  Xt = Rayuela.read_dataset(dataset_name, ntrain, V)
  Xb = Rayuela.read_dataset(dataset_name * "_base", nbase, V)
  Xq = Rayuela.read_dataset(dataset_name * "_query", nquery, V)[:,1:nquery]
  gt = Rayuela.read_dataset(dataset_name * "_groundtruth", nquery, V)
  if dataset_name == "SIFT1M" || dataset_name == "GIST1M"
    gt = gt .+ 1
  end

  if dataset_name != "Deep1M"
    gt = convert( Vector{UInt32}, gt[1,1:nquery] )
  else
    gt = convert( Vector{UInt32}, gt[1:nquery] )
  end

  return Xt, Xb, Xq, gt
end

# === demos query base ===
function run_demos_query_base(
  dataset_name="labelme",
  m::Integer=8,
  h::Integer=256,
  niter::Integer=5,
  sr_method::String="SR_D",
  ilsiter::Integer=8,
  icmiter::Integer=4,
  randord::Bool=true,
  npert::Integer=4,
  schedule::Integer=1,
  p::Float64=0.5)

  if !(sr_method in ["LSQ", "SR_D", "SR_C"])
    error("Unknown sr_method $sr_method")
  end

  nquery, nbase, knn = 0, 0, 0
  if dataset_name == "MNIST"
    ntrain, nquery, nbase, knn = Int(60e3), Int(10e3), Int(60e3), Int(1e3)
  elseif dataset_name == "labelme"
    ntrain, nquery, nbase, knn = Int(20e3), Int(2e3), Int(20019), Int(1e3)
    # ntrain, nquery, nbase, knn = Int(20019), Int(2e3), Int(20019), Int(1e3)
  else
    error("dataset unknown")
  end
  # @show ntrain, nquery, nbase, knn

  verbose = true

  Xt, Xb, Xq, gt = load_experiment_data(dataset_name, ntrain, nbase, nquery, verbose)
  @show size(Xt), size(Xb), size(Xq), size(gt)
  d, _ = size(Xt)

  nsplits_train = 1

  recall = [0]
  # trial = rand(1:10, 1)[1]
  for trial = 1:10
    if sr_method == "LSQ"
      C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
      C, B, R, train_error, recall = Rayuela.experiment_lsq_cuda_query_base(Xt, B, C, R, Xq, gt, m-1, h, niter, ilsiter,
          icmiter, randord, npert, knn, nsplits_train, verbose)
      save_results_lsq_query_base("./results/smac/$(lowercase(dataset_name))/lsq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
    elseif sr_method == "SR_D"
      C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
      C, B, R, train_error, recall = Rayuela.experiment_sr_cuda_query_base(Xt, B, C, R, Xq, gt, m-1, h, niter, ilsiter,
          icmiter, randord, npert, knn, nsplits_train, sr_method, schedule, p, verbose)
        save_results_lsq_query_base("./results/smac/$(lowercase(dataset_name))/srd_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
    elseif sr_method == "SR_C"
      C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
      C, B, R, train_error, recall = Rayuela.experiment_sr_cuda_query_base(Xt, B, C, R, Xq, gt, m-1, h, niter, ilsiter,
          icmiter, randord, npert, knn, nsplits_train, sr_method, schedule, p, verbose)
      save_results_lsq_query_base("./results/smac/$(lowercase(dataset_name))/src_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, recall)
    end
  end

end

function run_demos_train_query_base(
  dataset_name="SIFT1M",
  m::Integer=8,
  h::Integer=256,
  niter::Integer=5,
  sr_method::String="SR_D",
  ilsiter::Integer=8,
  icmiter::Integer=4,
  randord::Bool=true,
  npert::Integer=4,
  schedule::Integer=1,
  p::Float64=0.5)

  if !(sr_method in ["LSQ", "SR_D", "SR_C"])
    error("Unknown sr_method $sr_method")
  end

  nquery, nbase, knn = 0, 0, 0
  if dataset_name == "SIFT1M" || dataset_name == "Deep1M" || dataset_name == "Convnet1M"
    ntrain, nquery, nbase, knn = Int(1e5), Int(1e4), Int(1e6), Int(1e3)
  else
    error("dataset unknown")
  end
  # @show ntrain, nquery, nbase, knn

  verbose = true

  Xt, Xb, Xq, gt = load_experiment_data(dataset_name, ntrain, nbase, nquery, verbose)
  d, _ = size(Xt)

  nsplits_train = 1
  nsplits_base = 4

  recall = [0]
  # trial = rand(1:10, 1)[1]
  for trial = 1:10
    if sr_method == "LSQ"
      C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
      C, B, R, train_error, B_base, recall = Rayuela.experiment_lsq_cuda(Xt, B, C, R, Xb, Xq, gt, m-1, h, niter, ilsiter,
          icmiter, randord, npert, knn, nsplits_train, nsplits_base, verbose)
      save_results_lsq("./results/smac/$(lowercase(dataset_name))/lsq_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, B_base, recall)

    elseif sr_method == "SR_D"
      C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
      C, B, R, train_error, B_base, recall = Rayuela.experiment_sr_cuda(Xt, B, C, R, Xb, Xq, gt, m-1, h, niter, ilsiter,
          icmiter, randord, npert, knn, nsplits_train, nsplits_base, sr_method, schedule, p, verbose)
      save_results_lsq("./results/smac/$(lowercase(dataset_name))/srd_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, B_base, recall)

    elseif sr_method == "SR_C"
      C, B, R, chainq_error = load_chainq("./results/$(lowercase(dataset_name))/chainq_m$(m-1)_it$(niter).h5", m-1, trial)
      C, B, R, train_error, B_base, recall = Rayuela.experiment_sr_cuda(Xt, B, C, R, Xb, Xq, gt, m-1, h, niter, ilsiter,
          icmiter, randord, npert, knn, nsplits_train, nsplits_base, sr_method, schedule, p, verbose)
      save_results_lsq("./results/smac/$(lowercase(dataset_name))/src_m$(m-1)_it$(niter).h5", trial, C, B, R, train_error, chainq_error, B_base, recall)
    end
  end

end

# Query/base datasets
for niter = [100]
  # run_demos_query_base("labelme", 8,  256, niter, "SR_D", 9, 3, true,  1, 1, 0.43098784299895454)
  # # run_demos_query_base("labelme", 16, 256, niter, "SR_D", 8, 4, true,  4, 1, 0.5)  # No change here
  # run_demos_query_base("MNIST",   8,  256, niter, "SR_D", 9, 3, false, 5, 1, 0.18979255389609623)
  # run_demos_query_base("MNIST",   16, 256, niter, "SR_D", 8, 4, false, 4, 1, 0.8282107865533627)
end

# Train/Query/base datasets
# TODO change the function so it savez the results
for niter = [25, 100]
  run_demos_train_query_base("SIFT1M",    8,  256, niter, "SR_D", 8,  4, true,  4, 1, 0.6458745069743886)
  run_demos_train_query_base("SIFT1M",    16, 256, niter, "SR_D", 7,  4, true,  2, 1, 0.18722222602931293)

  # run_demos_train_query_base("Deep1M",    8,  256, niter, "SR_D", 8,  4, true,  4, 1, 0.5) # No change here
  run_demos_train_query_base("Deep1M",    16, 256, niter, "SR_C", 15, 2, true,  2, 1, 0.9534092523209057)

  run_demos_train_query_base("Convnet1M", 8,  256, niter, "SR_C", 8,  4, true,  4, 1, 0.7134116312190524)
  run_demos_train_query_base("Convnet1M", 16, 256, niter, "SR_C", 10, 3, false, 5, 1, 0.937363908221641)
end

end # module smac_util


text:
import Base.CoreLogging.LogState

"""
    LoggingStream(logger; level, id)

An `IO` object which collects incoming calls to `write` and writes them to the
Julia logging system via `logger`. Most useful when combined with
`LineBufferedIO`. The standard logging `_id` field will be set to `id`.
"""
struct LoggingStream{Logger <: AbstractLogger, Level} <: IO
    logger::Logger
    logstate::LogState
    id::Symbol
    level::Level
end

function LoggingStream(logger; id, level=Logging.Info)
    LoggingStream(logger, LogState(logger), Symbol(id), level)
end

function Base.unsafe_write(s::LoggingStream, p::Ptr{UInt8}, n::UInt)
    # Remove any trailing newline
    level = s.level
    # Unclear what to set this metadata to, it's not really available.
    id = s.id
    file = "unknown"
    line = 0
    group = :io
    _module = Logging
    # These tests capture the essence of the logging early-bailout in Base.CoreLogging.
    #
    # NOTE: These refer to internals of CoreLogging and will need to be kept in
    # sync if things change upstream.
    if #==# level >= Base.CoreLogging._min_enabled_level[] &&
            level >= s.logstate.min_enabled_level &&
            Logging.shouldlog(s.logger, level, _module, group, id)
        m = (n > 0 && unsafe_load(p, n) == UInt8('\n')) ? n-1 : n
        # Wrap in Text to force plain rather than markdown text
        message = Text(unsafe_string(p, m))
        Logging.handle_message(s.logger, level, message, _module, group, id, file, line)
    end
    # Could support some approximation of dynamic scoping with the following,
    # but it's a bit of a hack.
    # @logmsg level message _id=id _file=file _line=line _group=group _module=_module
    return n
end

function Base.write(s::LoggingStream, c::UInt8)
    Base.unsafe_write(s, Ref(c), 1)
end


text:
using FaceDatasets
using Base.Test

# Cootes dataset is bundled with a package, so it's safe to load it
load_images(CootesDataset)
load_shapes(CootesDataset)