text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct UpdateDataLakeStoreAccountProperties <: SwaggerModel
    defaultGroup::Any # spec type: Union{ Nothing, String } # spec name: defaultGroup
    encryptionConfig::Any # spec type: Union{ Nothing, UpdateEncryptionConfig } # spec name: encryptionConfig
    firewallRules::Any # spec type: Union{ Nothing, Vector{UpdateFirewallRuleWithAccountParameters} } # spec name: firewallRules
    virtualNetworkRules::Any # spec type: Union{ Nothing, Vector{UpdateVirtualNetworkRuleWithAccountParameters} } # spec name: virtualNetworkRules
    firewallState::Any # spec type: Union{ Nothing, String } # spec name: firewallState
    firewallAllowAzureIps::Any # spec type: Union{ Nothing, String } # spec name: firewallAllowAzureIps
    trustedIdProviders::Any # spec type: Union{ Nothing, Vector{UpdateTrustedIdProviderWithAccountParameters} } # spec name: trustedIdProviders
    trustedIdProviderState::Any # spec type: Union{ Nothing, String } # spec name: trustedIdProviderState
    newTier::Any # spec type: Union{ Nothing, String } # spec name: newTier

    function UpdateDataLakeStoreAccountProperties(;defaultGroup=nothing, encryptionConfig=nothing, firewallRules=nothing, virtualNetworkRules=nothing, firewallState=nothing, firewallAllowAzureIps=nothing, trustedIdProviders=nothing, trustedIdProviderState=nothing, newTier=nothing)
        o = new()
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("defaultGroup"), defaultGroup)
        setfield!(o, Symbol("defaultGroup"), defaultGroup)
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("encryptionConfig"), encryptionConfig)
        setfield!(o, Symbol("encryptionConfig"), encryptionConfig)
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("firewallRules"), firewallRules)
        setfield!(o, Symbol("firewallRules"), firewallRules)
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("virtualNetworkRules"), virtualNetworkRules)
        setfield!(o, Symbol("virtualNetworkRules"), virtualNetworkRules)
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("firewallState"), firewallState)
        setfield!(o, Symbol("firewallState"), firewallState)
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("firewallAllowAzureIps"), firewallAllowAzureIps)
        setfield!(o, Symbol("firewallAllowAzureIps"), firewallAllowAzureIps)
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("trustedIdProviders"), trustedIdProviders)
        setfield!(o, Symbol("trustedIdProviders"), trustedIdProviders)
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("trustedIdProviderState"), trustedIdProviderState)
        setfield!(o, Symbol("trustedIdProviderState"), trustedIdProviderState)
        validate_property(UpdateDataLakeStoreAccountProperties, Symbol("newTier"), newTier)
        setfield!(o, Symbol("newTier"), newTier)
        o
    end
end # type UpdateDataLakeStoreAccountProperties

const _property_map_UpdateDataLakeStoreAccountProperties = Dict{Symbol,Symbol}(Symbol("defaultGroup")=>Symbol("defaultGroup"), Symbol("encryptionConfig")=>Symbol("encryptionConfig"), Symbol("firewallRules")=>Symbol("firewallRules"), Symbol("virtualNetworkRules")=>Symbol("virtualNetworkRules"), Symbol("firewallState")=>Symbol("firewallState"), Symbol("firewallAllowAzureIps")=>Symbol("firewallAllowAzureIps"), Symbol("trustedIdProviders")=>Symbol("trustedIdProviders"), Symbol("trustedIdProviderState")=>Symbol("trustedIdProviderState"), Symbol("newTier")=>Symbol("newTier"))
const _property_types_UpdateDataLakeStoreAccountProperties = Dict{Symbol,String}(Symbol("defaultGroup")=>"String", Symbol("encryptionConfig")=>"UpdateEncryptionConfig", Symbol("firewallRules")=>"Vector{UpdateFirewallRuleWithAccountParameters}", Symbol("virtualNetworkRules")=>"Vector{UpdateVirtualNetworkRuleWithAccountParameters}", Symbol("firewallState")=>"String", Symbol("firewallAllowAzureIps")=>"String", Symbol("trustedIdProviders")=>"Vector{UpdateTrustedIdProviderWithAccountParameters}", Symbol("trustedIdProviderState")=>"String", Symbol("newTier")=>"String")
Base.propertynames(::Type{ UpdateDataLakeStoreAccountProperties }) = collect(keys(_property_map_UpdateDataLakeStoreAccountProperties))
Swagger.property_type(::Type{ UpdateDataLakeStoreAccountProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpdateDataLakeStoreAccountProperties[name]))}
Swagger.field_name(::Type{ UpdateDataLakeStoreAccountProperties }, property_name::Symbol) =  _property_map_UpdateDataLakeStoreAccountProperties[property_name]

const _allowed_UpdateDataLakeStoreAccountProperties_firewallState = ["Enabled", "Disabled"]

const _allowed_UpdateDataLakeStoreAccountProperties_firewallAllowAzureIps = ["Enabled", "Disabled"]

const _allowed_UpdateDataLakeStoreAccountProperties_trustedIdProviderState = ["Enabled", "Disabled"]

const _allowed_UpdateDataLakeStoreAccountProperties_newTier = ["Consumption", "Commitment_1TB", "Commitment_10TB", "Commitment_100TB", "Commitment_500TB", "Commitment_1PB", "Commitment_5PB"]

function check_required(o::UpdateDataLakeStoreAccountProperties)
    true
end

function validate_property(::Type{ UpdateDataLakeStoreAccountProperties }, name::Symbol, val)
    if name === Symbol("firewallState")
        Swagger.validate_param(name, "UpdateDataLakeStoreAccountProperties", :enum, val, _allowed_UpdateDataLakeStoreAccountProperties_firewallState)
    end
    if name === Symbol("firewallAllowAzureIps")
        Swagger.validate_param(name, "UpdateDataLakeStoreAccountProperties", :enum, val, _allowed_UpdateDataLakeStoreAccountProperties_firewallAllowAzureIps)
    end
    if name === Symbol("trustedIdProviderState")
        Swagger.validate_param(name, "UpdateDataLakeStoreAccountProperties", :enum, val, _allowed_UpdateDataLakeStoreAccountProperties_trustedIdProviderState)
    end
    if name === Symbol("newTier")
        Swagger.validate_param(name, "UpdateDataLakeStoreAccountProperties", :enum, val, _allowed_UpdateDataLakeStoreAccountProperties_newTier)
    end
end


text:
using Distributed, SharedArrays
addprocs(3)
@everywhere using KitBase

begin
    vars = Dict{Symbol,Any}()
    vars[:matter] = "gas"
    vars[:case] = "sod"
    vars[:space] = "1d0f0v"
    vars[:flux] = "kfvs"
    vars[:collision] = "bgk"
    vars[:nSpecies] = 1
    vars[:interpOrder] = 1
    vars[:limiter] = "vanleer"
    vars[:boundary] = "fix"
    vars[:cfl] = 0.5
    vars[:maxTime] = 0.2
    vars[:x0] = 0.0
    vars[:x1] = 1.0
    vars[:nx] = 2000
    vars[:pMeshType] = "uniform"
    vars[:nxg] = 0
    vars[:knudsen] = 0.001
    vars[:mach] = 0.0
    vars[:prandtl] = 1.0
    vars[:inK] = 0.0
    vars[:omega] = 0.81
    vars[:alphaRef] = 1.0
    vars[:omegaRef] = 0.5
end

set = KitBase.set_setup(vars)
pSpace = KitBase.set_geometry(vars)
vSpace = KitBase.set_velocity(vars)
gas = KitBase.set_property(vars)
ib = KitBase.set_ib(vars, set, vSpace, gas)
folder = @__DIR__
ks = KitBase.SolverSet(set, pSpace, vSpace, gas, ib, folder)

dt = ks.pSpace.dx[1] / (5.0 + KitBase.sound_speed(ks.ib.primL, ks.gas.γ))
nt = floor(ks.set.maxTime / dt) |> Int

#--- parallel ---#
wp = SharedArray{Float64}((ks.pSpace.nx, 3), init = A -> (A = zeros(ks.pSpace.nx, 3)))
for i = 1:ks.pSpace.nx
    if i <= ks.pSpace.nx ÷ 2
        wp[i, :] .= ks.ib.wL
    else
        wp[i, :] .= ks.ib.wR
    end
end
fwp = SharedArray{Float64}(
    (ks.pSpace.nx + 1, 3),
    init = A -> (A = zeros(ks.pSpace.nx + 1, 3)),
)

@time for iter = 1:nt÷3
    @sync @distributed for i = 2:ks.pSpace.nx
        flux = @view fwp[i, :]
        KitBase.flux_gks!(
            flux,
            wp[i-1, :],
            wp[i, :],
            ks.gas.γ,
            ks.gas.K,
            ks.gas.μᵣ,
            ks.gas.ω,
            dt,
            0.5 * ks.pSpace.dx[i-1],
            0.5 * ks.pSpace.dx[i],
        )
    end

    @sync @distributed for i = 2:ks.pSpace.nx-1
        for j = 1:3
            wp[i, j] += (fwp[i, j] - fwp[i+1, j]) / ks.pSpace.dx[i]
        end
    end
end
"""~13.620491 seconds (2.26 M allocations: 101.219 MiB, 0.22% gc time)"""

#--- serial ---#
w = zeros(ks.pSpace.nx, 3)
for i = 1:ks.pSpace.nx
    if i <= ks.pSpace.nx ÷ 2
        w[i, :] .= ks.ib.wL
    else
        w[i, :] .= ks.ib.wR
    end
end
fw = zeros(ks.pSpace.nx + 1, 3)

@time for iter = 1:nt÷3
    for i = 2:ks.pSpace.nx
        flux = @view fw[i, :]
        KitBase.flux_gks!(
            flux,
            w[i-1, :],
            w[i, :],
            ks.gas.γ,
            ks.gas.K,
            ks.gas.μᵣ,
            ks.gas.ω,
            dt,
            0.5 * ks.pSpace.dx[i-1],
            0.5 * ks.pSpace.dx[i],
        )
    end

    for i = 2:ks.pSpace.nx-1
        for j = 1:3
            w[i, j] += (fw[i, j] - fw[i+1, j]) / ks.pSpace.dx[i]
        end
    end
end
"""~20.830331 seconds (323.96 M allocations: 24.472 GiB, 16.89% gc time)"""

ext: jl


text:
struct Landmark <: AbstractLandmark
    pos_::Vector{Float64}
    id::Int64
    function Landmark(pos::Vector{Float64}, id::Int64)
        new(pos, id)
    end
end

#=
function Base.copy(lm::Landmark)
    return Landmark(copy(lm.pos_), lm.id)
end
=#

function draw(mark::Landmark, p::Plot{T}) where {T}
    p = scatter!(
        [mark.pos_[1]],
        [mark.pos_[2]],
        markershape = :star,
        markersize = 10,
        color = "orange",
    )
    p = annotate!(mark.pos_[1] + 0.5, mark.pos_[2] + 0.5, text("id: $(mark.id)", 10))
end

mutable struct Map <: AbstractObject
    landmarks_::Vector{AbstractLandmark}
    function Map()
        new(Vector{AbstractLandmark}[])
    end
end

function Base.copy(src::Map)
    dst = Map()
    push!(dst, copy(src.landmarks_))
    return dst
end

function Base.push!(map::Map, landmark::AbstractLandmark)
    push!(map.landmarks_, landmark)
end

function Base.push!(map::Map, landmarks::Vector{T}) where {T<:AbstractLandmark}
    for landmark in landmarks
        push!(map.landmarks_, landmark)
    end
end

function draw(map::Map, p::Plot{T}) where {T}
    for landmark in map.landmarks_
        draw(landmark, p)
    end
end

function Base.getindex(map::Map, index::Int)
    return map.landmarks_[index]
end

mutable struct World <: AbstractWorld
    objects_::Vector{AbstractObject}
    xlim_::Vector{Float64}
    ylim_::Vector{Float64}
    debug_::Bool
    function World(xlim::Vector{Float64}, ylim::Vector{Float64}, debug = false)
        new(Vector{AbstractObject}[], [xlim[1], xlim[2]], [ylim[1], ylim[2]], debug)
    end
end

function Base.push!(world::World, obj::AbstractObject)
    push!(world.objects_, obj)
end

function draw(world::AbstractWorld, annota::String)
    p = plot(aspect_ratio = :equal, xlim = world.xlim_, ylim = world.ylim_)
    xpos = world.xlim_[1] + (world.xlim_[2] - world.xlim_[1]) * 0.2
    ypos = world.ylim_[2] - 0.1 * (world.ylim_[2] - world.ylim_[1])
    if annota != nothing
        p = annotate!(xpos, ypos, text(annota, 10))
    end
    for obj in world.objects_
        draw(obj, p)
    end
    return p
end

function draw(world::AbstractWorld, annota::Nothing)
    p = plot(aspect_ratio = :equal, xlim = world.xlim_, ylim = world.ylim_)
    xpos = (world.xlim_[1] + world.xlim_[2]) / 2.0
    ypos = world.ylim_[2] - 0.1 * (world.ylim_[2] - world.ylim_[1])
    for obj in world.objects_
        draw(obj, p)
    end
    return p
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
@doc raw"""
    triag_rng(a, b, m, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Triangular(`a`, `b`, `m`) distribution. Optionally you can set a specific seed.

# Notes

The Triangular distribution is given by:

```math
f(x, a, b, m) = \begin{cases}
\frac{2(x-a)}{(m-a)(b-a)}  & \text{if  } a < x ≤ m, \\
\frac{2(b-x)}{(b-m)(b-a)}  & \text{if  } m < x ≤ c, \\
0 & \text{otherwise}
\end{cases}
```

# Examples

```julia-repl
julia> triag_rng()
1-element Vector{Float64}:
 0.3559088458688944

julia> triag_rng(0,7,2,5)
5-element Vector{Float64}:
 1.3758072115332673
 6.049452463477447
 6.042781317411027
 2.914959243260448
 4.454707036522528
```

# References

W. Stein and M. Keblis. A new method to simulate the triangular distribution. Mathematical and Computer Modelling, Volume 49, Issues 5–6, 2009.
"""
function triag_rng(a=0,b=1,m=.5,shape=1)
    u = get_std_uniform(shape)
    v = get_std_uniform(shape)
    t_min = min.(u,v)
    t_max = max.(u,v)
    c = (m-a)/(b-a)
    t = ((1-c) .* t_min) .+ (c .* t_max)
    X = a .+ (b-a) .*t
    return X
end


text:
# ==============================================================================
#                        size of model
# ==============================================================================
nz = 101; nx = 641; nstep=4000;
# grid cell size
 deltaZ = 10.0        ;  deltaX = 10.0;
hdeltaZ = deltaZ / 2.0; hdeltaX = deltaX / 2.0;
# time step size
f0  = 7.0; dt = 0.001;# source parameter
# ==============================================================================
#                        elastic model parameters
# ==============================================================================
vp = 3300.0; vs = vp / 1.732; density = 2800.0;
courantNumber = vp * dt * sqrt(1.0/deltaZ^2 + 1.0/deltaX^2)
if courantNumber > 1.0
   println("Courant number is: $courantNumber")
   error("time step is too large, simulation will be unstable")
end

# arrays for elastic model parameter
lambda = zeros(nz+2, nx+2);
mu     = zeros(nz+2, nx+2);
rho    = zeros(nz+2, nx+2);

for j = 2 : nx+1
    for i = 2 : nz+1
           rho[i,j] = density
            mu[i,j] = density * vs * vs
        lambda[i,j] = density * (vp*vp - 2.0*vs*vs)
    end
end

# half grid elastic model parameters
hLambdaZ    = 0.0
hLambda2MuZ = 0.0
hMuZ        = 0.0
hMuX        = 0.0
hRhoZX      = 0.0

# arrays for wave field variables
vz = zeros(nz+2, nx+2);
vx = zeros(nz+2, nx+2);
sigmazz = zeros(nz+2, nx+2);
sigmaxx = zeros(nz+2, nx+2);
sigmazx = zeros(nz+2, nx+2);

# temperal variables to save spatial derivative
vdvzdz = 0.0
vdvxdx = 0.0
vdvxdz = 0.0
vdvzdx = 0.0
vdsigmazzdz = 0.0
vdsigmazxdx = 0.0
vdsigmaxxdx = 0.0
vdsigmazxdz = 0.0

# memory variables for CPML, it can be improved to save memory !!!
mdvzdz = zeros(nz+2, nx+2)
mdvxdx = zeros(nz+2, nx+2)
mdvxdz = zeros(nz+2, nx+2)
mdvzdx = zeros(nz+2, nx+2)
mdsigmazzdz = zeros(nz+2, nx+2)
mdsigmazxdx = zeros(nz+2, nx+2)
mdsigmaxxdx = zeros(nz+2, nx+2)
mdsigmazxdz = zeros(nz+2, nx+2)

# ==============================================================================
#                      prepare for CPML damping profile
# ==============================================================================
# thickness of CPML in number of grid points
nPointPml = 10;
thickPmlZ = nPointPml * deltaZ;
thickPmlX = nPointPml * deltaX;

# parameter for CPML profile
nPower   = 2.0
kmax     = 1.0
alphaMax = pi * f0

# damping profile, need to be improved  !!!
# (k, d, alpha, a, b) at grid point, # (hk, hd, halpha, ha, hb) at half grid point along Z direction
 kz=ones(nz);  dz=zeros(nz);  alphaz=zeros(nz);  az=zeros(nz);  bz=zeros(nz);
hkz=ones(nz); hdz=zeros(nz); halphaz=zeros(nz); haz=zeros(nz); hbz=zeros(nz);

# along X direction
 kx=ones(nx);  dx=zeros(nx);  alphax=zeros(nx);  ax=zeros(nx);  bx=zeros(nx);
hkx=ones(nx); hdx=zeros(nx); halphax=zeros(nx); hax=zeros(nx); hbx=zeros(nx);

# specify the obsorbing boundary, origin specify the computation domain, z
reflectCoef=1e-3;
d0Z = -(nPower + 1) * vp * log(reflectCoef) / (2.0 * thickPmlZ);
d0X = -(nPower + 1) * vp * log(reflectCoef) / (2.0 * thickPmlX);
locationZ=0.0; locationX=0.0; abscissa=0.0; normalizedAbscissa=0.0;
originTop = thickPmlZ; originBottom=(nz-1)*deltaZ - thickPmlZ;

# vertical damping profile
for i = 1 : nz
    locationZ = deltaZ * (i-1)

    # upper profile at grid point
    abscissa  = originTop - locationZ
    if abscissa > 0.0
       normalizedAbscissa = abscissa / thickPmlZ
       dz[i] = d0Z * normalizedAbscissa^nPower
       kz[i] = 1.0
       alphaz[i] = alphaMax * (1.0 - normalizedAbscissa) + 0.1 * alphaMax
    end

    # upper profile at half grid point
    abscissa = originTop - (locationZ + hdeltaZ)
    if abscissa > 0.0
       normalizedAbscissa = abscissa / thickPmlZ
       hdz[i] = d0Z * normalizedAbscissa^nPower
       hkz[i] = 1.0
       halphaz[i] = alphaMax * (1.0 - normalizedAbscissa) + 0.1 * alphaMax
    end

    # bottom profile at grid point
    abscissa = locationZ - originBottom
    if abscissa > 0.0
       normalizedAbscissa = abscissa / thickPmlZ
       dz[i] = d0Z * normalizedAbscissa^nPower
       kz[i] = 1.0
       alphaz[i] = alphaMax * (1.0 - normalizedAbscissa) + 0.1 * alphaMax
    end

    # bottom profile at half grid point
    abscissa = (locationZ + hdeltaZ) - originBottom
    if abscissa > 0.0
       normalizedAbscissa = abscissa / thickPmlZ
       hdz[i] = d0Z * normalizedAbscissa^nPower
       hkz[i] = 1.0
       halphaz[i] = alphaMax * (1.0 - normalizedAbscissa) + 0.1 * alphaMax
    end

    bz[i] = exp(- ( dz[i]/ kz[i] +  alphaz[i]) * dt)
   hbz[i] = exp(- (hdz[i]/hkz[i] + halphaz[i]) * dt)

    if (abs( dz[i]) > eps(Float32))
        az[i] =  dz[i] / ( kz[i]*( dz[i]+ kz[i]* alphaz[i])) * ( bz[i]-1.0)
    end
    if (abs(hdz[i]) > eps(Float32))
       haz[i] = hdz[i] / (hkz[i]*(hdz[i]+hkz[i]*halphaz[i])) * (hbz[i]-1.0)
    end
end

# damping in x direction
thickPmlX = nPointPml * deltaX; originLeft= thickPmlX; originRight =(nx-1)*deltaX - thickPmlX;
for j = 1 : nx
    locationX = deltaX * (j-1)
    # left damping profile at grid point
    abscissa = originLeft - locationX
    if (abscissa > 0.0)
       normalizedAbscissa = abscissa / thickPmlX
       dx[j] = d0X * normalizedAbscissa^nPower
       kx[j] = 1.0
       alphax[j] = alphaMax * (1.0 - normalizedAbscissa) + 0.1 * alphaMax
    end

    # left profile at half grid point
    abscissa = originLeft - (locationX + hdeltaX)
    if abscissa > 0.0
       normalizedAbscissa = abscissa / thickPmlX
       hdx[j] = d0X * normalizedAbscissa^nPower
       hkx[j] = 1.0
       halphax[j] = alphaMax * (1.0 - normalizedAbscissa) + 0.1 * alphaMax
    end

    # right damping profile at grid point
    abscissa = locationX - originRight
    if abscissa > 0.0
       normalizedAbscissa = abscissa / thickPmlX
       dx[j] = d0X * normalizedAbscissa^nPower
       kx[j] = 1.0
       alphax[j] = alphaMax * (1.0 - normalizedAbscissa) + 0.1 * alphaMax
    end
    # bottom profile at half grid point
    abscissa = (locationX + hdeltaX) - originRight
    if abscissa > 0.0
       normalizedAbscissa = abscissa / thickPmlX
       hdx[j] = d0Z * normalizedAbscissa^nPower
       hkx[j] = 1.0
       halphax[j] = alphaMax * (1.0 - normalizedAbscissa) + 0.1 * alphaMax
    end

    bx[j] = exp(- ( dx[j]/ kx[j] +  alphax[j]) * dt)
   hbx[j] = exp(- (hdx[j]/hkx[j] + halphax[j]) * dt)

    if (abs( dx[j]) > eps(Float32))
        ax[j] =  dx[j] / ( kx[j]*( dx[j]+ kx[j]* alphax[j])) * ( bx[j]-1.0)
    end
    if (abs(hdx[j]) > eps(Float32))
       hax[j] = hdx[j] / (hkx[j]*(hdx[j]+hkx[j]*halphax[j])) * (hbx[j]-1.0)
    end

end

isz = 20; isx = floor(Int64, nx/2)+1;
# ==============================================================================
#                             start computation
# ==============================================================================
@time ccall((:onetimestepforward_, "/Users/wenlei/Documents/ElasticWave/test/test.so"),
      Void, (Ptr{Int64}, Ptr{Int64}, Ptr{Int64},
             Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
             Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}},
             Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
             Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}},
             Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
             Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}}, Ptr{Array{Float64,2}},
             Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
             Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}),
             &nz, &nx, &nstep,
             &deltaZ, &deltaX, &dt,
             lambda, mu, rho,
             &hLambdaZ, &hLambda2MuZ, &hMuZ, &hMuX, &hRhoZX,
             sigmazz, sigmaxx, sigmazx, vz, vx,
             &vdvzdz, &vdvxdx, &vdvxdz, &vdvzdx, &vdsigmazzdz, &vdsigmazxdx, &vdsigmaxxdx, &vdsigmazxdz,
              mdvzdz,  mdvxdx,  mdvxdz,  mdvzdx,  mdsigmazzdz,  mdsigmazxdx,  mdsigmaxxdx,  mdsigmazxdz,
             az, bz, kz, haz, hbz, hkz,
             ax, bx, kx, hax, hbx, hkx)

             m = 7; n=6;
             A = zeros(Int64, m+2, n+2);
             ccall((:assign_, "/Users/wenlei/Documents/ElasticWave/test/test.so"),
                   Void, (Ptr{Array{Int64,2}}, Ptr{Int64}, Ptr{Int64}),
                          A                  , &m        , &n        )

# crazy

max_stars_repo_name: WneleiGao/ElasticWave


text:
"""
    Base.length(N::T) where {T <: AbstractEcologicalNetwork}

The length of a network is the number of non-zero elements in it.
"""
function Base.length(N::T) where {T <: AbstractEcologicalNetwork}
    return count(!iszero, N.edges)
end

"""
    similar(N::T) where {T <: AbstractEcologicalNetwork}

Returns a network with the same species, but an empty interaction matrix. This
is useful if you want to generate a "blank slate" for some analyses.
"""
function Base.similar(N::T) where {T <: AbstractEcologicalNetwork}
   Y = sparse(zeros(_interaction_type(N), size(N)))
   return T(Y, _species_objects(N)...)
end

"""
    size(N::AbstractEcologicalNetwork)

Return the size of the adjacency matrix of an `AbstractEcologicalNetwork` object.
"""
function Base.size(N::AbstractEcologicalNetwork)
  size(N.edges)
end

"""
    size(N::AbstractEcologicalNetwork, i::Int64)

Return the size of the adjacency matrix of an `AbstractEcologicalNetwork` object.
"""
function Base.size(N::AbstractEcologicalNetwork, i::Int64)
  size(N.edges, i)
end

"""
    getindex(N::AbstractEcologicalNetwork, i::T, j::T)

Get the value of an interaction based on the position of the species.
"""
function Base.getindex(N::AbstractEcologicalNetwork, i::T, j::T) where {T <: Int}
  return N.edges[i,j]
end

"""
    getindex(N::AbstractEcologicalNetwork, i::T)

Get the value of an interaction based on the position of the species.
"""
function Base.getindex(N::AbstractEcologicalNetwork, i::T) where {T <: Int}
  return N.edges[i]
end

"""
    getindex(N::AbstractEcologicalNetwork, ::Colon, j::T)

Get the value of an interaction based on the position of the species.
"""
function Base.getindex(N::AbstractEcologicalNetwork, ::Colon, j::T) where {T <: Int}
  return N.edges[:,j]
end

"""
    getindex(N::AbstractEcologicalNetwork, i::T, ::Colon)

Get the value of an interaction based on the position of the species.
"""
function Base.getindex(N::AbstractEcologicalNetwork, i::T, ::Colon) where {T <: Int}
  return N.edges[i,:]
end

"""
    getindex{T}(N::AbstractEcologicalNetwork, s1::T, s2::T)

Get the value of an interaction based on the *name* of the species. This is the
recommended way to look for things in a network.
"""
function Base.getindex(N::AbstractEcologicalNetwork, s1::T, s2::T) where {T}
  @assert T == _species_type(N)
  s1_pos = findfirst(isequal(s1), species(N; dims=1))
  s2_pos = findfirst(isequal(s2), species(N; dims=2))
  return N[s1_pos, s2_pos]
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, ::Colon, sp::T)

Gets the predecessors (*i.e.* species that interacts with / consume) of a focal
species. This returns the list of species as a `Set` object, in which ordering
is unimportant.
"""
function Base.getindex(N::AbstractEcologicalNetwork, ::Colon, sp::T) where {T}
  @assert T == _species_type(N)
  j = findfirst(isequal(sp), species(N; dims=2))
  i = findall(!iszero, N[:,j])
  return Set(species(N; dims=1)[i])
end

"""
    getindex{T}(N::AbstractEcologicalNetwork, sp::T, ::Colon)

Gets the successors (*i.e.* species that are interacted with / consumed) of a
focal species. This returns the list of species as a `Set` object, in which
ordering is unimportant.
"""
function Base.getindex(N::AbstractEcologicalNetwork, sp::T, ::Colon) where {T}
  @assert T == _species_type(N)
  i = findfirst(isequal(sp), species(N; dims=1))
  j = findall(!iszero, N[i,:])
  return Set(species(N; dims=2)[j])
end

"""
    getindex{T}(N::AbstractUnipartiteNetwork, sp::Array{T})

Induce a unipartite network based on a list of species, all of which must be in
the original network. This function takes a single argument (as opposed to two
arrays, or an array and a colon) to ensure that the returned network is
unipartite.

The network which is returned by this function may not have the species in the
order specified by the user for performance reasons.
"""
function Base.getindex(N::AbstractUnipartiteNetwork, sp::Vector{T}) where {T}
  @assert T == _species_type(N)
  sp_pos = indexin(sp, species(N))
  n_sp = species(N)[sp_pos]
  n_int = N.edges[sp_pos, sp_pos]
  return typeof(N)(n_int, n_sp)
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, ::Colon, sp::Array{T})

TODO
"""
function Base.getindex(N::AbstractBipartiteNetwork, ::Colon, sp::Vector{T}) where {T}
  @assert T == _species_type(N)
  sp_pos = Int.(something(indexin(sp, species(N; dims=2))))
  n_t = N.T
  n_b = N.B[sp_pos]
  n_int = N.edges[:, sp_pos]
  return typeof(N)(n_int, n_t, n_b)
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, sp::Array{T}, ::Colon)

TODO
"""
function Base.getindex(N::AbstractBipartiteNetwork, sp::Vector{T}, ::Colon) where {T}
  @assert T == _species_type(N)
  sp_pos = Int.(something(indexin(sp, species(N; dims=1))))
  n_t = N.T[sp_pos]
  n_b = N.B
  n_int = N.edges[sp_pos,:]
  return typeof(N)(n_int, n_t, n_b)
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, sp1::Array{T}, sp2::Array{T})

TODO
"""
function Base.getindex(N::AbstractBipartiteNetwork, sp1::Vector{T}, sp2::Vector{T}) where {T <: Integer}
   @assert all(sp1 .<= richness(N; dims=1))
   @assert all(sp2 .<= richness(N; dims=2))
   n_t = N.T[sp1]
   n_b = N.B[sp2]
   n_int = N.edges[sp1,sp2]
   return typeof(N)(n_int, n_t, n_b)
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, sp1::Array{T}, sp2::Array{T})

TODO
"""
function Base.getindex(N::AbstractBipartiteNetwork, sp1::Vector{T}, sp2::Vector{T}) where {T}
  @assert T == _species_type(N)
  sp1_pos = Int.(something(indexin(sp1, species(N; dims=1))))
  sp2_pos = Int.(something(indexin(sp2, species(N; dims=2))))
  return N[sp1_pos, sp2_pos]
end

"""
    setindex!(N::T, A::Any, i::E, j::E) where {T <: AbstractEcologicalNetwork, E}

Changes the value of the interaction at the specificied position, where `i` and
`j` are species *names*. Note that this operation **changes the network**.
"""
function Base.setindex!(N::T, A::Any, i::E, j::E) where {T <: AbstractEcologicalNetwork, E}
  @assert E == _species_type(N)
  @assert i ∈ species(N; dims=1)
  @assert j ∈ species(N; dims=2)
  i_pos = something(findfirst(isequal(i), species(N; dims=1)),0)
  j_pos = something(findfirst(isequal(j), species(N; dims=2)),0)
  N.edges[i_pos, j_pos] = A
end

"""
    setindex!(N::T, A::K, i::E, j::E) where {T <: AbstractEcologicalNetwork, K <: _interaction_type(N), E <: Int}

Changes the value of the interaction at the specificied position, where `i` and
`j` are species *positions*. Note that this operation **changes the network**.
"""
function Base.setindex!(N::T, A::Any, i::E, j::E) where {T <: AbstractEcologicalNetwork, E <: Int}
  @assert i ≤ richness(N; dims=1)
  @assert j ≤ richness(N; dims=2)
  N.edges[i, j] = A
end


text:
const m00 = -0.00006103098165
const n00 = 1.0
const m01 = 5.33967643357688
const n01 = 22.96302109010794
const m10 = -0.40661990365427
const n10 = -0.48466536361620
const m02 = 3.25023425332360
const n02 = -0.77268824532468
const m11 = -36.19405221599028
const n11 = -1.34102279982050
const m20 = 0.08975394404851
const n20 = 0.43027619553168
const m03 = 83.84593224417796
const n03 = -5.70531500645109
const m12 = 41.21772632732834
const n12 = 2.45782574294244
const m21 = 3.83815885394565
const n21 = -0.04763802358853
const m30 = -0.21619763215668
const n30 = -0.03326944290044

function impliedVolatilitySqrtTimeRationalGuess(x::Float64, c::Float64)::Float64
    num = m00 + x * (m10 + x * (m20 + x * m30)) + c * (m01 + m11 * x + m21 * x^2) + c^2 * (m02 + m12 * x) + c^3 * m03
    den = n00 + x * (n10 + x * (n20 + x * n30)) + c * (n01 + n11 * x + n21 * x^2) + c^2 * (n02 + n12 * x) + c^3 * n03
    return num / den
end

function impliedVolatilityLiRationalGuess(
    isCall::Bool,
    price::Float64,
    forward::Float64,
    strike::Float64,
    timeToExpiry::Float64,
    df::Float64,
)::Float64
    c, ex = normalizePrice(isCall,price,forward,strike,df)
    if c >= 1 / ex || c >= 1
        throw(DomainError(c, string("price higher than intrinsic value ", 1 / ex)))
    elseif c <= 0
        throw(DomainError(c, "price is negative"))
    end
    x = log(ex)
    return impliedVolatilitySqrtTimeRationalGuess(x, c) / sqrt(timeToExpiry), nil
end

abstract type SORSolver end
struct SORTS <: SORSolver end
struct SORDR <: SORSolver end

function impliedVolatilityLiSOR(
    isCall::Bool,
    price::T,
    forward::T,
    strike::T,
    timeToExpiry::T,
    df::T,
    impliedVolGuess::T,
    tolerance::T,
    maxIterations::Int,
    solver::SORSolver)::T where {T}
    c, ex = normalizePrice(isCall,price,forward,strike,df)
    if c > 1 / ex
        return 0
    end
    x = log(ex)
    v0 = T(0)
    if impliedVolGuess == 0
        # if abs(x) < 3 && c > 0.0005 && c < 0.9995 #domain of validity of the rational approximation
        #     v0 = impliedVolatilitySqrtTimeRationalGuess(x, c)
        #     #        println("guess R ",v0," ",v0/sqrtte)
        # else
            v0 = T(impliedVolatilitySRGuessUndiscountedCall(Float64(c), Float64(ex), Float64(x)))
            #    println("guess SR ",v0," ",v0/sqrtte,x,c)
        # end
    else
        v0 = impliedVolGuess * sqrt(timeToExpiry)
    end
    return computeLiSOR(c, x, ex, v0, sqrt(timeToExpiry), tolerance, maxIterations, solver)
end

function computeLiSOR(
    c::T,
    x::T,
    ex::T,
    v0::T,
    sqrttte::T,
    tolerance::T,
    maxIterations::Int,
    solver::SORSolver,
)::T where {T}
    v = v0
    if maxIterations <= 0
        return v
    end
    exinvsqrt = 1 / sqrt(ex)
    h = x / v
    t = v / 2
    sqrt2 = sqrt(T(2))
    Np = erfcx(-(h + t)/sqrt2)
    Nm = erfcx(-(h - t)/sqrt2)
    norm = exinvsqrt /2 * exp(-(h^2 + t^2)/2)
    cEstimate = norm * (Np - Nm)
    if cEstimate < 0
        cEstimate = T(0)
    end

    phi = (v^2 + 2 * x) / (v^2 - 2 * x)
    omega = omegaSOR(solver, phi)
    iterations = 0
    vOld = T(0)
    vTolerance = 64 * eps(T) * sqrttte #accuracy of inversion is limited, do not loop forever
    while abs(v - vOld) > vTolerance && abs(c - cEstimate) > tolerance && iterations < maxIterations
        vOld = v
        phi = (v^2 + 2 * x) / (v^2 - 2 * x)
        alpha = (1 + omega) / (1 + phi)
        F = c + norm * (Nm + omega * Np)
        Fom = F / (1 + omega)
        if Fom >= 1
            Fom = F / 2
        end
        Nm1 = T(0)
        #try
            Nm1 = norminv(Fom)
        #catch y
        #    #println(iterations, " ", v, " ", c, " ", F, " ", Fom, " ", omega)
        #    return v0 / sqrttte
        #end
        G = Nm1 + sqrt(Nm1^2 - 2 * x)
        v = iterateSOR(solver, v, G, alpha)
        omega = omegaSOR(solver, phi)
        h = x / v
        t = v / 2

        Np = erfcx(-(h + t)/sqrt2)
        Nm = erfcx(-(h - t)/sqrt2)
        norm = exinvsqrt * exp(-(h^2 + t^2)/2)/2
        # slower alternative:
        # Np = normcdf(h+t)
        # Nm = normcdf(h-t)/ex
        # norm = 1.0
        cEstimate = norm * (Np - Nm)
        if cEstimate < 0
            cEstimate = 0
        end
        iterations += 1
        #    println("SORTS ",iterations, " ",abs(v-vOld)," ",abs(c-cEstimate))
    end
    sigma = v / sqrttte
    return sigma
end

function iterateSOR(solver::SORTS, v::T, G::T, alpha::T)::T where {T}
    return alpha * G + (1 - alpha) * v
end
function iterateSOR(solver::SORDR, v::T, G::T, alpha::T)::T where {T}
    return G
end

function omegaSOR(solver::SORTS, phi::T)::T where {T}
    return 1
end

function omegaSOR(solver::SORDR, phi::T)::T where {T}
    return phi
end


text:
using XLStats

data_dir="./data"
xml_data="$data_dir/hitsDetail.xls"
#xml_data="$data_dir/salbiii_hitsDetail.dat"
links = read_all(xml_file=xml_data,
                 topolink_log="$data_dir/salbiii_topolink.log",
                 topolink_input="$data_dir/topolink.inp",
                 xic_file_name="$data_dir/salbiii_xic.dat",
                 domain=2:134)


text:
using TupleTools, Base.Cartesian
export loop_einsum, loop_einsum!, allow_loops

"""
    loop_einsum(::EinCode, xs, size_dict)

evaluates the eincode specified by `EinCode` and the tensors `xs` by looping
over all possible indices and calculating the contributions ot the result.
Scales exponentially in the number of distinct index-labels.
"""
function loop_einsum(code::EinCode{ixs, iy}, xs::NTuple{N, AbstractArray{<:Any,M} where M},
                size_dict) where {N, ixs, iy}
    size = getindex.(Ref(size_dict), iy)
    loop_einsum!(code, xs, get_output_array(xs, size; has_repeated_indices=!allunique(iy)), size_dict)
end

"""
    loop_einsum!(::EinCode, xs, y, size_dict)

inplace-version of `loop_einsum`, saving the result in a preallocated tensor
of correct size `y`.
"""
function loop_einsum!(code::EinCode{ixs, iy},
                xs::NTuple{N, AbstractArray{<:Any,M} where M},
                y::AbstractArray{T,L}, size_dict) where {N,L,T,IT <: Union{AbstractChar,Integer}, ixs, iy}
    ALLOW_LOOPS[] || @error "using `loop_einsum` to evaluate" code size.(xs) size(y)
    A = einarray(code, xs, size_dict)
    reduce_einarray!(A, y)
end

function reduce_einarray!(A::EinArray{T}, y) where T
    @inbounds for ind_y in A.OCIS
        iy = subindex(A.y_indexer,ind_y)
        yi = zero(T)
        for ind_x in A.ICIS
            ind = TupleTools.vcat(ind_x.I,ind_y.I)
            yi += map_prod(A.xs, ind, A.x_indexers)
        end
        y[iy] = yi
    end
    y
end

@inline function get_output_array(xs::NTuple{N, AbstractArray{<:Any,M} where M}, size; has_repeated_indices=true) where N
    if has_repeated_indices
        zeros(promote_type(map(eltype,xs)...), size...)
    else
        Array{promote_type(map(eltype,xs)...)}(undef, size...)
    end
end

const ALLOW_LOOPS = Ref(true)

"""
    allow_loops(flag::Bool)

Setting this to `false` will cause OMEinsum to log an error if it falls back to
`loop_einsum` evaluation, instead of calling specialised kernels. The default is `true`.
"""
function allow_loops(flag::Bool)
    ALLOW_LOOPS[] = flag
end


text:
"""
*Find the optimal parameters for quantum Bell violations.*

The `Nonlocality.jl` module provides tools to optimize quantum non-locality in
Bell scenarios.

### Module Exports:
* [`optimize_measurement`](@ref): Finds the quantum measurement which violates
    a specified Bell inequality maximally.
"""
module Nonlocality

using ..BellScenario

using QBase, LinearAlgebra

using Convex, SCS

include("./optimize_measurements.jl")

end


text:
function lotomania()
  i = 0; jogo = 0; sorteio = 1;

  while !(sorteio ⊆ jogo)
    jogo = sortear(100,50)
    sorteio = sortear(100,20)
    i += 1;
  end

  println("Sorteio: ", sorteio)
  println("Jogo: ", jogo)
  println("Numero de Jogos: ", i)
  return
end

function lotomania(t)
  i = 0; jogo = 0; sorteio = 1;

  if t == 0
    while 0 < length(∩(sorteio,jogo)) < 20
      jogo = sortear(100,50);
      sorteio = sortear(100,20);
      i += 1;
    end
  elseif t == 20
    while length(∩(sorteio,jogo)) < 20
      jogo = sortear(100,50);
      sorteio = sortear(100,20);
      i += 1;
    end
  else
    while 0 < length(∩(sorteio,jogo)) < t
      jogo = sortear(100,50);
      sorteio = sortear(100,20);
      i += 1;
    end
  end

  println("Sorteio: ", sorteio)
  println("Jogo de acerto: ", jogo, " (", length(∩(jogo,sorteio)), " acertos!)");
  println("Numero de Jogos: ", i)
  return
end

max_stars_repo_name: fabiom/loteria

max_issues_repo_name: fabiom/loteria

max_forks_repo_name: fabiom/loteria


text:
using Plots
using StaticArrays

# include("gp.jl")
# include("errors.jl")

function train_validate_test(𝒟_train, 𝒟_validate, 𝒟_test, problem; log_γs=-1.0:0.1:1.0, distances=[euclidean_distance, derivative_distance, antiderivative_distance],descriptor="")
    # Train GP on the filenames in train;
    # Optimize hyperparameter values by testing on filenames in validate;
    # Compute error on the filenames in test.

    nd = length(distances)
    nk = 4
    min_logγs       = zeros(nd,nk)
    validate_errors = zeros(nd,nk)
    test_errors     = zeros(nd,nk)

    for k in 1:nk, (i, d) in enumerate(distances)
        min_logγ, min_error_validate, test_error = get_min_gamma(k, d, 𝒟_train, 𝒟_validate, 𝒟_test; log_γs=log_γs)
        min_logγs[i,k]       = min_logγ
        validate_errors[i,k] = min_error_validate
        test_errors[i,k]     = test_error
    end

    # # for rational quadratic kernel, have 2 hyperparameters to optimize
    # k=5
    # for (i, d) in enumerate(distances)
    #     min_logγ, min_error_validate, test_error = get_min_gamma_alpha(k, d, 𝒟_train, 𝒟_validate, 𝒟_test; log_γs=log_γs)
    #     min_logγs[i,5]       = min_logγ
    #     validate_errors[i,5] = min_error_validate
    #     test_errors[i,5]     = test_error
    # end

    println("MIN LOG γs")
    println(min_logγs)
    println("VALIDATE Mean Error")
    println(validate_errors)
    println("TEST Mean Error")
    println(test_errors)

    a = argmin(test_errors)
    d = a[1]
    k = a[2]
    logγ    = min_logγs[d,k]
    kernel  = get_kernel(k, logγ, 0.0, distances[d])
    ℳ      = model(𝒟_train; kernel=kernel)
    anim    = animate_profile_and_model_output(ℳ, 𝒟_test)
    gif(anim, "$(descriptor)_$(typeof(problem))_$(problem.type)_kernel_$(k)_gamma_$(logγ).gif");

    println("===============")
    println("-- kernel ............. $(k)")
    println("-- norm ............... $(distances[d])")
    println("-- logγ ............... $(logγ)")
    println("-- validate error ..... $(validate_errors[d, k])")
    println("-- test error ......... $(test_errors[d, k])")

    return (min_logγs, validate_errors, test_errors)
end

function get_min_gamma(k, distance, 𝒟_train, 𝒟_validate, 𝒟_test; log_γs=-0.3:0.1:0.3)

    errors_validate = zeros(length(log_γs))

    for (i, logγ) in enumerate(log_γs)

        kernel = get_kernel(k, logγ, 0.0, distance)
        ℳ = model(𝒟_train; kernel=kernel)

        # -----compute mean error for true check----
        errors_validate[i] = get_me_true_check(ℳ, 𝒟_validate)
    end

    i                   = argmin(errors_validate)
    min_logγ            = log_γs[i]
    min_error_validate  = errors_validate[i]

    # using the log_γ value that minimizes the error on the validation set,
    # see how the model performs on the test set.
    kernel = get_kernel(k, min_logγ, 0.0, distance)
    ℳ = model(𝒟_train; kernel=kernel);
    error_test = get_me_true_check(ℳ, 𝒟_test)

    return (min_logγ, min_error_validate, error_test)
end

# function get_min_gamma(k::Int64, 𝒟::ProfileData, distance, log_γs)
#     # returns the gamma value that minimizes the mean error on the true check
#     # - only tests the gamma values listed in the log_γs parameter
#
#     mets  = zeros(length(log_γs)) # mean error for each gamma (true check)
#     for (i, logγ) in enumerate(log_γs)
#
#         kernel = get_kernel(k, logγ, 0.0, distance)
#         ℳ = model(𝒟; kernel=kernel);
#
#         # -----compute mean error for true check----
#         mets[i] = get_me_true_check(ℳ, 𝒟)
#     end
#
#     i = argmin(mets)
#     min_logγ = log_γs[i]
#     min_error = mets[i]
#
#     return (min_logγ, min_error)
# end

function get_min_gamma_alpha(k, distance, 𝒟_train, 𝒟_validate, 𝒟_test; log_γs=-0.3:0.1:0.3, log_αs=-0.3:0.1:0.3)
    # returns the gamma value that minimizes the mean error on the true check
    # only tests the gamma values listed in log_γs parameter

    errors_validate = @MArray zeros(length(log_γs), length(log_αs))

    for i in eachindex(log_γs), j in eachindex(log_αs)

        kernel = get_kernel(k, log_γs[i], 0.0, distance; logα=log_αs[j])
        ℳ = model(𝒟_train; kernel=kernel);

        # -----compute mean error for true check----
        errors_validate[i,j] = get_me_true_check(ℳ, 𝒟_validate)
    end

    m = argmin(errors_validate)
    min_logγ = log_γs[m[1]]
    min_logα = log_αs[m[2]]
    min_error_validate = errors_validate[m]

    # using the log_γ value that minimizes the error on the validation set,
    # see how the model performs on the test set.
    kernel = get_kernel(k, min_logγ, 0.0, distance; logα=min_logα)
    ℳ = model(𝒟_train; kernel=kernel);
    error_test = get_me_true_check(ℳ, 𝒟_test)

    return (min_logγ, min_error_validate, error_test)
end

function plot_landscapes_compare_error_metrics(k::Int64, 𝒟::ProfileData, distance, log_γs)
    # Compare mean log marginal likelihood with
    #    mean error on greedy check and
    #    mean error on true check

    mlls = zeros(length(log_γs)) # mean log marginal likelihood
    mes  = zeros(length(log_γs)) # mean error (greedy check)
    mets  = zeros(length(log_γs)) # mean error (true check)

    for i in 1:length(log_γs)

        kernel = get_kernel(k, log_γs[i], 0.0, distance)
        ℳ = model(𝒟; kernel=kernel)

        # -----compute mll loss----
        mlls[i] = -1*mean_log_marginal_loss(𝒟.y_train, ℳ, add_constant=false)

        # -----compute mean error for greedy check (same as in plot log error)----
        mes[i] = get_me_greedy_check(ℳ, 𝒟)

        # -----compute mean error for true check----
        mets[i] = get_me_true_check(ℳ, 𝒟)

    end

    ylims = ( minimum([minimum(mets), minimum(mes)]) , maximum([maximum(mets), maximum(mes)]) )

    mll_plot = plot(log_γs, mlls, xlabel="log(γ)", title="negative mean log marginal likelihood, P(y|X)", legend=false, yscale=:log10) # 1D plot: mean log marginal loss vs. γ
    vline!([log_γs[argmin(mlls)]])
    mes_plot  = plot(log_γs, mes,  xlabel="log(γ)", title="ME on greedy check, min = $(round(minimum(mes);digits=7))", legend=false, yscale=:log10, ylims=ylims)  # 1D plot: mean error vs. γ
    vline!([log_γs[argmin(mes)]])
    met_plot  = plot(log_γs, mets,  xlabel="log(γ)", title="ME on true check, min = $(round(minimum(mets);digits=7))", legend=false, yscale=:log10, ylims=ylims)  # 1D plot: mean error vs. γ
    vline!([log_γs[argmin(mets)]])

    layout = @layout [a; b; c]
    return plot(mll_plot, mes_plot, met_plot, layout = layout)
end

function plot_landscapes_compare_files_me(filenames, k::Int64, distance, log_γs, problem; D=16, N=4)
    # visual comparison of the mean error on true check for every file in filenames

    function get_me(filename)
        𝒟 = data(file, problem; D=D, N=N)

        mes  = zeros(length(log_γs))
        for i in 1:length(log_γs)
            kernel = get_kernel(k, log_γs[i], 0.0, distance)
            ℳ = model(𝒟; kernel=kernel)
            mes[i] = get_me_true_check(ℳ, 𝒟)
        end

        return mes
    end

    results = Dict(file => get_me(file) for file in filenames)

    # put all the data into one array for plotting
    for r in results
        all = hcat(all, r[file])
    end

    layout = (length(filenames), 1)
    ylims = (minimum(all),maximum(all))

    # minimizing γ values
    argmin_logγ = vcat([log_γs[argmin(results[file])]
                for file in filenames])

    titles = ["$(file), log(γ)=$(argmin_logγ[i]), min = $(round(minimum(results[filenames[i]]);digits=5))"
             for i in eachindex(filenames)]

    p = plot(log_γs, xlabel="log(γ)", ylabel="ME, true check", title=titles, legend=false, yscale=:log10, ylims=ylims, layout=layout)  # 1D plot: mean error vs. γ

    vline!(argmin_γ')

    return p
end

function plot_error_histogram(ℳ, 𝒟::ProfileData, time_index)
    # mean error for true check
    gpr_prediction = predict(ℳ, 𝒟; postprocessed=true)
    n = 𝒟.Nt-1

    gpr_error = zeros(n-1)
    for i in 1:n-1
        exact    = 𝒟.vavg[i+1]
        predi    = gpr_prediction[i+1]
        gpr_error[i] = euclidean_distance(exact, predi) # euclidean distance
    end
    mean_error = sum(gpr_error) / n

    error_plot_log = histogram(log.(gpr_error), title = "log(error) at each timestep of the full evolution", xlabel="log(Error)", ylabel="Frequency",ylims=(0,250), label="frequency")
    vline!([log(mean_error)], line = (4, :dash, 0.8), label="mean error")
    vline!([log(gpr_error[time_index])], line = (1, :solid, 0.6), label="error at t=$(time_index)")
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
"""
    DunklXuDisk(β)

are the orthogonal polynomials on the unit disk with respect to (1-x^2-y^2)^β, defined as `P_{n-k}^{k+β+1/2,k+β+1/2}(x)*(1-x^2)^{k/2}*P_k^{β,β}(y/sqrt{1-x^2})`.
"""
struct DunklXuDisk{T, V} <: BivariateOrthogonalPolynomial{T}
    β::V
end

DunklXuDisk(β::T) where T = DunklXuDisk{float(T), T}(β)
DunklXuDisk() = DunklXuDisk(0)

==(D1::DunklXuDisk, D2::DunklXuDisk) = D1.β == D2.β

axes(P::DunklXuDisk{T}) where T = (Inclusion(UnitDisk{T}()),blockedrange(oneto(∞)))

copy(A::DunklXuDisk) = A

Base.summary(io::IO, P::DunklXuDisk) = print(io, "DunklXuDisk($(P.β))")

"""
    DunklXuDiskWeight(β)

is a quasi-vector representing `(1-x^2-y^2)^β`.
"""
struct DunklXuDiskWeight{T, V} <: Weight{T}
    β::V
end

DunklXuDiskWeight(β::T) where T = DunklXuDiskWeight{float(T),T}(β)

axes(P::DunklXuDiskWeight{T}) where T = (Inclusion(UnitDisk{T}()),)

Base.summary(io::IO, P::DunklXuDiskWeight) = print(io, "(1-x^2-y^2)^$(P.β) on the unit disk")

const WeightedDunklXuDisk{T} = WeightedBasis{T,<:DunklXuDiskWeight,<:DunklXuDisk}

WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β)

@simplify function *(Dx::PartialDerivative{1}, P::DunklXuDisk)
    β = P.β
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockBroadcastArray(hcat,
        ((k .+ (β - 1)) ./ (2k .+ (2β - 1)) .* (n .- k .+ 1)),
        Zeros((axes(n,1),)),
        (((k .+ 2β) .* (k .+ (2β + 1))) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)) .* (n .+ k .+ 2β) ./ 2)
        )
    DunklXuDisk(β+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,2))
end

@simplify function *(Dy::PartialDerivative{2}, P::DunklXuDisk)
    β = P.β
    k = mortar(Base.OneTo.(oneto(∞)))
    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
    DunklXuDisk(β+1) * _BandedBlockBandedMatrix(((k .+ T(2β)) ./ 2)', axes(k,1), (-1,1), (-1,1))
end

@simplify function *(Dx::PartialDerivative{1}, w_P::WeightedDunklXuDisk)
    wP, P = w_P.args
    @assert P.β == wP.β
    β = P.β
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockBroadcastArray(hcat,
        (-4 .* (k .+ (β - 1)).*(n .- k .+ 1)./(2k .+ (2β - 1))),
        Zeros((axes(n,1),)),
        (-k .* (k .+ 1) ./ ((2k .+ (2β - 1)) .* (k .+ β)) .* (n .+ k .+ 2β))
        )
    WeightedDunklXuDisk(β-1) * _BandedBlockBandedMatrix(dat', axes(k,1), (1,-1), (2,0))
end

@simplify function *(Dy::PartialDerivative{2}, w_P::WeightedDunklXuDisk)
    wP, P = w_P.args
    @assert P.β == wP.β
    k = mortar(Base.OneTo.(oneto(∞)))
    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
    WeightedDunklXuDisk(P.β-1) * _BandedBlockBandedMatrix((T(-2).*k)', axes(k,1), (1,-1), (1,-1))
end

# P^{(β)} ↗ P^{(β+1)}
function dunklxu_raising(β)
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockHcat(
        BlockBroadcastArray(hcat,
            (-(k .+ (β - 1)) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (4n .+ 4β))), # n-2, k-2
            Zeros((axes(n,1),)), # n-2, k-1
            (-(k .+ 2β) .* (k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)) .* (2n .+ (2β-1)) ./ (8n .+ 8β))), # n-2, k
        Zeros((axes(n,1), Base.OneTo(3))),
        BlockBroadcastArray(hcat,
            (2 .* (k .+ (β - 1)) .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2k .+ (2β-1)) .* (2n .+ 2β) .* (2n .+ (2β + 1)))), # n, k-2
            Zeros((axes(n,1),)), # n, k-1
            ((k .+ 2β) .* (k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)) .* (n .+ k .+ 2β) .* (n .+ k .+ (2β + 1)) ./ ((2n .+ 2β) .* (2n .+ (2β + 1)))) # n, k
        ))
    _BandedBlockBandedMatrix(dat', axes(k,1), (0,2), (0,2))
end

function \(A::DunklXuDisk, B::DunklXuDisk)
    if A.β == B.β
        Eye((axes(B,2),))
    elseif A.β == B.β + 1
        dunklxu_raising(B.β)
    else
        error("not implemented for $A and $B")
    end
end

# (1-x^2-y^2) P^{(β)} ↘ P^{(β-1)}
function dunklxu_lowering(β)
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockHcat(
        BlockBroadcastArray(hcat,
            ((2n .+ (2β - 1)) ./ (n .+ β) .* (k .+ (β - 1)) ./ (2k .+ (2β - 1))), # n, k
            Zeros((axes(n,1),)), # n, k+1
            ((n .+ (β - 1/2)) ./ (n .+ β) .* k .* (k .+ 1) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)))), # n, k+2
        Zeros((axes(n,1), Base.OneTo(3))),
        BlockBroadcastArray(hcat,
            (-8 .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2n .+ 2β) .* (2n .+ (2β + 1))) .* (k .+ (β - 1)) ./(2k .+ (2β - 1))), # n+2, k
            Zeros((axes(n,1),)), # n+2, k+1
            (-4 .* (n .+ k .+ 2β) .* (n .+ k .+ (2β + 1)) ./ ((2n .+ 2β) .* (2n .+ (2β + 1))) .* k .* (k .+ 1) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β))) # n+2, k+2
        ))
    _BandedBlockBandedMatrix(dat', axes(k,1), (2,0), (2,0))
end

function \(w_A::WeightedDunklXuDisk, w_B::WeightedDunklXuDisk)
    wA,A = w_A.args
    wB,B = w_B.args

    @assert wA.β == A.β
    @assert wB.β == B.β

    if A.β == B.β
        Eye((axes(B,2),))
    elseif A.β + 1 == B.β
        dunklxu_lowering(B.β)
    else
        error("not implemented for $A and $B")
    end
end

\(w_A::DunklXuDisk, w_B::WeightedDunklXuDisk) =
    (DunklXuDiskWeight(0) .* w_A) \ w_B

# Actually Jxᵀ
function jacobimatrix(::Val{1}, P::DunklXuDisk)
    β = P.β
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockHcat(
        ((2n .+ (2β - 1)) ./ (4n .+ 4β)), # n-1, k
        Zeros((axes(n,1),)), # n, k
        ((n .- k .+ 1) .* (n .+ k .+ 2β) ./ ((n .+ β) .* (2n .+ (2β + 1)))), # n+1, k
        )
    _BandedBlockBandedMatrix(dat', axes(k,1), (1,1), (0,0))
end

# Actually Jyᵀ
function jacobimatrix(::Val{2}, P::DunklXuDisk)
    β = P.β
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockHcat(
        BlockBroadcastArray(hcat,
            ((k .+ (β - 1)) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (2n .+ 2β))), # n-1, k-1
            Zeros((axes(n,1),)), # n-1, k
            (-k .* (k .+ 2β) .* (2n .+ (2β - 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β) .* (4n .+ 4β)))), # n-1, k+1
        Zeros((axes(n,1),Base.OneTo(3))),
        BlockBroadcastArray(hcat,
            (-(2k .+ (2β - 2)) .* (n .- k .+ 1) .* (n .- k .+ 2) ./ ((2k .+ (2β - 1)) .* (n .+ β) .* (2n .+ (2β + 1)))), # n+1, k-1
            Zeros((axes(n,1),)), # n+1, k
            (k .* (k .+ 2β) .* (n .+ k .+ 2β) .* (n .+ k .+ (2β + 1)) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β) .* (n .+ β) .* (2n .+ (2β + 1)))))) # n+1, k+1
    _BandedBlockBandedMatrix(dat', axes(k,1), (1,1), (1,1))
end

@simplify function *(A::AngularMomentum, P::DunklXuDisk)
    β = P.β
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockBroadcastArray(hcat,
        (2 .* (k .+ (β - 1)) .* (n .- k .+ 1) ./ (2k .+ (2β - 1))), # n, k-1
        Zeros((axes(n,1),)), # n, k
        (-k .* (k .+ 2β) .* (n .+ k .+ 2β) ./ ((2k .+ (2β - 1)) .* (2k .+ 2β)))) # n, k+1

    DunklXuDisk(β) * _BandedBlockBandedMatrix(dat', axes(k,1), (0,0), (1,1))
end


text:
] add SumOfSquares
] add DynamicPolynomials
] add Mosek

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
using Decimals
using Base.Test

# Normalization
@test Decimal(1, 151100, -4) == Decimal(1, 1511, -2)
@test Decimal(0, 100100, -5) == Decimal(0, 1001, -3)
@test norm(Decimal(1, 151100, -4)) == Decimal(1, 1511, -2)
@test norm(Decimal(0, 100100, -5)) == Decimal(0, 1001, -3)
@test decimal("3.0") == Decimal(0, 3, 0)
@test decimal("3.0") == Decimal(0, 30, -1)
@test decimal("3.1400") == Decimal(0, 314, -2)
@test decimal("1234") == Decimal(0, 1234, 0)


text:
if Pkg.installed("Gurobi") != nothing
    import Gurobi
end

immutable Razaviyayn2013_MinMaxWMMSEState
    U::Array{Matrix{Complex128},1}
    W::Array{Hermitian{Complex128},1}
    V::Array{Matrix{Complex128},1}
end

function Razaviyayn2013_MinMaxWMMSE(channel, network)
    assignment = get_assignment(network)

    # The implementation is currently limited in the respects below. This is in
    # order to simplify the Gurobi optimization variable indexing. With equal
    # number of antennas and streams, it is very easy to calculate the variable
    # offsets required. (See functions v_ind and mu_ind.)
    require_equal_num_BS_antennas(network)
    require_equal_num_streams(network)
    require_equal_num_MSs_per_cell(assignment)

    K = get_num_MSs(network)
    Ps = get_transmit_powers(network)
    sigma2s = get_receiver_noise_powers(network)
    ds = get_num_streams(network); max_d = maximum(ds)
    alphas = get_user_priorities(network)

    aux_params = get_aux_precoding_params(network)
    @defaultize_param! aux_params "Razaviyayn2013_MinMaxWMMSE:Gurobi_verbose" false
    @defaultize_param! aux_params "Razaviyayn2013_MinMaxWMMSE:Gurobi_PSDTol" 1e-5

    state = Razaviyayn2013_MinMaxWMMSEState(
        Array(Matrix{Complex128}, K),
        Array(Hermitian{Complex128}, K),
        initial_precoders(channel, Ps, sigma2s, ds, assignment, aux_params))
    objective = Float64[]
    logdet_rates = zeros(Float64, K, max_d, aux_params["max_iters"])
    MMSE_rates = zeros(Float64, K, max_d, aux_params["max_iters"])
    weighted_logdet_rates = zeros(Float64, K, max_d, aux_params["max_iters"])
    weighted_MMSE_rates = zeros(Float64, K, max_d, aux_params["max_iters"])
    allocated_power = zeros(Float64, K, max_d, aux_params["max_iters"])

    if Pkg.installed("Gurobi") == nothing
        iters = 1; objective = 0.
        Lumberjack.warn("Razaviyayn2013_MinMaxWMMSE exiting because Gurobi.jl is not installed.")
    else
        iters = 0; conv_crit = Inf
        while iters < aux_params["max_iters"]
            update_MSs!(state, channel, sigma2s, assignment)
            iters += 1

            # Results after this iteration
            logdet_rates[:,:,iters] = calculate_logdet_rates(state)
            push!(objective, minimum(sum(logdet_rates[:,:,iters], 2)))
            MMSE_rates[:,:,iters] = calculate_MMSE_rates(state)
            weighted_logdet_rates[:,:,iters] = calculate_weighted_logdet_rates(state, alphas)
            weighted_MMSE_rates[:,:,iters] = calculate_weighted_MMSE_rates(state, alphas)
            allocated_power[:,:,iters] = calculate_allocated_power(state)

            # Check convergence
            if iters >= 2
                conv_crit = abs(objective[end] - objective[end-1])/abs(objective[end-1])
                if conv_crit < aux_params["stop_crit"]
                    Lumberjack.debug("Razaviyayn2013_MinMaxWMMSE converged.",
                        @compat Dict(
                            :num_iters => iters,
                            :final_objective => objective[end],
                            :conv_crit => conv_crit,
                            :stop_crit => aux_params["stop_crit"],
                            :max_iters => aux_params["max_iters"]))
                    break
                end
            end

            # Begin next iteration, unless the loop will end
            if iters < aux_params["max_iters"]
                update_BSs!(state, channel, Ps, sigma2s, assignment, aux_params)
            end
        end
        if iters == aux_params["max_iters"]
            Lumberjack.debug("Razaviyayn2013_MinMaxWMMSE did NOT converge.",
                @compat Dict(
                    :num_iters => iters,
                    :final_objective => objective[end],
                    :conv_crit => conv_crit,
                    :stop_crit => aux_params["stop_crit"],
                    :max_iters => aux_params["max_iters"]))
        end
    end

    results = PrecodingResults()
    if aux_params["output_protocol"] == :all_iterations
        results["objective"] = objective
        results["logdet_rates"] = logdet_rates
        results["MMSE_rates"] = MMSE_rates
        results["weighted_logdet_rates"] = weighted_logdet_rates
        results["weighted_MMSE_rates"] = weighted_MMSE_rates
        results["allocated_power"] = allocated_power
    elseif aux_params["output_protocol"] == :final_iteration
        results["objective"] = objective[iters]
        results["logdet_rates"] = logdet_rates[:,:,iters]
        results["MMSE_rates"] = MMSE_rates[:,:,iters]
        results["weighted_logdet_rates"] = weighted_logdet_rates[:,:,iters]
        results["weighted_MMSE_rates"] = weighted_MMSE_rates[:,:,iters]
        results["allocated_power"] = allocated_power[:,:,iters]
    end
    return results
end

function update_MSs!(state::Razaviyayn2013_MinMaxWMMSEState,
    channel::SinglecarrierChannel, sigma2s, assignment)

    ds = [ size(state.V[k], 2) for k = 1:channel.K ]

    for i = 1:channel.I; for k in served_MS_ids(i, assignment)
        Phi = complex(sigma2s[k]*eye(channel.Ns[k]))
        for j = 1:channel.I; for l in served_MS_ids(j, assignment)
            F = channel.H[k,j]*state.V[l]
            Phi += F*F'
        end; end

        # MMSE receiver and optimal MSE weight
        F = channel.H[k,i]*state.V[k]
        state.U[k] = Hermitian(Phi)\F
        Emmse = UniformScaling(1.) - state.U[k]'*F
        state.W[k] = inv(Hermitian((Emmse + Emmse')/2))
    end; end
end

# Some of the vectors describing the power constraints and rate constraints
# could be prebuilt and reused in the iterations. Probably no major performance
# boost would come from that though, based on experience with the
# JointPrecodingMCSSelection implementation.
function update_BSs!(state::Razaviyayn2013_MinMaxWMMSEState,
    channel::SinglecarrierChannel, Ps, sigma2s, assignment, aux_params)

    # The following is OK, since we have checked this at the top of
    # the outer function.
    K = channel.K; I = channel.I; Kc = @compat Int(channel.K/channel.I)
    N = channel.Ns[1]; M = channel.Ms[1]; d = size(state.W[1], 2)

    # Bookkeeping for the Gurobi model
    num_t_vars = 1
    num_v_vars = 2*K*M*d
    num_lin_vars_per_rate_constraint = 2*M*d + 1 # including t variable
    num_quad_vars_per_rate_constraint = K*d*(2*M)^2
    num_vars_per_power_constraint = 2*Kc*M*d

    # Index functions for optimization variables
    t_ind = 1
    v_ind(k, m, n, imag::Bool) = (num_t_vars + (k-1)*M*d*2 + (m-1)*d*2 + (n-1)*2 + 1 + @compat Int(imag))

    # Gurobi environment
    env = Gurobi.Env()
    if !aux_params["Razaviyayn2013_MinMaxWMMSE:Gurobi_verbose"]
        Gurobi.setparam!(env, "OutputFlag", 0)
    end
    Gurobi.setparam!(env, "PSDTol", aux_params["Razaviyayn2013_MinMaxWMMSE:Gurobi_PSDTol"])
    model = Gurobi.Model(env, "Razaviyayn2013_MinMaxWMMSE", finalize_env=true)
    Gurobi.set_sense!(model, :maximize)

    # Objective
    Gurobi.add_cvar!(model, 1., 0, Inf) # t

    # Other variables not in objective
    Gurobi.add_cvars!(model, zeros(Float64, num_v_vars))

    # Need to update model before adding constraints
    Gurobi.update_model!(model)

    # Rate constraints
    rate_constraint_lind = Array(Int, num_lin_vars_per_rate_constraint)
    rate_constraint_lval = Array(Float64, num_lin_vars_per_rate_constraint)
    rate_constraint_qr = Array(Int, num_quad_vars_per_rate_constraint)
    rate_constraint_qc = Array(Int, num_quad_vars_per_rate_constraint)
    rate_constraint_qv = Array(Float64, num_quad_vars_per_rate_constraint)
    for i = 1:I; for k = served_MS_ids(i, assignment)
        # Linear part
        g = (channel.H[k,i]'*state.U[k])*state.W[k]
        l_ind = 1
        for m = 1:M; for n = 1:d
            real_ind = v_ind(k, m, n, false)
            rate_constraint_lind[l_ind] = real_ind
            rate_constraint_lval[l_ind] = -2.0*real(g[m,n])
            l_ind += 1

            imag_ind = v_ind(k, m, n, true)
            rate_constraint_lind[l_ind] = imag_ind
            rate_constraint_lval[l_ind] = -2.0*imag(g[m,n])
            l_ind += 1
        end; end
        rate_constraint_lind[l_ind] = t_ind
        rate_constraint_lval[l_ind] = 1.

        # Quadratic part
        # Enforce PSDness explicitly with get_Gr and get_Gi.
        # Hopefully these are inlined.
        q_ind = 1
        for j = 1:I
            G = (channel.H[k,j]'*state.U[k])*state.W[k]*(state.U[k]'*channel.H[k,j])

            # symmetric
            Gr = real(G)
            get_Gr(mm1, mm2) = ( mm1 >= mm2 ? Gr[mm1,mm2] :  Gr[mm2,mm1] )

            # skew-symmetric
            Gi = imag(G)
            get_Gi(mm1, mm2) = ( mm1 >= mm2 ? Gi[mm1,mm2] : -Gi[mm2,mm1] )

            for l = served_MS_ids(j, assignment)
                for nn = 1:d
                    for mm1 = 1:M
                        real_ind1 = v_ind(l, mm1, nn, false)
                        imag_ind1 = v_ind(l, mm1, nn, true)
                        for mm2 = 1:M
                            real_ind2 = v_ind(l, mm2, nn, false)
                            imag_ind2 = v_ind(l, mm2, nn, true)

                            # real(v)'*real(Q)*real(v)
                            rate_constraint_qr[q_ind] = real_ind1
                            rate_constraint_qc[q_ind] = real_ind2
                            rate_constraint_qv[q_ind] = get_Gr(mm1, mm2)
                            q_ind += 1

                            # imag(v)'*real(Q)*imag(v)
                            rate_constraint_qr[q_ind] = imag_ind1
                            rate_constraint_qc[q_ind] = imag_ind2
                            rate_constraint_qv[q_ind] = get_Gr(mm1, mm2)
                            q_ind += 1

                            # -real(v)'*imag(Q)*imag(v)
                            rate_constraint_qr[q_ind] = real_ind1
                            rate_constraint_qc[q_ind] = imag_ind2
                            rate_constraint_qv[q_ind] = -get_Gi(mm1, mm2)
                            q_ind += 1

                            # imag(v)'*imag(Q)*real(v)
                            rate_constraint_qr[q_ind] = imag_ind1
                            rate_constraint_qc[q_ind] = real_ind2
                            rate_constraint_qv[q_ind] = get_Gi(mm1, mm2)
                            q_ind += 1
                        end
                    end
                end
            end
        end

        # Constant part
        rhs = logabsdet(state.W[k]) - abs(trace(state.W[k]*(UniformScaling(1.) + sigma2s[k]*state.U[k]'*state.U[k]))) + d

        Gurobi.add_qconstr!(model,
            rate_constraint_lind,
            rate_constraint_lval,
            rate_constraint_qr,
            rate_constraint_qc,
            rate_constraint_qv,
            '<',
            rhs)
    end; end

    # Power constraints
    power_constraint_qr = Array(Int, num_vars_per_power_constraint)
    power_constraint_qc = Array(Int, num_vars_per_power_constraint)
    power_constraint_qv = Array(Float64, num_vars_per_power_constraint)
    for i = 1:I
        q_ind = 1
        for k = served_MS_ids(i, assignment)
            for m = 1:M; for n = 1:d
                real_ind = v_ind(k, m, n, false)
                power_constraint_qr[q_ind] = real_ind
                power_constraint_qc[q_ind] = real_ind
                power_constraint_qv[q_ind] = 1.
                q_ind += 1

                imag_ind = v_ind(k, m, n, true)
                power_constraint_qr[q_ind] = imag_ind
                power_constraint_qc[q_ind] = imag_ind
                power_constraint_qv[q_ind] = 1.
                q_ind += 1
            end; end
        end

        Gurobi.add_qconstr!(model,
            [], # lind
            [], # lval
            power_constraint_qr,
            power_constraint_qc,
            power_constraint_qv,
            '<',
            Ps[i])
    end

    # Solve it!
    Gurobi.update_model!(model)
    Gurobi.optimize(model)

    if Gurobi.get_status(model) == :optimal || Gurobi.get_status(model) == :suboptimal
        if Gurobi.get_status(model) == :suboptimal
            Lumberjack.warn("Suboptimal precoder solution found in Razaviyayn2013_MinMaxWMMSE")
        end

        sol = Gurobi.get_solution(model)

        # Put into state variable
        for i = 1:channel.I
            for k in served_MS_ids(i, assignment)
                for m = 1:M
                    for n = 1:d
                        real_ind = v_ind(k, m, n, false)
                        imag_ind = v_ind(k, m, n, true)
                        state.V[k][m,n] = sol[real_ind] + im*sol[imag_ind]
                    end
                end
            end
        end
    else
        Lumberjack.warn("Numerical problems with Gurobi in Razaviyayn2013_MinMaxWMMSE", @compat Dict(:gurobi_output => string(Gurobi.get_status(model))))
    end
end


text:
# The set-up was designed for the
# 8th International Cloud Modelling Workshop
# (ICMW, Muhlbauer et al., 2013, case 1, doi:10.1175/BAMS-D-12-00188.1)
#
# See chapter 2 in Arabas et al 2015 for setup details:
#@Article{gmd-8-1677-2015,
#AUTHOR = {Arabas, S. and Jaruga, A. and Pawlowska, H. and Grabowski, W. W.},
#TITLE = {libcloudph++ 1.0: a single-moment bulk, double-moment bulk,
#         and particle-based warm-rain microphysics library in C++},
#JOURNAL = {Geoscientific Model Development},
#VOLUME = {8},
#YEAR = {2015},
#NUMBER = {6},
#PAGES = {1677--1707},
#URL = {https://www.geosci-model-dev.net/8/1677/2015/},
#DOI = {10.5194/gmd-8-1677-2015}
#}

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

using ClimateMachine
ClimateMachine.init()
using ClimateMachine.Atmos
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.Grids
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Thermodynamics:
    gas_constants,
    PhaseEquil,
    PhasePartition_equil,
    PhasePartition,
    internal_energy,
    q_vap_saturation,
    relative_humidity,
    TemperatureSHumEquil,
    TemperatureSHumNonEquil,
    air_temperature
using ClimateMachine.Microphysics
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.VariableTemplates
using ClimateMachine.VTK

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cp_d, cv_d, cv_v, T_0, e_int_v0, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import ClimateMachine.DGMethods:
    BalanceLaw,
    DGModel,
    LocalGeometry,
    vars_state_conservative,
    vars_state_auxiliary,
    vars_state_gradient,
    vars_state_gradient_flux,
    init_state_conservative!,
    init_state_auxiliary!,
    update_auxiliary_state!,
    nodal_update_auxiliary_state!,
    flux_first_order!,
    flux_second_order!,
    wavespeed,
    boundary_state!,
    source!

struct KinematicModelConfig{FT}
    xmax::FT
    ymax::FT
    zmax::FT
    wmax::FT
    θ_0::FT
    p_0::FT
    p_1000::FT
    qt_0::FT
    z_0::FT
end

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

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

    @assert param_set ≠ nothing
    @assert init_state_conservative ≠ nothing

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

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

vars_state_gradient(m::KinematicModel, FT) = @vars()

vars_state_gradient_flux(m::KinematicModel, FT) = @vars()

function init_state_auxiliary!(
    m::KinematicModel,
    aux::Vars,
    geom::LocalGeometry,
)

    FT = eltype(aux)
    x, y, z = geom.coord
    dc = m.data_config

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

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

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

function init_state_conservative!(
    m::KinematicModel,
    state::Vars,
    aux::Vars,
    coords,
    t,
    args...,
)
    m.init_state_conservative(m, state, aux, coords, t, args...)
end

function update_auxiliary_state!(
    dg::DGModel,
    m::KinematicModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    nodal_update_auxiliary_state!(
        kinematic_model_nodal_update_auxiliary_state!,
        dg,
        m,
        Q,
        t,
        elems,
    )
    return true
end

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

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

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

function config_kinematic_eddy(
    FT,
    N,
    resolution,
    xmax,
    ymax,
    zmax,
    wmax,
    θ_0,
    p_0,
    p_1000,
    qt_0,
    z_0,
)
    # Choose explicit solver
    ode_solver = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    kmc = KinematicModelConfig(
        FT(xmax),
        FT(ymax),
        FT(zmax),
        FT(wmax),
        FT(θ_0),
        FT(p_0),
        FT(p_1000),
        FT(qt_0),
        FT(z_0),
    )

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

    config = ClimateMachine.AtmosLESConfiguration(
        "KinematicModel",
        N,
        resolution,
        FT(xmax),
        FT(ymax),
        FT(zmax),
        param_set,
        init_kinematic_eddy!;
        solver_type = ode_solver,
        model = model,
    )

    return config
end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
using LinearAlgebra
import Adapt

@testset "constructors" begin
  xs = oneArray{Int}(undef, 2, 3)
  @test collect(oneArray([1 2; 3 4])) == [1 2; 3 4]
  @test testf(vec, rand(5,3))
  @test Base.elsize(xs) == sizeof(Int)
  @test oneArray{Int, 2}(xs) === xs

  @test_throws ArgumentError Base.unsafe_convert(Ptr{Int}, xs)
  @test_throws ArgumentError Base.unsafe_convert(Ptr{Float32}, xs)

  @test collect(oneAPI.zeros(2, 2)) == zeros(2, 2)
  @test collect(oneAPI.ones(2, 2)) == ones(2, 2)

  @test collect(oneAPI.fill(0, 2, 2)) == zeros(2, 2)
  @test collect(oneAPI.fill(1, 2, 2)) == ones(2, 2)
end

@testset "adapt" begin
  A = rand(Float32, 3, 3)
  dA = oneArray(A)
  @test Adapt.adapt(Array, dA) == A
  @test Adapt.adapt(oneArray, A) isa oneArray
  @test Array(Adapt.adapt(oneArray, A)) == A
end

@testset "reshape" begin
  A = [1 2 3 4
       5 6 7 8]
  gA = reshape(oneArray(A),1,8)
  _A = reshape(A,1,8)
  _gA = Array(gA)
  @test all(_A .== _gA)
  A = [1,2,3,4]
  gA = reshape(oneArray(A),4)
end

@testset "fill(::SubArray)" begin
  xs = oneAPI.zeros(3)
  fill!(view(xs, 2:2), 1)
  @test Array(xs) == [0,1,0]
end


text:
using ForwardDiff
using ForwardDiff: Dual

function test_derivative(f, x::T) where T
  buf = CuArray(zeros(T))

  function kernel(a, x)
    a[] = ForwardDiff.derivative(f, x)
    return
  end
  @cuda kernel(buf, x)
  return CUDA.@allowscalar buf[]
end

testf(cuf, f, x) = test_derivative(cuf, x) ≈ ForwardDiff.derivative(f, x)

@testset "UNARY" begin
  fs = filter(x->x[1] ==:CUDA && x[3] == 1, keys(ForwardDiff.DiffRules.DEFINED_DIFFRULES))

  nonneg = [:log, :log1p, :log10, :log2, :sqrt, :acosh]

  for (m, fn, _) ∈ fs
    cuf = @eval $m.$fn
    f = @eval $fn

    x32 = rand(Float32)
    x64 = rand(Float64)
    nx32 = -x32
    nx64 = -x64

    if fn == :acosh
      x32 += 1
      x64 += 1
    end

    @test testf(cuf, f, x32)
    @test testf(cuf, f, x64)

    if fn ∉ nonneg
      @test testf(cuf, f, nx32)
      @test testf(cuf, f, nx64)
    end
  end
end

@testset "POW" begin
  x32 = rand(Float32)
  x64 = rand(Float64)
  y32 = rand(Float32)
  y64 = rand(Float64)
  y = Int32(7)

  @test testf(x->CUDA.pow(x, Int32(7)), x->x^y, x32)
  @test testf(x->CUDA.pow(x, y), x->x^y, x64)
  @test testf(x->CUDA.pow(x, y32), x->x^y32, x32)
  @test testf(x->CUDA.pow(x, y64), x->x^y64, x64)

  @test testf(y->CUDA.pow(x32, y), y->x32^y, y32)
  @test testf(y->CUDA.pow(x64, y), y->x64^y, y64)

  @test testf(x->CUDA.pow(x, x), x->x^x, x32)
  @test testf(x->CUDA.pow(x, x), x->x^x, x64)
end

@testset "LITERAL_POW" begin
  for x in [rand(Float32, 10), rand(Float64, 10)],
      p in [1, 2, 3, 4, 5]
    @test ForwardDiff.gradient(_x -> sum(_x .^ p), x) ≈ p .* (x .^ (p - 1))
  end
end

@testset "Broadcast Fix" begin
  if CUDA.has_cudnn()
    using NNlib

    f(x) = logσ.(x)
    ds = Dual.(rand(5),1)
    @test f(ds) ≈ collect(f(CuArray(ds)))
  end
end


text:
include("variable.jl")
include("constraint.jl")

function MOIU.allocate(optimizer::SOItoMOIBridge, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    # To be sure that it is done before load(optimizer, ::ObjectiveFunction, ...), we do it in allocate
    optimizer.objsign = sense == MOI.MIN_SENSE ? -1 : 1
end
function MOIU.allocate(::SOItoMOIBridge, ::MOI.ObjectiveFunction, ::Union{MOI.SingleVariable, MOI.ScalarAffineFunction}) end

function MOIU.load(::SOItoMOIBridge, ::MOI.ObjectiveSense, ::MOI.OptimizationSense) end
# Loads objective coefficient α * vi
function load_objective_term!(optimizer::SOItoMOIBridge, α, vi::MOI.VariableIndex)
    for (blk, i, j, coef, shift) in varmap(optimizer, vi)
        if !iszero(blk)
            # in SDP format, it is max and in MPB Conic format it is min
            setobjectivecoefficient!(optimizer.sdoptimizer, optimizer.objsign * coef * α, blk, i, j)
        end
        optimizer.objshift += α * shift
    end
end
function MOIU.load(optimizer::SOItoMOIBridge, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction)
    obj = MOIU.canonical(f)
    optimizer.objconstant = f.constant
    for t in obj.terms
        if !iszero(t.coefficient)
            load_objective_term!(optimizer, t.coefficient, t.variable_index)
        end
    end
end
function MOIU.load(optimizer::SOItoMOIBridge{T}, ::MOI.ObjectiveFunction, f::MOI.SingleVariable) where T
    load_objective_term!(optimizer, one(T), f.variable)
end

function MOIU.allocate_variables(optimizer::SOItoMOIBridge{T}, nvars) where T
    optimizer.free = BitSet(1:nvars)
    optimizer.varmap = Vector{Vector{Tuple{Int, Int, Int, T, T}}}(undef, nvars)
    VI.(1:nvars)
end

function MOIU.load_variables(optimizer::SOItoMOIBridge, nvars)
    @assert nvars == length(optimizer.varmap)
    loadfreevariables!(optimizer)
    init!(optimizer.sdoptimizer, optimizer.blockdims, optimizer.nconstrs)
end


text:
mutable struct RandWalker
  i::Int
  x::Float64
  y::Float64
  ymin::Float64
  reporter
  function RandWalker(reporter = nothing)
    new(0, 0.0, 0.0, Inf, reporter)
  end
end

function step!(rw::RandWalker)
  rw.i += 1
  rw.x += randn()
  rw.y += (randn() - 0.05) # Slight downward trend to make more fitness improvement happen...
  if rw.y < rw.ymin
    println("New min found: ($(rw.i),$(rw.y))")
    rw.ymin = rw.y
    if rw.reporter != nothing
      rw.reporter(rw.i, rw.ymin)
    end
  end
end

function rand_walk!(rw::RandWalker)
  while true
    print("."); flush(stdout);
    step!(rw)
    sleep(0.4 * rand())
  end
end

using Gadfly

function fitness_plot(iterations, fitnesses)
  plot(x = iterations, y = fitnesses,
    Geom.point, Geom.line,
    Guide.xlabel("Iteration"), Guide.ylabel("Fitness"),
    Stat.step,
    #Scale.y_log10,
  ) |> drawing(10inch, 6inch)
end

function main(window)
  push!(window.assets, "widgets")

  # If we don't already have a fitness history we create one and start a
  # random walker that updates it.
  if !haskey(task_local_storage(), :fhistory)
    println("Starting random walker...")
    fitness_history = Input( (Int, Float64)[] )

    function report_new_min(i::Int, newmin::Float64)
      push!(fitness_history, push!(value(fitness_history), (i, newmin)))
    end

    # Create random walker, start task for walking
    rw = RandWalker(report_new_min)
    t = @async rand_walk!(rw)

    task_local_storage(:fhistory, fitness_history)
    task_local_storage(:rwtask, t)
  end

  local content

  lift(task_local_storage(:fhistory)) do hist
    println("Changed! len = ", length(hist))

    if length(hist) >= 1
      its  = map(first, hist)
      fits = map(t -> t[2], hist)
      content = vbox(
        hbox("Best fitness:       ", hskip(1em), @sprintf("%.3f", fits[end]) |> emph),
        hbox("Found at iteration: ", hskip(1em), string(its[end]) |> emph),
        fitness_plot(its, fits),
      ) |> packacross(center)
    else
      content = hbox("No min found yet...")
    end

    task_started = task_running = false
    if haskey(task_local_storage(), :rwtask)
      t = task_local_storage(:rwtask)
      task_running = !istaskdone(t)
      #task_started = istaskstarted(t) # Only in 0.4
    end

    vbox(
      h1("Fitness progress") |> emph,
      content,
      #hbox("Task started: ", string(task_started)), # Only in 0.4
      #hbox("Task running: ", hskip(1em), string(task_running)),
    ) |> packacross(center)
  end
end


text:
-(x::Single32) = Core.Intrinsics.neg_float(x)

+(x::Single32, y::Single32) = Core.Intrinsics.add_float(x, y)
-(x::Single32, y::Single32) = Core.Intrinsics.sub_float(x, y)
*(x::Single32, y::Single32) = Core.Intrinsics.mul_float(x, y)
/(x::Single32, y::Single32) = Core.Intrinsics.div_float(x, y)

muladd(x::Single32, y::Single32, z::Single32) = muladd_float(x, y, z)

+(x::Single32, y::Float64) = Core.Intrinsics.add_float(x, Single32(y))
-(x::Single32, y::Float64) = Core.Intrinsics.sub_float(x, Single32(y))
*(x::Single32, y::Float64) = Core.Intrinsics.mul_float(x, Single32(y))
/(x::Single32, y::Float64) = Core.Intrinsics.div_float(x, Single32(y))

+(x::Single32, y::Float32) = Core.Intrinsics.add_float(x, Single32(y))
-(x::Single32, y::Float32) = Core.Intrinsics.sub_float(x, Single32(y))
*(x::Single32, y::Float32) = Core.Intrinsics.mul_float(x, Single32(y))
/(x::Single32, y::Float32) = Core.Intrinsics.div_float(x, Single32(y))

+(x::Float64, y::Single32) = Core.Intrinsics.add_float(Single32(x), y)
-(x::Float64, y::Single32) = Core.Intrinsics.sub_float(Single32(x), y)
*(x::Float64, y::Single32) = Core.Intrinsics.mul_float(Single32(x), y)
/(x::Float64, y::Single32) = Core.Intrinsics.div_float(Single32(x), y)

+(x::Float32, y::Single32) = Core.Intrinsics.add_float(Single32(x), y)
-(x::Float32, y::Single32) = Core.Intrinsics.sub_float(Single32(x), y)
*(x::Float32, y::Single32) = Core.Intrinsics.mul_float(Single32(x), y)
/(x::Float32, y::Single32) = Core.Intrinsics.div_float(Single32(x), y)

muladd(x::Single32, y::Single32, z::Float64) = muladd_float(x, y, Single32(z))
muladd(x::Single32, y::Float64, z::Single32) = muladd_float(x, Single32(y), z)
muladd(x::Float64, y::Single32, z::Single32) = muladd_float(Single32(x), y, z)

muladd(x::Single32, y::Float64, z::Float64) = muladd_float(x, Single32(y), Single32(z))
muladd(x::Float64, y::Single32, z::Float64) = muladd_float(Single32(x), y, Single32(z))
muladd(x::Float64, y::Float64, z::Single32) = muladd_float(Single32(x), Single32(y), z)

muladd(x::Single32, y::Single32, z::Float32) = muladd_float(x, y, Single32(z))
muladd(x::Single32, y::Float32, z::Single32) = muladd_float(x, Single32(y), z)
muladd(x::Float32, y::Single32, z::Single32) = muladd_float(Single32(x), y, z)

muladd(x::Single32, y::Float32, z::Float32) = muladd_float(x, Single32(y), Single32(z))
muladd(x::Float32, y::Single32, z::Float32) = muladd_float(Single32(x), y, Single32(z))
muladd(x::Float32, y::Float32, z::Single32) = muladd_float(Single32(x), Single32(y), z)

# unary
for Op in (:+, :cbrt, :rad2deg, :deg2rad, :mod2pi, :rem2pi, :abs, :abs2, :sign)
    @eval $Op(x::Single32) = Single32($Op(Float64(x)))
end
for Op in (:exp2, :exp10, :expm1, :log2, :log10, :log1p)
    @eval $Op(x::Single32) = Single32($Op(Float64(x)))
end
# binary
for Op in (:\, :^, :div, :rem, :fld, :mod, :cld, :hypot, :min, :max, :minmax, :divrem, :fldmod, :copysign, :flipsign)
    @eval $Op(x::Single32, y::Single32) = Single32($Op(Float64(x), Float64(y)))
end
# trinary
for Op in (:clamp,)
    @eval $Op(x::Single32, y::Single32, z::Single32) = Single32($Op(Float64(x), Float64(y), Float64(z)))
end

Base.fma(x::Single32, y::Single32, z::Single32) = Single32(fma(Float64(x), Float64(y), Float64(z)))
Base.fma(x::Single32, y::Single32, z::T) where {T<:IEEEFloat} = Single32(fma(Float64(x), Float64(y), Float64(z)))
Base.fma(x::Single32, y::T, z::Single32) where {T<:IEEEFloat} = Single32(fma(Float64(x), Float64(y), Float64(z)))
Base.fma(x::T, y::Single32, z::Single32) where {T<:IEEEFloat} = Single32(fma(Float64(x), Float64(y), Float64(z)))
Base.fma(x::Single32, y::T, z::T) where {T<:IEEEFloat} = Single32(fma(Float64(x), Float64(y), Float64(z)))
Base.fma(x::T, y::Single32, z::T) where {T<:IEEEFloat} = Single32(fma(Float64(x), Float64(y), Float64(z)))
Base.fma(x::T, y::T, z::Single32) where {T<:IEEEFloat} = Single32(fma(Float64(x), Float64(y), Float64(z)))

for Op in (:sqrt, :exp, :log,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth)
    @eval $Op(x::Single32) = Single32($Op(Float64(x)))
end


text:
"""
"""
function getlabelstring2intmap(
        levels::AbstractVector,
        index::Integer,
        )
    result = Dict()
    if length(levels) == 0
        error("length(levels) == 0")
    end
    for i = 1:length(levels)
        result[levels[i]] = i - 1 + index
    end
    result = fix_type(result)
    return result
end

"""
"""
function getlabelint2stringmap(
        levels::AbstractVector,
        index::Integer,
        )
    result = Dict()
    if length(levels) == 0
        error("length(levels) == 0")
    end
    for i = 1:length(levels)
        result[i - 1 + index] = levels[i]
    end
    result = fix_type(result)
    return result
end


text:
@enum TensorLocation::UInt8 DRAM PMEM

"""
Location information for the input and output tensors of a node.

Parameterized by length of inputs and outputs because I want to be really hard on the Julia
compiler apparently.
"""
struct IOConfig{N, M}
    inputs::NTuple{N, TensorLocation}
    outputs::NTuple{M, TensorLocation}
end

Base.iterate(io::IOConfig, args...) = iterate(vflatten(io.inputs, io.outputs), args...)
Base.length(io::IOConfig{N,M}) where {N,M} = N + M

"""
$(SIGNATURES)

Return the `idx` entry of `io`, beginning with the inputs and ending with the outputs.
"""
function Base.getindex(io::IOConfig{N,M}, idx::Integer) where {N,M}
    if 1 <= idx <= N
        return io.inputs[idx]
    elseif idx <= length(io)
        return io.outputs[idx - N]
    else
        throw(BoundsError(io, idx))
    end
end
Base.getindex(io::IOConfig, I::Vector{Int}) = [io[i] for i in I]

"""
$(SIGNATURES)

Return a copy of `io` with the `idx` entry changed to `x`.
"""
function Base.setindex(io::IOConfig{N,M}, idx::Integer, x::TensorLocation) where {N,M}
    if idx <= N
        inputs = ntuple(i -> i == idx ? x : io.inputs[i], N)
        outputs = io.outputs
    elseif idx <= length(io)
        idx = idx - length(io.inputs)
        inputs = io.inputs
        outputs = ntuple(i -> i == idx ? x : io.outputs[i], M)
    end
    return IOConfig{N,M}(inputs, outputs)
end

function Base.isless(a::IOConfig{N,M}, b::IOConfig{N,M}) where {N,M}
    return (a.inputs < b.inputs) || ((a.inputs == b.inputs) && a.outputs < b.outputs)
end

function Base.show(io::IO, config::IOConfig{N,M}) where {N,M}
    f = x -> (x == DRAM) ? "DRAM" : "PMEM"
    print(io, "IOConfig{$N,$M}: ")
    print(io, "(", join(f.(config.inputs), ", "), ") -- ")
    print(io, "(", join(f.(config.outputs), ", "), ")")
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
export  initPopReport,
        initPopAgeReport

"""
```
initPopReport( mpSim::MPsim )
```
This function creates a composition report on the initial population of the manpower simulation `mpSim`.

This function returns a `DataFrame`, with the nodes to which the members of the initial belong and the number in each node. If there is no population, this function returns an empty `DataFrame`.
"""
function initPopReport( mpSim::MPsim )::DataFrame

    queryCmd = string( "SELECT targetNode, count( targetNode ) amount FROM `",
        mpSim.transDBname , "`",
        "\n    WHERE transition IS 'Init'",
        "\n    GROUP BY targetNode" )
    return DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )

end  # function initPopReport( mpSim )

"""
```
initPopAgeReport(
    mpSim::MPsim,
    ageRes::Real,
    ageType::Symbol )
```
This function creates an age distribution report on the initial population of the manpower simulation `mpSim`. The report is generated on a grid of ages with resoluion `ageRes`.

The type of age that is reported is determined by `ageType`, which can take these valid values:
* `:age` for the actual age of the personnel members,
* `:tenure` for the tenure of the personnel members, and
* `:timeInNode` for the time the personnel members are in their current base node.

This function returns a `DataFrame`, with the number of personnel members in each age bracket. The last 5 entries in the `DataFrame` are the summary statistics (mean, median, standard deviation, minimum, and maximum) of the ages.
"""
function initPopAgeReport( mpSim::MPsim, ageRes::Real,
    ageType::Symbol )::DataFrame

    if ageRes <= 0
        @warn "Age resolution must be > 0, cannot generate report."
        return DataFrame()
    end  # if ageRes <= 0

    if ageType ∉ [:age, :tenure, :timeInNode]
        @warn "Unknown type of data requested, cannot generate report."
        return DataFrame()
    end  # if ageType ∉ [:age, :tenure, :timeInNode]

    subpopulation = Subpopulation( "Init" )
    addSubpopulationCondition!( subpopulation,
        MPcondition( "had transition", ==, "Init" ) )
    ages = getAgesAtTime( mpSim, 0.0, [subpopulation], ageType )[1]

    if isempty( ages )
        return DataFrame()
    end  # if isempty( ages )

    return processInitPopAges( ages, Float64( ageRes ) )

end  # initPopAgeReport( mpSim, ageRes, ageType )

include( joinpath( repPrivPath, "initpop.jl" ) )


text:
function savedataset!(path, storage::FileStorage, block::Image{2}, data)
    ensuredir(storage, path)

    # Create saving jobs
    npad = ceil(Int, log10(nobs(data)))
    savejobs = mapobs(1:nobs(data)) do i
        img = getobs(data, i)
        ext = get(storage.preferences, :imageformat, "jpg")
        filename = "$(lpad(string(i), npad, "0")).$(ext)"
        save(string(joinpath(path, filename)), img)
        return "Done"
    end

    # Run export
    progress = Progress(
        nobs(data),
        enabled = storage.verbose,
        desc = "Writing images: ",
        showspeed = true,
    )
    iter = if get(storage.preferences, :parallel, true)
        data -> eachobsparallel(data, buffered = false)
    else
        eachobs
    end
    for res in iter(savejobs)
        next!(progress)
    end
end

function loaddataset(path, ::FileStorage, ::Image{2})
    return mapobs(loadfile, FileDataset(path))
end


text:
# Result type.

"""
A `Result{T}` is returned from every API function.
It encapsulates the value, HTTP response, and thrown exception of the call.
"""
struct Result{T}
    val::Union{T, Nothing}
    resp::Union{HTTP.Response, Nothing}
    ex::Union{Exception, Nothing}
    st::StackTraces.StackTrace
end

Result(val::T, resp::HTTP.Response) where T = Result{T}(val, resp, nothing, [])
Result{T}(
    e::Exception,
    bt::Vector{Union{Ptr{Nothing}, Base.InterpreterIP}},
    resp::Union{HTTP.Response, Nothing}=nothing,
) where T = Result{T}(nothing, resp, e, stacktrace(bt))

"""
    value(::Result{T}) -> Union{T, Nothing}

Returns the result's value, if any exists.
"""
value(r::Result) = r.val

"""
    response(::Result) -> Union{HTTP.Response, Nothing}

Returns the result's HTTP response, if any exists.
"""
response(r::Result) = r.resp

"""
    exception(::Result{T}) -> Union{Tuple{Exception, Vector{StackFrame}}, Nothing}

Returns the result's thrown exception and stack trace, if any exists.
"""
exception(r::Result) = r.ex === nothing ? nothing : (r.ex, r.st)

# Response postprocessing.

"""
Determines the behaviour of [`postprocess`](@ref).
"""
abstract type PostProcessor end

"""
Does nothing and always returns `nothing`.
"""
struct DoNothing <: PostProcessor end

"""
    JSON(f::Function=identity) -> JSON

Parses a JSON response into a given type and runs `f` on that object.
"""
struct JSON <: PostProcessor
    f::Function
end
JSON() = JSON(identity)

"""
    DoSomething(::Function) -> DoSomething

Runs a user-defined postprocessor.
"""
struct DoSomething <: PostProcessor
    f::Function
end

"""
    postprocess(::PostProcessor, ::HTTP.Response, ::Type{T})

Computes a value from an HTTP response.
This is what is returned by [`value`](@ref).
"""
postprocess(::DoNothing, ::HTTP.Response, ::Type) = nothing
postprocess(p::JSON, r::HTTP.Response, ::Type{T}) where T =
    p.f(JSON2.read(IOBuffer(r.body), T))
postprocess(p::DoSomething, r::HTTP.Response, ::Type) = p.f(r)

# Requests.

"""
    request(
        f::Forge, fun::Function, ep::Endpoint;
        headers::Vector{<:Pair}=HTTP.Header[],
        query::AbstractDict=Dict(),
        request_opts=Dict(),
        kwargs...,
    ) -> Result{T}

Make an HTTP request and return a [`Result`](@ref).
`T` is determined by [`into`](@ref).

## Arguments
- `f::Forge` A [`Forge`](@ref) subtype.
- `fun::Function`: The API function being called.
- `ep::Endpoint`: The endpoint information.

## Keywords
- `query::AbstractDict=Dict()`: Query string parameters to add to the request.
- `headers::Vector{<:Pair}=HTTP.Header[]`: Headers to add to the request.
- `request_opts=Dict()`: Keywords passed into `HTTP.request`.

Trailing keywords are sent as a JSON body for `PATCH`, `POST`, and `PUT` requests.
For other request types, the keywords are sent as query string parameters.

!!! note
    Every API function passes its keyword arguments into this function.
    Therefore, to customize behaviour for a single request, pass the above keywords to the API function.
"""
function request(
    f::Forge, fun::Function, ep::Endpoint;
    headers::Vector{<:Pair}=HTTP.Header[],
    query::AbstractDict=Dict(),
    request_opts=Dict(),
    kwargs...,
)
    T = into(f, fun)

    if has_rate_limits(f, fun) && rate_limit_check(f, fun)
        orl = on_rate_limit(f, fun)
        if orl === ORL_RETURN
            return Result{T}(RateLimited(rate_limit_period(f, fun)), backtrace())
        elseif orl === ORL_WAIT
            rate_limit_wait(f, fun)
        else
            @warn "Ignoring unknown rate limit behaviour $orl"
        end
    end

    url = base_url(f) * ep.url
    headers = vcat(request_headers(f, fun), ep.headers, headers)
    query = merge(request_query(f, fun), ep.query, query)
    opts = merge(request_kwargs(f, fun), Dict(pairs(request_opts)))
    body = if ep.method in (:PATCH, :POST, :PUT)
        JSON2.write(Dict(kwargs))
    else
        merge!(query, Dict(kwargs))
        HTTP.nobody
    end

    resp = try
        HTTP.request(
            ep.method, url, headers, body;
            # Never throw status exceptions, we'll handle the status ourselves.
            query=query, opts..., status_exception=false,
        )
    catch e
        return Result{T}(e, catch_backtrace())
    end

    has_rate_limits(f, fun) && rate_limit_update!(f, fun, resp)

    resp.status >= 300 && !(resp.status == 404 && ep.allow_404) &&
        return Result{T}(HTTP.StatusError(resp.status, resp), backtrace(), resp)

    val = try
        postprocess(postprocessor(f, fun), resp, T)
    catch e
        return Result{T}(e, catch_backtrace(), resp)
    end

    return Result(val, resp)
end


text:
@testset "CPLayerGraph()" begin
    trailer = SeaPearl.Trailer()
    model = SeaPearl.CPModel(trailer)

    x = SeaPearl.IntVar(2, 3, "x", trailer)
    y = SeaPearl.IntVar(2, 3, "y", trailer)
    z = SeaPearl.IntVarViewOpposite(y, "-y")
    SeaPearl.addVariable!(model, x)
    SeaPearl.addVariable!(model, y)
    SeaPearl.addConstraint!(model, SeaPearl.Equal(x, y, trailer))
    SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x, y, trailer))
    SeaPearl.addConstraint!(model, SeaPearl.NotEqual(z, y, trailer))

    g = SeaPearl.CPLayerGraph(model)

    true_idToNode = [
        SeaPearl.ConstraintVertex(model.constraints[1]),
        SeaPearl.ConstraintVertex(model.constraints[2]),
        SeaPearl.ConstraintVertex(model.constraints[3]),
        SeaPearl.VariableVertex(z),
        SeaPearl.VariableVertex(x),
        SeaPearl.VariableVertex(y),
        SeaPearl.ValueVertex(2),
        SeaPearl.ValueVertex(3)
    ]

    @test g.cpmodel == model
    @test g.idToNode == true_idToNode
    @test length(keys(g.nodeToId)) == 8

    for i in 1:7
        @test g.nodeToId[true_idToNode[i]] == i
    end

    @test Matrix(LightGraphs.LinAlg.adjacency_matrix(g.fixedEdgesGraph)) == [
        0 0 0 0 1 1
        0 0 0 0 1 1
        0 0 0 1 0 1
        0 0 1 0 0 1
        1 1 0 0 0 0
        1 1 1 1 0 0
    ]

    @test g.numberOfConstraints == 3
    @test g.numberOfVariables == 3
    @test g.numberOfValues == 2
    @test g.totalLength == 8

    empty_g = SeaPearl.CPLayerGraph()
    @test isnothing(empty_g.cpmodel)
end


text:
@testset "Test simple inputs" begin

    parameters1 = BAT.NamedTupleDist(
        p1 = -20..20,
        p2 = -10..10,
    )

    function testfunc1(params)
        c = [20.12, 5.56, 325.556]
        return c[1] * params.p1^2 + c[2] * params.p1 * params.p2 + c[3] * params.p2^2
    end

    measurements1 = (
        meas1 = Measurement(testfunc1, 111.1,
        uncertainties = (unc1=10.1, unc2=12.2, unc3=13.3), active=true),

        meas2 = Measurement(Observable(testfunc1, min=0, max=1000), 222.2,
        uncertainties = (unc1=20.1, unc2=20.2, unc3=23.3), active=false),

        meas3 = Measurement(Observable(testfunc1, min=0, max=1000), 333.3,
        uncertainties = (unc1=30.1, unc2=30.2, unc3=30.3), active=true),
    )

    correlations1 = (
        unc1 = NoCorrelation(active=true),

        unc2 = Correlation([1.0 0.5 0.7;
                            0.5 1.0 0.6;
                            0.7 0.6 1.0], active=false),

        unc3 = Correlation([1.0 0.2 0.3;
                            0.2 1.0 0.4;
                            0.3 0.4 1.0])
    )

    model1 = EFTfitterModel(parameters1, measurements1, correlations1)

    @testset "Test EFTfitterModel" begin
        @test model1.measurements.meas1.observable == Observable(testfunc1, -Inf, Inf)

        @test length(model1.measurements) == 2
        @test keys(model1.measurements) == (:meas1, :meas3)
        @test model1.measurements.meas1 == Measurement(Observable(testfunc1, -Inf, Inf), 111.1, (unc1=10.1, unc3=13.3), true)
        @test model1.measurements.meas3 == Measurement(Observable(testfunc1, min=0, max=1000), 333.3, (unc1=30.1, unc3=30.3), true)

        @test get_observables(model1) == (testfunc1 = Observable(testfunc1, -Inf, Inf),)
    end

    @testset "Test EFTfitterDensity" begin
        eftfitter_density = PosteriorDensity(model1).likelihood
        @test isa(eftfitter_density, EFTfitterDensity)

        @test eftfitter_density.measured_values == [111.1, 333.3]
        @test eftfitter_density.observable_functions == [testfunc1, testfunc1]
        @test eftfitter_density.observable_mins == [-Inf, 0]
        @test eftfitter_density.observable_maxs == [Inf, 1000.0]
        @test eftfitter_density.invcov ≈ [ 0.00369157   -0.000244669; -0.000244669   0.000564432] rtol=0.01
        @test eftfitter_density.check_bounds == true

        # test evaluation of correct likelihood value at a few points:
        @test BAT.eval_logval_unchecked(eftfitter_density, (p1=0.0, p2=0.0)) ≈ -45.07398 rtol=0.01
        @test BAT.eval_logval_unchecked(eftfitter_density, (p1=10.0, p2=-30.0)) ≈ -Inf
        @test BAT.eval_logval_unchecked(eftfitter_density, (p1=1.5, p2=-0.9)) ≈ -68.6575 rtol=0.01
        @test BAT.eval_logval_unchecked(eftfitter_density, (p1=-1.2, p2=-0.8)) ≈ -37.1857 rtol=0.01
    end

end


text:
#!/home/dave/.local/bin/julia --startup-file=no
readline()
println("id name davidsguterbot")
println("id author David Weingut")
println("uciok")
using Pkg
Pkg.activate(".")
using LichessBot
LichessBot.uci_listen(;id=false)


text:
using Sobol

using ResumableFunctions

@resumable function sobols(n)
    it = SobolSeq(n)
    while true
        for p in next!(it)
            @yield p
        end
    end
end