text:
# Formatted Output (B.17)

export @cuprintf

@generated function promote_c_argument(arg)
    # > When a function with a variable-length argument list is called, the variable
    # > arguments are passed using C's old ``default argument promotions.'' These say that
    # > types char and short int are automatically promoted to int, and type float is
    # > automatically promoted to double. Therefore, varargs functions will never receive
    # > arguments of type char, short int, or float.

    if arg == Cchar || arg == Cshort || arg == Cuchar || arg == Cushort
        return :(Cint(arg))
    elseif arg == Cfloat
        return :(Cdouble(arg))
    else
        return :(arg)
    end
end

"""
    @cuprintf("%Fmt", args...)

Print a formatted string in device context on the host standard output.

Note that this is not a fully C-compliant `printf` implementation; see the CUDA
documentation for supported options and inputs.

Also beware that it is an untyped, and unforgiving `printf` implementation. Type widths need
to match, eg. printing a 64-bit Julia integer requires the `%ld` formatting string.
"""
macro cuprintf(fmt::String, args...)
    fmt_val = Val(Symbol(fmt))

    return :(_cuprintf($fmt_val, $(map(arg -> :(promote_c_argument($arg)), esc.(args))...)))
end

@generated function _cuprintf(::Val{fmt}, argspec...) where {fmt}
    Context() do ctx
        arg_exprs = [:( argspec[$i] ) for i in 1:length(argspec)]
        arg_types = [argspec...]

        T_void = LLVM.VoidType(ctx)
        T_int32 = LLVM.Int32Type(ctx)
        T_pint8 = LLVM.PointerType(LLVM.Int8Type(ctx))

        # create functions
        param_types = LLVMType[convert(LLVMType, typ; ctx) for typ in arg_types]
        llvm_f, _ = create_function(T_int32, param_types)
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            str = globalstring_ptr!(builder, String(fmt))

            # construct and fill args buffer
            if isempty(argspec)
                buffer = LLVM.PointerNull(T_pint8)
            else
                argtypes = LLVM.StructType("printf_args"; ctx)
                elements!(argtypes, param_types)

                args = alloca!(builder, argtypes)
                for (i, param) in enumerate(parameters(llvm_f))
                    p = struct_gep!(builder, args, i-1)
                    store!(builder, param, p)
                end

                buffer = bitcast!(builder, args, T_pint8)
            end

            # invoke vprintf and return
            vprintf_typ = LLVM.FunctionType(T_int32, [T_pint8, T_pint8])
            vprintf = LLVM.Function(mod, "vprintf", vprintf_typ)
            chars = call!(builder, vprintf, [str, buffer])

            ret!(builder, chars)
        end

        call_function(llvm_f, Int32, Tuple{arg_types...}, arg_exprs...)
    end
end

## print-like functionality

export @cuprint, @cuprintln

# simple conversions, defining an expression and the resulting argument type. nothing fancy,
# `@cuprint` pretty directly maps to `@cuprintf`; we should just support `write(::IO)`.
const cuprint_conversions = Dict(
    Float32         => (x->:(Float64($x)),                  Float64),
    Ptr{<:Any}      => (x->:(convert(Ptr{Cvoid}, $x)),      Ptr{Cvoid}),
    LLVMPtr{<:Any}  => (x->:(reinterpret(Ptr{Cvoid}, $x)),  Ptr{Cvoid}),
    Bool            => (x->:(Int32($x)),                    Int32),
)

# format specifiers
const cuprint_specifiers = Dict(
    # integers
    Int16       => "%hd",
    Int32       => "%d",
    Int64       => Sys.iswindows() ? "%lld" : "%ld",
    UInt16      => "%hu",
    UInt32      => "%u",
    UInt64      => Sys.iswindows() ? "%llu" : "%lu",

    # floating-point
    Float64     => "%f",

    # other
    Cchar       => "%c",
    Ptr{Cvoid}  => "%p",
    Cstring     => "%s",
)

@generated function _cuprint(parts...)
    fmt = ""
    args = Expr[]

    for i in 1:length(parts)
        part = :(parts[$i])
        T = parts[i]

        # put literals directly in the format string
        if T <: Val
            fmt *= string(T.parameters[1])
            continue
        end

        # try to convert arguments if they are not supported directly
        if !haskey(cuprint_specifiers, T)
            for Tmatch in keys(cuprint_conversions)
                if T <: Tmatch
                    conv, T = cuprint_conversions[Tmatch]
                    part = conv(part)
                    break
                end
            end
        end

        # render the argument
        if haskey(cuprint_specifiers, T)
            fmt *= cuprint_specifiers[T]
            push!(args, part)
        elseif T <: Tuple
            fmt *= "("
            for (j, U) in enumerate(T.parameters)
                if haskey(cuprint_specifiers, U)
                    fmt *= cuprint_specifiers[U]
                    push!(args, :($part[$j]))
                    if j < length(T.parameters)
                        fmt *= ", "
                    elseif length(T.parameters) == 1
                        fmt *= ","
                    end
                else
                    @error("@cuprint does not support values of type $U")
                end
            end
            fmt *= ")"
        elseif T <: String
            @error("@cuprint does not support non-literal strings")
        else
            @error("@cuprint does not support values of type $T")
        end
    end

    quote
        Base.@_inline_meta
        @cuprintf($fmt, $(args...))
    end
end

"""
    @cuprint(xs...)
    @cuprintln(xs...)

Print a textual representation of values `xs` to standard output from the GPU. The
functionality builds on `@cuprintf`, and is intended as a more use friendly alternative of
that API. However, that also means there's only limited support for argument types, handling
16/32/64 signed and unsigned integers, 32 and 64-bit floating point numbers, `Cchar`s and
pointers. For more complex output, use `@cuprintf` directly.

Limited string interpolation is also possible:

```julia
    @cuprint("Hello, World ", 42, "\\n")
    @cuprint "Hello, World \$(42)\\n"
```
"""
macro cuprint(parts...)
    args = Union{Val,Expr,Symbol}[]

    parts = [parts...]
    while true
        isempty(parts) && break

        part = popfirst!(parts)

        # handle string interpolation
        if isa(part, Expr) && part.head == :string
            parts = vcat(part.args, parts)
            continue
        end

        # expose literals to the generator by using Val types
        if isbits(part) # literal numbers, etc
            push!(args, Val(part))
        elseif isa(part, QuoteNode) # literal symbols
            push!(args, Val(part.value))
        elseif isa(part, String) # literal strings need to be interned
            push!(args, Val(Symbol(part)))
        else # actual values that will be passed to printf
            push!(args, part)
        end
    end

    quote
        _cuprint($(map(esc, args)...))
    end
end

@doc (@doc @cuprint) ->
macro cuprintln(parts...)
    esc(quote
        CUDA.@cuprint($(parts...), "\n")
    end)
end

export @cushow

"""
    @cushow(ex)

GPU analog of `Base.@show`. It comes with the same type restrictions as [`@cuprintf`](@ref).

```julia
@cushow threadIdx().x
```
"""
macro cushow(exs...)
    blk = Expr(:block)
    for ex in exs
        push!(blk.args, :(CUDA.@cuprintln($(sprint(Base.show_unquoted,ex)*" = "),
                                          begin local value = $(esc(ex)) end)))
    end
    isempty(exs) || push!(blk.args, :value)
    blk
end


text:
configs = Dict()
for (i, tick_authority_policy) in enumerate(0:20:300)
    configs[i] = EpidemicModel.Config(
        n_nodes=300,
        n0=3,
        k=3,
        agentcount=1000,
        transmission_prob=0.05,
        test_prob=0.8,
        quarantine_duration=140,
        behavior_reduction_factor=0.25,
        behavior_activated=true,
        tick_authority_recommendation=1001,
        tick_authority_policy=tick_authority_policy,
        n_public_spaces=100,
        initial_infectives=5,
        n_iter=1000
    )
end

ext: jl


text:
# This file is a part of JuliaPackageTemplate.jl, licensed under the MIT License (MIT).

__precompile__(true)

"""
    JuliaPackageTemplate

Template for Julia packages.
"""
module JuliaPackageTemplate

include("hello_world.jl")

end # module


text:
@testset "Is integral" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 1
  A = AlgAss(f)

  @test Hecke.isintegral(A[1]) == true
  @test Hecke.isintegral(fmpq(1, 2)*A[1]) == false

  B = AlgGrp(FlintQQ, small_group(2, 1))
  @test Hecke.isintegral(B[1]) == true
  @test Hecke.isintegral(fmpq(1, 2)*B[1]) == false

  C = AlgMat(FlintQQ, B, 2)
  @test Hecke.isintegral(C[1]) == true
  @test Hecke.isintegral(fmpq(1, 2)*C[1]) == false
end


text:
function CachedOperator(::Type{BandedMatrix},op::Operator;padding::Bool=false)
    l,u=bandwidths(op)
    padding && (u+=l)
    data = BandedMatrix{eltype(op)}(undef, (0,0), (l,u))
    CachedOperator(op,data,size(data),domainspace(op),rangespace(op),(-l,u),padding)
end

## Grow cached operator

function resizedata!(B::CachedOperator{T,<:BandedMatrix{T}},n::Integer,m_in::Integer) where T<:Number
    m = max(m_in,n+B.data.u)
    N,M = size(B)
    n = min(n, N)

    if n > B.datasize[1]
        pad!(B.data,min(N,2n),m)

        kr=B.datasize[1]+1:n
        jr=max(B.datasize[1]+1-B.data.l,1):min(n+B.data.u,M)
        BLAS.axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))

        B.datasize = (n,m)
    end

    B
end

resizedata!(B::CachedOperator{T,<:BandedMatrix{T}},n::Integer,::Colon) where T<:Number =
    resizedata!(B, n, n+B.data.u)

## Grow QR

function QROperator(R::CachedOperator{T,<:BandedMatrix{T}}) where T
    M = R.data.l+1   # number of diag+subdiagonal bands
    H = Array{T}(undef,M,100)
    QROperator(R,H,0)
end

function resizedata!(QR::QROperator{<:CachedOperator{T,<:BandedMatrix{T},
                                                  MM,DS,RS,BI}},
         ::Colon,col) where {T,MM,DS,RS,BI}
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R_cache
    W=QR.H

    R=MO.data
    M=R.l+1   # number of diag+subdiagonal bands

    if col+M-1 ≥ MO.datasize[1]
        resizedata!(MO,(col+M-1)+100,:)  # double the last rows
    end

    if col > size(W,2)
        W=QR.H=unsafe_resize!(W,:,2col)
    end

    for k=QR.ncols+1:col
        W[:,k] = view(R.data,R.u+1:R.u+R.l+1,k) # diagonal and below
        wp=view(W,:,k)
        W[1,k]+= flipsign(norm(wp),W[1,k])
        normalize!(wp)

        # scale banded entries
        for j=k:k+R.u
            dind=R.u+1+k-j
            v=view(R.data,dind:dind+M-1,j)
            dt=dot(wp,v)
            LinearAlgebra.axpy!(-2*dt,wp,v)
        end

        # scale banded/filled entries
        for j=k+R.u+1:k+R.u+M-1
            p=j-k-R.u
            v=view(R.data,1:M-p,j)  # shift down each time
            wp2=view(wp,p+1:M)
            dt=dot(wp2,v)
            LinearAlgebra.axpy!(-2*dt,wp2,v)
        end
    end
    QR.ncols=col
    QR
end

# BLAS versions, requires BlasFloat

function resizedata!(QR::QROperator{<:CachedOperator{T,<:BandedMatrix{T},
                                       MM,DS,RS,BI}},
                     ::Colon,col) where {T<:BlasFloat,MM,DS,RS,BI}
    if col ≤ QR.ncols
        return QR
    end

    col = min(col, size(QR,2))

    MO=QR.R_cache
    W=QR.H

    R=MO.data
    M=R.l+1   # number of diag+subdiagonal bands

    if col+M-1 ≥ MO.datasize[1]
        resizedata!(MO,(col+M-1)+100,:)  # double the last rows
    end

    if col > size(W,2)
        W=QR.H=unsafe_resize!(W,:,2col)
    end

    m,n=size(R)
    w=pointer(W)
    r=pointer(R.data)
    sz=sizeof(T)
    st=stride(R.data,2)
    stw=stride(W,2)

    for k=QR.ncols+1:col
        v=r+sz*(R.u + (k-1)*st)    # diagonal entry
        wp=w+stw*sz*(k-1)          # k-th column of W
        BLAS.blascopy!(M,v,1,wp,1)
        W[1,k]+= flipsign(BLAS.nrm2(M,wp,1),W[1,k])
        normalize!(M,wp)

        for j=k:k+R.u
            v=r+sz*(R.u + (k-1)*st + (j-k)*(st-1))
            dt = BandedMatrices.dot(M,wp,1,v,1)
            BLAS.axpy!(M,-2*dt,wp,1,v,1)
        end

        for j=k+R.u+1:k+R.u+M-1
            p=j-k-R.u
            v=r+sz*((j-1)*st)  # shift down each time
            dt = BandedMatrices.dot(M-p,wp+p*sz,1,v,1)
            BLAS.axpy!(M-p,-2*dt,wp+p*sz,1,v,1)
        end
    end
    QR.ncols=col
    QR
end

## back substitution
# loop to avoid ambiguity with AbstractTRiangular
for ArrTyp in (:AbstractVector, :AbstractMatrix)
    @eval function ldiv!(U::UpperTriangular{T, <:SubArray{T, 2, <:BandedMatrix{T}, Tuple{UnitRange{Int}, UnitRange{Int}}, false}},
                             u::$ArrTyp{T}) where T
        n = size(u,1)
        n == size(U,1) || throw(DimensionMismatch())

        V = parent(U)
        @assert first(parentindices(V)[1]) == 1
        @assert first(parentindices(V)[2]) == 1

        A = parent(V)

        b=bandwidth(A,2)

        for c=1:size(u,2)
            for k=n:-1:1
                @simd for j=k+1:min(n,k+b)
                    @inbounds u[k,c] = muladd(-A.data[k-j+A.u+1,j],u[j,c],u[k,c])
                end

                @inbounds u[k,c] /= A.data[A.u+1,k]
            end
        end
        u
    end
end


text:
# Tests for Profit General Direct Approach DEA Models
@testset "ProfitGDADEAAModel" begin

    # Test using Book data
    X = [2; 4; 8; 12; 6; 14; 14; 9.412];
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353];
    W = [1; 1; 1; 1; 1; 1; 1; 1];
    P = [2; 2; 2; 2; 2; 2; 2; 2];

    # Profit GDA
    profitgda = deaprofitgda(X, Y, W, P, :ERG)

    @test typeof(profitgda) == ProfitGDADEAModel

    @test nobs(profitgda) == 8
    @test ninputs(profitgda) == 1
    @test noutputs(profitgda) == 1
    @test ismonetary(profitgda) == false

    @test profitgda.measure == :ERG
    @test efficiency(profitgda) == efficiency(profitgda, :Economic)
    @test efficiency(profitgda, :Economic)   ≈ [4.0; 0.5; 0.0; 1/6; 0.8; 4/7; 2/7; 0.949449] atol = 1e-5
    @test efficiency(profitgda, :Technical)  ≈ [0.0; 0.0; 0.0; 0.0; 0.6; 11/21; 1/7; 0.8] atol = 1e-5
    @test efficiency(profitgda, :Allocative) ≈ [4.0; 0.5; 0.0; 1/6; 0.2; 1/21; 1/7; 0.149449] atol = 1e-5
    @test normfactor(profitgda) ≈ [2.0; 4.0; 8.0; 12.0; 10.0; 14.0; 14.0; 13.3825]

    # With monetary = true
    profitgdamon = deaprofitgda(X, Y, W, P, :ERG, monetary = true)

    @test ismonetary(profitgdamon) == true

    @test efficiency(profitgdamon) ≈ efficiency(profitgda) .* normfactor(profitgda)
    @test efficiency(profitgdamon, :Technical) ≈ efficiency(profitgda, :Technical) .* normfactor(profitgda)
    @test efficiency(profitgdamon, :Allocative) ≈ efficiency(profitgda, :Allocative) .* normfactor(profitgda)

    # Print
    show(IOBuffer(), profitgda)

    # Test errors
    @test_throws DimensionMismatch deaprofitgda([1; 2 ; 3], [4; 5], [1; 1; 1], [1; 1], :ERG) #  Different number of observations
    @test_throws DimensionMismatch deaprofitgda([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], [1; 1; 1], :ERG) # Different number of observation in input prices
    @test_throws DimensionMismatch deaprofitgda([1; 2; 3], [4; 5; 6], [1; 1; 1], [1; 2; 3; 4], :ERG) # Different number of observation in output prices
    @test_throws DimensionMismatch deaprofitgda([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], [1; 1; 1], :ERG) # Different number of input prices and inputs
    @test_throws DimensionMismatch deaprofitgda([1; 2; 3 ], [4; 5; 6], [1; 1; 1], [1 1 1; 2 2 2; 3 3 3], :ERG) # Different number of output prices and outputs
    @test_throws ArgumentError deaprofitgda([1; 2; 3], [4; 5; 6], [1; 2; 3], [4; 5; 6], :Error) # Invalid associated efficiency measure

end


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

mutable struct ExpressRouteCircuitPeeringPropertiesFormat <: SwaggerModel
    peeringType::Any # spec type: Union{ Nothing, String } # spec name: peeringType
    state::Any # spec type: Union{ Nothing, String } # spec name: state
    azureASN::Any # spec type: Union{ Nothing, Int32 } # spec name: azureASN
    peerASN::Any # spec type: Union{ Nothing, Int32 } # spec name: peerASN
    primaryPeerAddressPrefix::Any # spec type: Union{ Nothing, String } # spec name: primaryPeerAddressPrefix
    secondaryPeerAddressPrefix::Any # spec type: Union{ Nothing, String } # spec name: secondaryPeerAddressPrefix
    primaryAzurePort::Any # spec type: Union{ Nothing, String } # spec name: primaryAzurePort
    secondaryAzurePort::Any # spec type: Union{ Nothing, String } # spec name: secondaryAzurePort
    sharedKey::Any # spec type: Union{ Nothing, String } # spec name: sharedKey
    vlanId::Any # spec type: Union{ Nothing, Int32 } # spec name: vlanId
    microsoftPeeringConfig::Any # spec type: Union{ Nothing, ExpressRouteCircuitPeeringConfig } # spec name: microsoftPeeringConfig
    stats::Any # spec type: Union{ Nothing, ExpressRouteCircuitStats } # spec name: stats
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    gatewayManagerEtag::Any # spec type: Union{ Nothing, String } # spec name: gatewayManagerEtag
    lastModifiedBy::Any # spec type: Union{ Nothing, String } # spec name: lastModifiedBy
    routeFilter::Any # spec type: Union{ Nothing, RouteFilter } # spec name: routeFilter
    ipv6PeeringConfig::Any # spec type: Union{ Nothing, Ipv6ExpressRouteCircuitPeeringConfig } # spec name: ipv6PeeringConfig

    function ExpressRouteCircuitPeeringPropertiesFormat(;peeringType=nothing, state=nothing, azureASN=nothing, peerASN=nothing, primaryPeerAddressPrefix=nothing, secondaryPeerAddressPrefix=nothing, primaryAzurePort=nothing, secondaryAzurePort=nothing, sharedKey=nothing, vlanId=nothing, microsoftPeeringConfig=nothing, stats=nothing, provisioningState=nothing, gatewayManagerEtag=nothing, lastModifiedBy=nothing, routeFilter=nothing, ipv6PeeringConfig=nothing)
        o = new()
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("peeringType"), peeringType)
        setfield!(o, Symbol("peeringType"), peeringType)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("azureASN"), azureASN)
        setfield!(o, Symbol("azureASN"), azureASN)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("peerASN"), peerASN)
        setfield!(o, Symbol("peerASN"), peerASN)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("primaryPeerAddressPrefix"), primaryPeerAddressPrefix)
        setfield!(o, Symbol("primaryPeerAddressPrefix"), primaryPeerAddressPrefix)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("secondaryPeerAddressPrefix"), secondaryPeerAddressPrefix)
        setfield!(o, Symbol("secondaryPeerAddressPrefix"), secondaryPeerAddressPrefix)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("primaryAzurePort"), primaryAzurePort)
        setfield!(o, Symbol("primaryAzurePort"), primaryAzurePort)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("secondaryAzurePort"), secondaryAzurePort)
        setfield!(o, Symbol("secondaryAzurePort"), secondaryAzurePort)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("sharedKey"), sharedKey)
        setfield!(o, Symbol("sharedKey"), sharedKey)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("vlanId"), vlanId)
        setfield!(o, Symbol("vlanId"), vlanId)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("microsoftPeeringConfig"), microsoftPeeringConfig)
        setfield!(o, Symbol("microsoftPeeringConfig"), microsoftPeeringConfig)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("stats"), stats)
        setfield!(o, Symbol("stats"), stats)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("gatewayManagerEtag"), gatewayManagerEtag)
        setfield!(o, Symbol("gatewayManagerEtag"), gatewayManagerEtag)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("lastModifiedBy"), lastModifiedBy)
        setfield!(o, Symbol("lastModifiedBy"), lastModifiedBy)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("routeFilter"), routeFilter)
        setfield!(o, Symbol("routeFilter"), routeFilter)
        validate_property(ExpressRouteCircuitPeeringPropertiesFormat, Symbol("ipv6PeeringConfig"), ipv6PeeringConfig)
        setfield!(o, Symbol("ipv6PeeringConfig"), ipv6PeeringConfig)
        o
    end
end # type ExpressRouteCircuitPeeringPropertiesFormat

const _property_map_ExpressRouteCircuitPeeringPropertiesFormat = Dict{Symbol,Symbol}(Symbol("peeringType")=>Symbol("peeringType"), Symbol("state")=>Symbol("state"), Symbol("azureASN")=>Symbol("azureASN"), Symbol("peerASN")=>Symbol("peerASN"), Symbol("primaryPeerAddressPrefix")=>Symbol("primaryPeerAddressPrefix"), Symbol("secondaryPeerAddressPrefix")=>Symbol("secondaryPeerAddressPrefix"), Symbol("primaryAzurePort")=>Symbol("primaryAzurePort"), Symbol("secondaryAzurePort")=>Symbol("secondaryAzurePort"), Symbol("sharedKey")=>Symbol("sharedKey"), Symbol("vlanId")=>Symbol("vlanId"), Symbol("microsoftPeeringConfig")=>Symbol("microsoftPeeringConfig"), Symbol("stats")=>Symbol("stats"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("gatewayManagerEtag")=>Symbol("gatewayManagerEtag"), Symbol("lastModifiedBy")=>Symbol("lastModifiedBy"), Symbol("routeFilter")=>Symbol("routeFilter"), Symbol("ipv6PeeringConfig")=>Symbol("ipv6PeeringConfig"))
const _property_types_ExpressRouteCircuitPeeringPropertiesFormat = Dict{Symbol,String}(Symbol("peeringType")=>"String", Symbol("state")=>"String", Symbol("azureASN")=>"Int32", Symbol("peerASN")=>"Int32", Symbol("primaryPeerAddressPrefix")=>"String", Symbol("secondaryPeerAddressPrefix")=>"String", Symbol("primaryAzurePort")=>"String", Symbol("secondaryAzurePort")=>"String", Symbol("sharedKey")=>"String", Symbol("vlanId")=>"Int32", Symbol("microsoftPeeringConfig")=>"ExpressRouteCircuitPeeringConfig", Symbol("stats")=>"ExpressRouteCircuitStats", Symbol("provisioningState")=>"String", Symbol("gatewayManagerEtag")=>"String", Symbol("lastModifiedBy")=>"String", Symbol("routeFilter")=>"RouteFilter", Symbol("ipv6PeeringConfig")=>"Ipv6ExpressRouteCircuitPeeringConfig")
Base.propertynames(::Type{ ExpressRouteCircuitPeeringPropertiesFormat }) = collect(keys(_property_map_ExpressRouteCircuitPeeringPropertiesFormat))
Swagger.property_type(::Type{ ExpressRouteCircuitPeeringPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ExpressRouteCircuitPeeringPropertiesFormat[name]))}
Swagger.field_name(::Type{ ExpressRouteCircuitPeeringPropertiesFormat }, property_name::Symbol) =  _property_map_ExpressRouteCircuitPeeringPropertiesFormat[property_name]

const _allowed_ExpressRouteCircuitPeeringPropertiesFormat_peeringType = ["AzurePublicPeering", "AzurePrivatePeering", "MicrosoftPeering"]

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

function check_required(o::ExpressRouteCircuitPeeringPropertiesFormat)
    true
end

function validate_property(::Type{ ExpressRouteCircuitPeeringPropertiesFormat }, name::Symbol, val)
    if name === Symbol("peeringType")
        Swagger.validate_param(name, "ExpressRouteCircuitPeeringPropertiesFormat", :enum, val, _allowed_ExpressRouteCircuitPeeringPropertiesFormat_peeringType)
    end
    if name === Symbol("state")
        Swagger.validate_param(name, "ExpressRouteCircuitPeeringPropertiesFormat", :enum, val, _allowed_ExpressRouteCircuitPeeringPropertiesFormat_state)
    end
end


text:
using LinearAlgebra, ForwardDiff, Distributions, StaticArrays
using IterativeSolvers
using Plots
using Rotations

# noc
function κ_no(z)
    max.(0.0, z)
end

function Jκ_no(z)
    p = zero(z)
    for (i, pp) in enumerate(z)
        # println(pp)
        # println(i)
        if pp >= 0.0
            p[i] = 1.0
        end
    end
    return Diagonal(p)
end

# soc
function κ_soc(z)
    z1 = z[1:end-1]
    z2 = z[end]

    z_proj = zero(z)

    if norm(z1) <= z2
        z_proj = copy(z)
    elseif norm(z1) <= -z2
        z_proj = zero(z)
    else
        a = 0.5 * (1.0 + z2 / norm(z1))
        z_proj[1:end-1] = a * z1
        z_proj[end] = a * norm(z1)
    end
    return z_proj
end

function Jκ_soc(z)
    z1 = z[1:end-1]
    z2 = z[end]
    m = length(z)

    if norm(z1) <= z2
        return Diagonal(ones(m))
    elseif norm(z1) <= -z2
        return Diagonal(zeros(m))
    else
        D = zeros(m, m)
        for i = 1:m
            if i < m
                D[i, i] = 0.5 + 0.5 * z2 / norm(z1) - 0.5 * z2 * ((z1[i])^2.0) / norm(z1)^3.0
            else
                D[i, i] = 0.5
            end
            for j = 1:m
                if j > i
                    if j < m
                        D[i, j] = -0.5 * z2 * z1[i] * z1[j] / norm(z1)^3.0
                        D[j, i] = -0.5 * z2 * z1[i] * z1[j] / norm(z1)^3.0
                    elseif j == m
                        D[i, j] = 0.5 * z1[i] / norm(z1)
                        D[j, i] = 0.5 * z1[i] / norm(z1)
                    end
                end
            end
        end
        return D
    end
end

# impact problem
struct Box
    m # mass
    J # inertia
    μ # friction coefficient
    g # gravity

    r             # corner length
	n_corners     # number of corners
    corner_offset # precomputed corner offsets

    nq # configurations dimension
end

# Kinematics
r = 0.5
c8 = @SVector [r, r, r]
c2 = @SVector [r, r, -r]
c3 = @SVector [r, -r, r]
c4 = @SVector [r, -r, -r]
c5 = @SVector [-r, r, r]
c6 = @SVector [-r, r, -r]
c7 = @SVector [-r, -r, r]
c1 = @SVector [-r, -r, -r]

cr = @SVector [-r, -r, 0.0]
cl = @SVector [r, r, 0.0]

corner_offset = @SVector [cl, cr]
# corner_offset = @SVector [c1, c2, c3, c4, c5, c6, c7, c8]

model = Box(1.0, 1.0 / 12.0 * 1.0 * ((2.0 * r)^2 + (2.0 * r)^2),
 			1.0, 9.81,
            r, 2, corner_offset,
            6)

# Methods
function mass_matrix(model::Box)
	Diagonal(@SVector [model.m, model.m, model.m,
		model.J, model.J, model.J])
end

function gravity(model::Box)
	@SVector [0., 0., model.m * model.g, 0., 0., 0.]
end

function kinematics(model::Box, q)
    p = view(q, 1:3)
    rot = view(q, 4:6)

    R = MRP(rot...)

    # SVector{24}([(p + R * model.corner_offset[1])...,
    #           (p + R * model.corner_offset[2])...,
    #           (p + R * model.corner_offset[3])...,
    #           (p + R * model.corner_offset[4])...,
    #           (p + R * model.corner_offset[5])...,
    #           (p + R * model.corner_offset[6])...,
    #           (p + R * model.corner_offset[7])...,
    #           (p + R * model.corner_offset[8])...])
	SVector{6}([(p + R * model.corner_offset[1])...,
		(p + R * model.corner_offset[2])...])

end

function jacobian(model::Box, q)
	k(z) = kinematics(model, z)
	ForwardDiff.jacobian(k, q)
end

function ϕ_func(model::Box, q)
    p = view(q, 1:3)
    rot = view(q, 4:6)

    R = MRP(rot...)

    # @SVector [(p + R * model.corner_offset[1])[3],
    #           (p + R * model.corner_offset[2])[3],
    #           (p + R * model.corner_offset[3])[3],
    #           (p + R * model.corner_offset[4])[3],
    #           (p + R * model.corner_offset[5])[3],
    #           (p + R * model.corner_offset[6])[3],
    #           (p + R * model.corner_offset[7])[3],
    #           (p + R * model.corner_offset[8])[3]]
	SVector{2}([(p + R * model.corner_offset[1])[3],
		(p + R * model.corner_offset[2])[3]])
end

# dynamics
function dynamics(model, q1, q2, q3, λ, h)
      nq = model.nq
      SVector{nq}(mass_matrix(model) * (2.0 * q2 - q1 - q3) / h
            - h * gravity(model)
            + h * jacobian(model, q3)' * λ)
end

mrp = MRP(UnitQuaternion(RotY(0.0) * RotX(0.0)))

q1 = [0.0; 0.0; 0.5; mrp.x; mrp.y; mrp.z]
q2 = [0.0; 0.0; 0.5; mrp.x; mrp.y; mrp.z]
h = 0.1

ϕ_func(model, q2)

function r_action(z)
    q3 = z[1:6]
    ϕ = ϕ_func(model, q3)
	# λ = z[7:8]
	# f = [0.0; 0.0; λ[1]; 0.0; 0.0; λ[2]]
    λ = z[7:8]
	s = z[9:10]
	f = zeros(eltype(z), 6)

	for i = 1:model.n_corners
		f[i * 3] = λ[i]
	end

    [dynamics(model, q1, q2, q3, f, h);
	 s - ϕ
     s - κ_no(s - λ)]
end

function solve()
    z = rand(10)
	z[1:6] = q2

    extra_iters = 0

    for i = 1:10
        _F = r_action(z)
        _J = ForwardDiff.jacobian(r_action, z)
        Δ = gmres(_J, 1.0 * _F, abstol = 1.0e-12, maxiter = i + extra_iters)
        # Δ = (_J' * _J + 1.0e-5 * I) \ (_J' * _F)
        iter = 0
        α = 1.0
        while norm(r_action(z - α * Δ))^2.0 >= (1.0 - 0.001 * α) * norm(_F)^2.0 && α > 1.0e-4
            α = 0.5 * α
            # println("   α = $α")
            iter += 1
            if iter > 100
                @error "line search fail"

                return z
            end
        end

        if α <= 1.0e-4
            extra_iters += 1
        end

        println("iter ($i) - norm: $(norm(r_action(z)))")

        z .-= α * Δ
    end

    return z
end

z_sol = solve()
ϕ_func(model, z_sol[1:6])

# impact and friction
function solve(q1, q2, h; init = 0.0, step = :ls)

    z = init * rand(24)
    z[1:6] = copy(q2)

	function r_if(z)
		q3 = z[1:6]
		ϕ = ϕ_func(model, q3)
		λ = z[7:8]
		s = z[9:10]
		b̄ = z[11:16]
		μ = z[17:18]
		η = z[19:24]

		# velocities
		v = jacobian(model, q3) * (q3 - q2) ./ h

		# contact forces
		f = zeros(eltype(z), 6)

		for i = 1:model.n_corners
			f[i * 3] = λ[i]
			f[(i-1) * 3 .+ (1:2)] = b̄[(i-1) * 3 .+ (1:2)]
		end

		[dynamics(model, q1, q2, q3, f, h); # 1:6
		 s - ϕ;
		 s - κ_no(s - λ);
		 [v[1:2]; -μ[1]] - η[1:3];
		 λ[1] - b̄[3];
		 b̄[1:3] - κ_soc(b̄[1:3] - η[1:3]);
		 [v[4:5]; -μ[2]] - η[4:6];
	 	 λ[2] - b̄[6];
	 	 b̄[4:6] - κ_soc(b̄[4:6] - η[4:6])]
	end

	function R_if(z)
		_R = ForwardDiff.jacobian(r_if, z)

		# fix projection
		I2 = Diagonal(ones(2))
		I3 = Diagonal(ones(3))

		λ = z[7:8]
		s = z[9:10]
		b̄ = z[11:16]
		μ = z[17:18]
		η = z[19:24]

		Jno = Jκ_no(s - λ)
		Jsoc1 = Jκ_soc(b̄[1:3] - η[1:3])
		Jsoc2 = Jκ_soc(b̄[4:6] - η[4:6])

		_R[9:10, 9:10] = I2 - Jno
		_R[9:10, 7:8] = Jno
		_R[15:17, 11:13] = I3 - Jsoc1
		_R[15:17, 19:21] = Jsoc1
		_R[22:24, 14:16] = I3 - Jsoc2
		_R[22:24, 22:24] = Jsoc2

		return _R
	end

    extra_iters = 0

    for i = 1:500
        _F = r_if(z)
		if norm(_F) < 1.0e-8
			println("iter ($i) - norm: $(norm(r_if(z)))")
	        return z, true
		end

        _J = R_if(z)
		if step == :ls
        	Δ = (_J' * _J + 1.0e-6 * I) \ (_J' * _F)
		else
			Δ = gmres(_J, 1.0 * _F, abstol = 1.0e-12, maxiter = i + extra_iters)
		end
        iter = 0
        α = 1.0
        while norm(r_if(z - α * Δ))^2.0 >= (1.0 - 0.001 * α) * norm(_F)^2.0 && α > 1.0e-4
            α = 0.5 * α
            # println("   α = $α")
            iter += 1
            if iter > 100
                @error "line search fail"

                return z
            end
        end

        if α <= 1.0e-4
            extra_iters += 1
        end

        println("iter ($i) - norm: $(norm(r_if(z)))")

        z .-= α * Δ
    end

    if norm(r_if(z)) < 1.0e-8
        status = true
    else
        status = false
    end

    return z, status
end

h = 0.1
mrp = MRP(UnitQuaternion(RotY(pi / 3.0) * RotX(pi / 6.0)))

v1 = [1.0; -1.0; 0.0; 0.0; 0.0; 0.0]
q2 = [0.0; 0.0; 1.0; mrp.x; mrp.y; mrp.z]
q1 = q2 - h * v1

ϕ_func(model, q1)
z_sol, = solve(q1, q2, h, init = 0.0)
# z_sol[23:25]

function simulate(q1, q2, T, h)
    println("simulation")
    q = [q1, q2]
    y = [zeros(8)]
    b = [zeros(16)]
    for t = 1:T
        println("   t = $t")
        z_sol, status = solve(q[end-1], q[end], h, init = 0.001, step = :ls)

        if !status
            @warn "failed step (t = $t)"
			z_sol, status = solve(q[end-1], q[end], h, init = 0.001, step = :gmres)

			if !status
				@error "failed step (t = $t) [again]"
				println(ϕ_func(model, q[end]))
				println(ϕ_func(model, z_sol[1:6]))
            	return q, y, b
			end
        else
            nothing
        end
		push!(q, z_sol[1:6])
		# push!(y, z_sol[4])
		# push!(b, z_sol[5:6])
    end

    return q, y, b
end

q_sol, y_sol, b_sol = simulate(q1, q2, 100, h)
q_sol[end]
q_sol[end-1]
plot(hcat(q_sol...)[3:3, :]', xlabel = "")
plot(h .* hcat(y_sol...)', xlabel = "", linetype = :steppost)
plot(hcat([ϕ_func(model, q) for q in q_sol]...)', xlabel = "")
#
plot(hcat(q_sol...)[1:2, :]', xlabel = "")
plot(h * hcat(b_sol...)', linetype = :steppost)


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct QuantUniform <: AbstractDistDelayed
    low::NestedFloat
    high::NestedFloat
    q::NestedFloat
end

function quantuniform(low::Float64, high::Float64, q::Float64)::Float64
    return round(uniform(low, high) / q) * q
end


text:
"""
	$(SIGNATURES)

Bar graph. Simple wrapper for `bar` that applies default settings.
`dataV` can be a matrix. But then series are plotted on top of each other. This usually does not make sense.
"""
# There is no way to not show an xlabel (setting it to `nothing` shows "nothing")
function bar_graph(groupLabelV, dataV; kwargs...)
    @assert size(dataV, 1) == size(groupLabelV, 1);
    args = merge(bar_defaults(), kwargs);
    p = bar(groupLabelV, dataV;  args...);
    return p
end

function test_bar_graph()
    dataM = (1:10) .+ range(1.0, 1.5, length = 4)';
    p = bar_graph(string.(1 : 10), dataM,
        labels = ["Line 1"], xlabel = "x label",
        ylabel = "y label", margin = plot_setting(:subPlotMargin),
        linealpha = plot_setting(:lineAlpha));
    return p
end

"""
	$(SIGNATURES)

Grouped bar graph. Rows of `dataM` are groups.
"""
function grouped_bar_graph(xStrV, dataM; kwargs...)
    @assert size(xStrV, 1) == size(dataM, 1)  "Size mismatch: $(size(xStrV)), $(size(dataM))";
    args = merge(grouped_bar_defaults(), kwargs);
    p = groupedbar(xStrV, dataM; args...);
    return p
end

# ---------------


text:
using Test
using HTTP
using JSON2
using Carrot
using Mocking

Mocking.activate()

response = "{\"instrument\":\"EUR_USD\",\"granularity\":\"H1\",\"candles\":[{\"complete\":true,\"volume\":1013,\"time\":\"2020-06-26T04:00:00.000000000Z\",\"mid\":{\"o\":\"1.12227\",\"h\":\"1.12264\",\"l\":\"1.12211\",\"c\":\"1.12238\"}}]}"
data = JSON2.read(response)

@testset "fetch_data" begin
    @test create_time_series_metrix(data.candles) ==
          [1.12227 1.12264 1.12211 1.12238 1013.0]
    patch = @patch HTTP.request(a::String, b::String, c::Dict{String,String}; query="granularity=H1&count=1500", readtimeout=2) = HTTP.Response(200, body=Vector{UInt8}(response))
    apply(patch) do
        @test getCandles("EUR_GBP", "H1") == [1.12227 1.12264 1.12211 1.12238 1013.0]
    end

end


text:
# Set random seed for the random generator
#  so that our results are reproducible

function gibs(nrep=105000;
              seed=123,
              nb  =5000,
              ρ   =.5)
   # Seed - sets the random number generator seed
   #   so that it is easy to create reproducible results.
   # nrep - in the number of repetitions of the MCMC draws.
   # nb   - number of observations for the "Burn-in"
   # ρ    - Correlation between Y1 and Y2
   srand(seed)

   # Create empty vectors of length nrep
   yy₁ = yy₂ = spzeros(nrep,1)

   # Calculate standard deviation of Y2
   σ   = sqrt(1-ρ^2)

   # Initialize y1 with a single normal draw
   y₁ = randn()*σ

   # Now we start the Gibbs sampler loop
   for i in 1:nrep
     y₂ = randn()*σ + ρ*y₁
     y₁ = randn()*σ + ρ*y₂
     yy₁[i]=y₁
     yy₂[i]=y₂
   end

   return (yy₁[(nb+1):end,1],yy₂[(nb+1):end,1])
end
gibs()

# Set random seed for the random generator
#  so that our results are reproducible
function gibs2(nrep=105000;
               seed=123,
               nb  =5000,
               ρ   =.5)
   # Seed - sets the random number generator seed
   #   so that it is easy to create reproducible results.
   # nrep - in the number of repetitions of the MCMC draws.
   # nb   - number of observations for the "Burn-in"
   # ρ    - Correlation between Y1 and Y2
   srand(seed)

   # Calculate standard deviation of Y2
   sigma = sqrt(1-ρ^2)

   # Create empty vectors of length nrep
   y1 = randn(nrep)
   y2 = randn(nrep)

   # Now we start the Gibbs sampler loop
   for i in 2:nrep
     y2[i] += y1[i-1]*ρ
     y1[i] += y2[i]*ρ
   end

   return (y1[(nb+1):end,1],y2[(nb+1):end,1])
end
gibs2()

@elapsed yy = gibs()
# 0.063151467
@elapsed gibs(10^6)
# 0.479542057

@elapsed yy = gibs2()
# 0.010729382
@elapsed yy = gibs2(10^6)
# 0.065821774

max_stars_repo_name: EconometricsBySimulation/Ninja

max_issues_repo_name: EconometricsBySimulation/Ninja

max_forks_repo_name: EconometricsBySimulation/Ninja


text:
#=

  circuit_path in Julia ConstraintSolver.jl

  circuit_path is the same as circuit(model, x)
  where x is the cicuit. Here we also expose the
  z array which containts the path.

  The advantage of this is that one can impose
  constraints on the path as well as on the circuit.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# The circuit_path() constraint is in constraints_utils.jl
#

#
# The number of solutions for n is (n-1)!
#
# The following is almost identical as circuit_test in circuit.jl
#
function circuit_path_test(n,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                        ))

    @variable(model, 1 <= x[1:n] <= n, Int)
    @variable(model, 1 <= path[1:n] <= n, Int)

    circuit_path(model, x, path)

    # Example of constraint on path
    # @constraint(model, path[1] == 2)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        # println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                pathx = convert.(Integer,JuMP.value.(path; result=sol))
                println("x   :$xx")
                println("path:$pathx")
                println()

            end
        end
        println("num_sols:$num_sols\n")
    else
        println("status:$status")
    end

    return status
end

function circuit_path_test_range(from, to)
    print_solutions=true
    all_solutions=true
    for n in from:to

        println("n:$n")
        if n > 5
            print_solutions = false
        end
        if n > 8
            all_solutions = false
            print_solutions = true
        end
        @time status = circuit_path_test(n,print_solutions,all_solutions)
        if status != MOI.OPTIMAL
            break
        end
    end
end

circuit_path_test_range(2,39)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
export conjugates_init, isconstant, issquarefree, conjugates, angle, cos,
       sin, abs, abs2, sqrt

function isconstant(f::PolyElem)
  return f.length<2
end

function issquarefree(f::PolyElem)
  return isconstant(gcd(f, derivative(f)))
end

function conjugates_init(f::Union{fmpz_poly, fmpq_poly})
  if typeof(f) == fmpq_poly
    f = f*denominator(f)
    g = Array{fmpz}(undef, length(f))
    for i = 1:f.length
      g[i] = FlintZZ(numerator(coeff(f, i-1)))
    end
    g = PolynomialRing(FlintZZ, string(var(parent(f))), cached = false)[1](g)
    f = g
  end
  isconstant(gcd(f, derivative(f))) || error("poly should be square-free")
  c = roots_ctx()
  c.f = f
  r = _roots(f, 100)

  r_d = Array{BigComplex}(undef, 0)
  c_d = Array{BigComplex}(undef, 0)
  for i = 1:length(r)
    rr = BigComplex(r[i])
    if Base.abs2(imag(rr)) < 1e-20
      push!(r_d, rr)
      continue
    end
    if imag(rr) > 0
      push!(c_d, rr)
      continue
    end
  end
  @assert (length(r_d) + 2*length(c_d) == length(r))
  c.r1 = length(r_d)
  c.r2 = length(c_d)
  sort!(r_d, lt = function(a,b) return real(a) < real(b); end)
  sort!(c_d, lt = function(a,b) return angle(a) < angle(b); end)
  c.r_d = vcat(r_d, c_d)

  c.r = Array{BigComplex}(undef, 0)

  old = precision(BigFloat)
  setprecision(53)
  for i = 1:length(c.r_d)
    push!(c.r, c.r_d[i])
  end
  setprecision(old)

  return c
end

function evaluate(f::fmpq_poly, r::BigComplex)
  #Horner - not elegant, but workable
  l = f.length-1
  s = BigComplex(BigFloat(coeff(f, l)))
  for i =l-1:-1:0
    s = s*r+BigComplex(BigFloat(coeff(f, i)))
  end
  return s
end

function evaluate(f::fmpq_poly, r::T) where T <: RingElem
  R = parent(r)
  if iszero(f)
    return zero(R)
  end
  l = length(f) - 1
  s = R(coeff(f, l))
  for i in l-1:-1:0
    s = s*r + R(coeff(f, i))
  end
  return s
end

function evaluate(f::fmpz_poly, r::BigComplex)
  #Horner - not elegant, but workable
  l = f.length-1
  s = BigComplex(coeff(f, l))
  for i =l-1:-1:0
    s = s*r+BigComplex(coeff(f, i))
  end
  return s
end

function hensel_lift(f::fmpz_poly, r::BigComplex)
  return r - evaluate(f, r)/evaluate(derivative(f), r)
end

function conjugates(K::AnticNumberField, p::Int)
  return conjugates(roots_ctx(K), p)
end

function conjugates(c::roots_ctx, p::Int)
  prec = precision(c.r[1])
  while prec < 2*p
    prec *= 2
    for i = 1:length(c.r)
      function _h()
        return hensel_lift(c.f, c.r[i])
      end
      c.r[i] = setprecision(_h, prec)
    end
  end
  return setprecision(c.r, p)
end

function Base.setprecision(a::BigComplex, p::Int)
  return BigComplex(setprecision(a.re, p), setprecision(a.im, p))
end

function Base.setprecision(a::Array{BigComplex, 1}, p::Int)
  b = Array{BigComplex}(undef, 0);
  for i = 1:length(a)
    push!(b, setprecision(a[i], p))
  end
  return b
end

function minkowski(a::nf_elem, p::Int)
  c = roots_ctx(parent(a))
  x = conjugates_arb(a, p)
  old = precision(BigFloat)
  setprecision(BigFloat, p)
  m = Array{BigFloat}(undef, 0)
  for i=1:c.r1
    v = BigFloat(real(x[i]))
    push!(m, v)
  end
  s2 = sqrt(BigFloat(2))
  for i=1:c.r2
    z = x[c.r1+i]
    push!(m, s2*BigFloat(real(z)))
    push!(m, s2*BigFloat(imag(z)))
  end
  setprecision(BigFloat, old)
  return m
end

function length(a::nf_elem, p::Int = 50)
  m = minkowski(a, p)
  return sum([x*x for x in m])
end

function setprecision!(x::BigFloat, p::Int)
  ccall((:mpfr_prec_round, :libmpfr), Nothing, (Ref{BigFloat}, Clong, Int32), x, p, Base.MPFR.ROUNDING_MODE[])
end

function Base.setprecision(x::BigFloat, p::Int)
  setprecision(BigFloat, p) do
    y = BigFloat()
    ccall((:mpfr_set, :libmpfr), Nothing, (Ref{BigFloat}, Ref{BigFloat}, Int32), y, x, Base.MPFR.ROUNDING_MODE[])
    return y
  end
end

function minkowski_matrix(K::AnticNumberField, p::Int = 50)
  c = roots_ctx(K)

  if isdefined(c, :minkowski_matrix)
    if c.minkowski_mat_p == p
      return c.minkowski_matrix
    elseif c.minkowski_mat_p >= p
      return map(x->setprecision(x, p), c.minkowski_matrix)
    end
  end

  old = precision(BigFloat)
  setprecision(p)

  g = gen(K)
  n = degree(K)
  mm = vcat([minkowski(g^i, p) for i=0:n-1])
  m = Array{BigFloat}(undef, n, n)
  for i=1:n
    for j=1:n
      m[i,j] = mm[i][j]
    end
  end
  c.minkowski_matrix = m
  c.minkowski_mat_p = p
  setprecision(old)
  return m
end

function *(a::fmpz_mat, b::Array{BigFloat, 2})
  s = Base.size(b)
  ncols(a) == s[1] || error("dimensions do not match")

  c = Array{BigFloat}(undef, nrows(a), s[2])
  return mult!(c, a, b)
end

for (s,f) in ((:trunc, Base.trunc), (:round, Base.round), (:ceil, Base.ceil), (:floor, Base.floor))
  @eval begin
    function ($s)(a::Array{BigFloat, 2})
      s = Base.size(a)
      m = zero_matrix(FlintZZ, s[1], s[2])
      for i = 1:s[1]
        for j = 1:s[2]
          m[i,j] = FlintZZ(BigInt(($f)(a[i,j])))
        end
      end
      return m
    end
  end
end


text:
# Julia wrapper for header: cheader/tox.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function tox_version_major()
    ccall((:tox_version_major,libtoxcore),Uint32,())
end

function tox_version_minor()
    ccall((:tox_version_minor,libtoxcore),Uint32,())
end

function tox_version_patch()
    ccall((:tox_version_patch,libtoxcore),Uint32,())
end

function tox_version_is_compatible(major::Uint32,minor::Uint32,patch::Uint32)
    ccall((:tox_version_is_compatible,libtoxcore),Bool,(Uint32,Uint32,Uint32),major,minor,patch)
end

function tox_options_default(options::Ptr{Tox_Options})
    ccall((:tox_options_default,libtoxcore),Void,(Ptr{Tox_Options},),options)
end

function tox_options_new(error::Ptr{TOX_ERR_OPTIONS_NEW})
    ccall((:tox_options_new,libtoxcore),Ptr{Tox_Options},(Ptr{TOX_ERR_OPTIONS_NEW},),error)
end

function tox_options_free(options::Ptr{Tox_Options})
    ccall((:tox_options_free,libtoxcore),Void,(Ptr{Tox_Options},),options)
end

function tox_new(options::Ptr{Tox_Options},error::Ptr{TOX_ERR_NEW})
    ccall((:tox_new,libtoxcore),Ptr{Tox},(Ptr{Tox_Options},Ptr{TOX_ERR_NEW}),options,error)
end

function tox_kill(tox::Ptr{Tox})
    ccall((:tox_kill,libtoxcore),Void,(Ptr{Tox},),tox)
end

function tox_get_savedata_size(tox::Ptr{Tox})
    ccall((:tox_get_savedata_size,libtoxcore),Csize_t,(Ptr{Tox},),tox)
end

function tox_get_savedata(tox::Ptr{Tox},savedata::Ptr{Uint8})
    ccall((:tox_get_savedata,libtoxcore),Void,(Ptr{Tox},Ptr{Uint8}),tox,savedata)
end

function tox_bootstrap(tox::Ptr{Tox},address::Ptr{Uint8},port::Uint16,public_key::Ptr{Uint8},error::Ptr{TOX_ERR_BOOTSTRAP})
    ccall((:tox_bootstrap,libtoxcore),Bool,(Ptr{Tox},Ptr{Uint8},Uint16,Ptr{Uint8},Ptr{TOX_ERR_BOOTSTRAP}),tox,address,port,public_key,error)
end

function tox_add_tcp_relay(tox::Ptr{Tox},address::Ptr{Uint8},port::Uint16,public_key::Ptr{Uint8},error::Ptr{TOX_ERR_BOOTSTRAP})
    ccall((:tox_add_tcp_relay,libtoxcore),Bool,(Ptr{Tox},Ptr{Uint8},Uint16,Ptr{Uint8},Ptr{TOX_ERR_BOOTSTRAP}),tox,address,port,public_key,error)
end

function tox_self_get_connection_status(tox::Ptr{Tox})
    ccall((:tox_self_get_connection_status,libtoxcore),TOX_CONNECTION,(Ptr{Tox},),tox)
end

function tox_callback_self_connection_status(tox::Ptr{Tox},callback::Ptr{tox_self_connection_status_cb},user_data::Ptr{Void})
    ccall((:tox_callback_self_connection_status,libtoxcore),Void,(Ptr{Tox},Ptr{tox_self_connection_status_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_iteration_interval(tox::Ptr{Tox})
    ccall((:tox_iteration_interval,libtoxcore),Uint32,(Ptr{Tox},),tox)
end

function tox_iterate(tox::Ptr{Tox})
    ccall((:tox_iterate,libtoxcore),Void,(Ptr{Tox},),tox)
end

function tox_self_get_address(tox::Ptr{Tox},address::Ptr{Uint8})
    ccall((:tox_self_get_address,libtoxcore),Void,(Ptr{Tox},Ptr{Uint8}),tox,address)
end

function tox_self_set_nospam(tox::Ptr{Tox},nospam::Uint32)
    ccall((:tox_self_set_nospam,libtoxcore),Void,(Ptr{Tox},Uint32),tox,nospam)
end

function tox_self_get_nospam(tox::Ptr{Tox})
    ccall((:tox_self_get_nospam,libtoxcore),Uint32,(Ptr{Tox},),tox)
end

function tox_self_get_public_key(tox::Ptr{Tox},public_key::Ptr{Uint8})
    ccall((:tox_self_get_public_key,libtoxcore),Void,(Ptr{Tox},Ptr{Uint8}),tox,public_key)
end

function tox_self_get_secret_key(tox::Ptr{Tox},secret_key::Ptr{Uint8})
    ccall((:tox_self_get_secret_key,libtoxcore),Void,(Ptr{Tox},Ptr{Uint8}),tox,secret_key)
end

function tox_self_set_name(tox::Ptr{Tox},name::Ptr{Uint8},length::Csize_t,error::Ptr{TOX_ERR_SET_INFO})
    ccall((:tox_self_set_name,libtoxcore),Bool,(Ptr{Tox},Ptr{Uint8},Csize_t,Ptr{TOX_ERR_SET_INFO}),tox,name,length,error)
end

function tox_self_get_name_size(tox::Ptr{Tox})
    ccall((:tox_self_get_name_size,libtoxcore),Csize_t,(Ptr{Tox},),tox)
end

function tox_self_get_name(tox::Ptr{Tox},name::Ptr{Uint8})
    ccall((:tox_self_get_name,libtoxcore),Void,(Ptr{Tox},Ptr{Uint8}),tox,name)
end

function tox_self_set_status_message(tox::Ptr{Tox},status_message::Ptr{Uint8},length::Csize_t,error::Ptr{TOX_ERR_SET_INFO})
    ccall((:tox_self_set_status_message,libtoxcore),Bool,(Ptr{Tox},Ptr{Uint8},Csize_t,Ptr{TOX_ERR_SET_INFO}),tox,status_message,length,error)
end

function tox_self_get_status_message_size(tox::Ptr{Tox})
    ccall((:tox_self_get_status_message_size,libtoxcore),Csize_t,(Ptr{Tox},),tox)
end

function tox_self_get_status_message(tox::Ptr{Tox},status_message::Ptr{Uint8})
    ccall((:tox_self_get_status_message,libtoxcore),Void,(Ptr{Tox},Ptr{Uint8}),tox,status_message)
end

function tox_self_set_status(tox::Ptr{Tox},status::TOX_USER_STATUS)
    ccall((:tox_self_set_status,libtoxcore),Void,(Ptr{Tox},TOX_USER_STATUS),tox,status)
end

function tox_self_get_status(tox::Ptr{Tox})
    ccall((:tox_self_get_status,libtoxcore),TOX_USER_STATUS,(Ptr{Tox},),tox)
end

function tox_friend_add(tox::Ptr{Tox},address::Ptr{Uint8},message::Ptr{Uint8},length::Csize_t,error::Ptr{TOX_ERR_FRIEND_ADD})
    ccall((:tox_friend_add,libtoxcore),Uint32,(Ptr{Tox},Ptr{Uint8},Ptr{Uint8},Csize_t,Ptr{TOX_ERR_FRIEND_ADD}),tox,address,message,length,error)
end

function tox_friend_add_norequest(tox::Ptr{Tox},public_key::Ptr{Uint8},error::Ptr{TOX_ERR_FRIEND_ADD})
    ccall((:tox_friend_add_norequest,libtoxcore),Uint32,(Ptr{Tox},Ptr{Uint8},Ptr{TOX_ERR_FRIEND_ADD}),tox,public_key,error)
end

function tox_friend_delete(tox::Ptr{Tox},friend_number::Uint32,error::Ptr{TOX_ERR_FRIEND_DELETE})
    ccall((:tox_friend_delete,libtoxcore),Bool,(Ptr{Tox},Uint32,Ptr{TOX_ERR_FRIEND_DELETE}),tox,friend_number,error)
end

function tox_friend_by_public_key(tox::Ptr{Tox},public_key::Ptr{Uint8},error::Ptr{TOX_ERR_FRIEND_BY_PUBLIC_KEY})
    ccall((:tox_friend_by_public_key,libtoxcore),Uint32,(Ptr{Tox},Ptr{Uint8},Ptr{TOX_ERR_FRIEND_BY_PUBLIC_KEY}),tox,public_key,error)
end

function tox_friend_exists(tox::Ptr{Tox},friend_number::Uint32)
    ccall((:tox_friend_exists,libtoxcore),Bool,(Ptr{Tox},Uint32),tox,friend_number)
end

function tox_self_get_friend_list_size(tox::Ptr{Tox})
    ccall((:tox_self_get_friend_list_size,libtoxcore),Csize_t,(Ptr{Tox},),tox)
end

function tox_self_get_friend_list(tox::Ptr{Tox},friend_list::Ptr{Uint32})
    ccall((:tox_self_get_friend_list,libtoxcore),Void,(Ptr{Tox},Ptr{Uint32}),tox,friend_list)
end

function tox_friend_get_public_key(tox::Ptr{Tox},friend_number::Uint32,public_key::Ptr{Uint8},error::Ptr{TOX_ERR_FRIEND_GET_PUBLIC_KEY})
    ccall((:tox_friend_get_public_key,libtoxcore),Bool,(Ptr{Tox},Uint32,Ptr{Uint8},Ptr{TOX_ERR_FRIEND_GET_PUBLIC_KEY}),tox,friend_number,public_key,error)
end

function tox_friend_get_last_online(tox::Ptr{Tox},friend_number::Uint32,error::Ptr{TOX_ERR_FRIEND_GET_LAST_ONLINE})
    ccall((:tox_friend_get_last_online,libtoxcore),Uint64,(Ptr{Tox},Uint32,Ptr{TOX_ERR_FRIEND_GET_LAST_ONLINE}),tox,friend_number,error)
end

function tox_friend_get_name_size(tox::Ptr{Tox},friend_number::Uint32,error::Ptr{TOX_ERR_FRIEND_QUERY})
    ccall((:tox_friend_get_name_size,libtoxcore),Csize_t,(Ptr{Tox},Uint32,Ptr{TOX_ERR_FRIEND_QUERY}),tox,friend_number,error)
end

function tox_friend_get_name(tox::Ptr{Tox},friend_number::Uint32,name::Ptr{Uint8},error::Ptr{TOX_ERR_FRIEND_QUERY})
    ccall((:tox_friend_get_name,libtoxcore),Bool,(Ptr{Tox},Uint32,Ptr{Uint8},Ptr{TOX_ERR_FRIEND_QUERY}),tox,friend_number,name,error)
end

function tox_callback_friend_name(tox::Ptr{Tox},callback::Ptr{tox_friend_name_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_name,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_name_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_friend_get_status_message_size(tox::Ptr{Tox},friend_number::Uint32,error::Ptr{TOX_ERR_FRIEND_QUERY})
    ccall((:tox_friend_get_status_message_size,libtoxcore),Csize_t,(Ptr{Tox},Uint32,Ptr{TOX_ERR_FRIEND_QUERY}),tox,friend_number,error)
end

function tox_friend_get_status_message(tox::Ptr{Tox},friend_number::Uint32,status_message::Ptr{Uint8},error::Ptr{TOX_ERR_FRIEND_QUERY})
    ccall((:tox_friend_get_status_message,libtoxcore),Bool,(Ptr{Tox},Uint32,Ptr{Uint8},Ptr{TOX_ERR_FRIEND_QUERY}),tox,friend_number,status_message,error)
end

function tox_callback_friend_status_message(tox::Ptr{Tox},callback::Ptr{tox_friend_status_message_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_status_message,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_status_message_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_friend_get_status(tox::Ptr{Tox},friend_number::Uint32,error::Ptr{TOX_ERR_FRIEND_QUERY})
    ccall((:tox_friend_get_status,libtoxcore),TOX_USER_STATUS,(Ptr{Tox},Uint32,Ptr{TOX_ERR_FRIEND_QUERY}),tox,friend_number,error)
end

function tox_callback_friend_status(tox::Ptr{Tox},callback::Ptr{tox_friend_status_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_status,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_status_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_friend_get_connection_status(tox::Ptr{Tox},friend_number::Uint32,error::Ptr{TOX_ERR_FRIEND_QUERY})
    ccall((:tox_friend_get_connection_status,libtoxcore),TOX_CONNECTION,(Ptr{Tox},Uint32,Ptr{TOX_ERR_FRIEND_QUERY}),tox,friend_number,error)
end

function tox_callback_friend_connection_status(tox::Ptr{Tox},callback::Ptr{tox_friend_connection_status_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_connection_status,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_connection_status_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_friend_get_typing(tox::Ptr{Tox},friend_number::Uint32,error::Ptr{TOX_ERR_FRIEND_QUERY})
    ccall((:tox_friend_get_typing,libtoxcore),Bool,(Ptr{Tox},Uint32,Ptr{TOX_ERR_FRIEND_QUERY}),tox,friend_number,error)
end

function tox_callback_friend_typing(tox::Ptr{Tox},callback::Ptr{tox_friend_typing_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_typing,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_typing_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_self_set_typing(tox::Ptr{Tox},friend_number::Uint32,typing::Bool,error::Ptr{TOX_ERR_SET_TYPING})
    ccall((:tox_self_set_typing,libtoxcore),Bool,(Ptr{Tox},Uint32,Bool,Ptr{TOX_ERR_SET_TYPING}),tox,friend_number,typing,error)
end

function tox_friend_send_message(tox::Ptr{Tox},friend_number::Uint32,_type::TOX_MESSAGE_TYPE,message::Ptr{Uint8},length::Csize_t,error::Ptr{TOX_ERR_FRIEND_SEND_MESSAGE})
    ccall((:tox_friend_send_message,libtoxcore),Uint32,(Ptr{Tox},Uint32,TOX_MESSAGE_TYPE,Ptr{Uint8},Csize_t,Ptr{TOX_ERR_FRIEND_SEND_MESSAGE}),tox,friend_number,_type,message,length,error)
end

function tox_callback_friend_read_receipt(tox::Ptr{Tox},callback::Ptr{tox_friend_read_receipt_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_read_receipt,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_read_receipt_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_callback_friend_request(tox::Ptr{Tox},callback::Ptr{tox_friend_request_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_request,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_request_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_callback_friend_message(tox::Ptr{Tox},callback::Ptr{tox_friend_message_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_message,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_message_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_hash(hash::Ptr{Uint8},data::Ptr{Uint8},length::Csize_t)
    ccall((:tox_hash,libtoxcore),Bool,(Ptr{Uint8},Ptr{Uint8},Csize_t),hash,data,length)
end

function tox_file_control(tox::Ptr{Tox},friend_number::Uint32,file_number::Uint32,control::TOX_FILE_CONTROL,error::Ptr{TOX_ERR_FILE_CONTROL})
    ccall((:tox_file_control,libtoxcore),Bool,(Ptr{Tox},Uint32,Uint32,TOX_FILE_CONTROL,Ptr{TOX_ERR_FILE_CONTROL}),tox,friend_number,file_number,control,error)
end

function tox_callback_file_recv_control(tox::Ptr{Tox},callback::Ptr{tox_file_recv_control_cb},user_data::Ptr{Void})
    ccall((:tox_callback_file_recv_control,libtoxcore),Void,(Ptr{Tox},Ptr{tox_file_recv_control_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_file_seek(tox::Ptr{Tox},friend_number::Uint32,file_number::Uint32,position::Uint64,error::Ptr{TOX_ERR_FILE_SEEK})
    ccall((:tox_file_seek,libtoxcore),Bool,(Ptr{Tox},Uint32,Uint32,Uint64,Ptr{TOX_ERR_FILE_SEEK}),tox,friend_number,file_number,position,error)
end

function tox_file_get_file_id(tox::Ptr{Tox},friend_number::Uint32,file_number::Uint32,file_id::Ptr{Uint8},error::Ptr{TOX_ERR_FILE_GET})
    ccall((:tox_file_get_file_id,libtoxcore),Bool,(Ptr{Tox},Uint32,Uint32,Ptr{Uint8},Ptr{TOX_ERR_FILE_GET}),tox,friend_number,file_number,file_id,error)
end

function tox_file_send(tox::Ptr{Tox},friend_number::Uint32,kind::Uint32,file_size::Uint64,file_id::Ptr{Uint8},filename::Ptr{Uint8},filename_length::Csize_t,error::Ptr{TOX_ERR_FILE_SEND})
    ccall((:tox_file_send,libtoxcore),Uint32,(Ptr{Tox},Uint32,Uint32,Uint64,Ptr{Uint8},Ptr{Uint8},Csize_t,Ptr{TOX_ERR_FILE_SEND}),tox,friend_number,kind,file_size,file_id,filename,filename_length,error)
end

function tox_file_send_chunk(tox::Ptr{Tox},friend_number::Uint32,file_number::Uint32,position::Uint64,data::Ptr{Uint8},length::Csize_t,error::Ptr{TOX_ERR_FILE_SEND_CHUNK})
    ccall((:tox_file_send_chunk,libtoxcore),Bool,(Ptr{Tox},Uint32,Uint32,Uint64,Ptr{Uint8},Csize_t,Ptr{TOX_ERR_FILE_SEND_CHUNK}),tox,friend_number,file_number,position,data,length,error)
end

function tox_callback_file_chunk_request(tox::Ptr{Tox},callback::Ptr{tox_file_chunk_request_cb},user_data::Ptr{Void})
    ccall((:tox_callback_file_chunk_request,libtoxcore),Void,(Ptr{Tox},Ptr{tox_file_chunk_request_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_callback_file_recv(tox::Ptr{Tox},callback::Ptr{tox_file_recv_cb},user_data::Ptr{Void})
    ccall((:tox_callback_file_recv,libtoxcore),Void,(Ptr{Tox},Ptr{tox_file_recv_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_callback_file_recv_chunk(tox::Ptr{Tox},callback::Ptr{tox_file_recv_chunk_cb},user_data::Ptr{Void})
    ccall((:tox_callback_file_recv_chunk,libtoxcore),Void,(Ptr{Tox},Ptr{tox_file_recv_chunk_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_friend_send_lossy_packet(tox::Ptr{Tox},friend_number::Uint32,data::Ptr{Uint8},length::Csize_t,error::Ptr{TOX_ERR_FRIEND_CUSTOM_PACKET})
    ccall((:tox_friend_send_lossy_packet,libtoxcore),Bool,(Ptr{Tox},Uint32,Ptr{Uint8},Csize_t,Ptr{TOX_ERR_FRIEND_CUSTOM_PACKET}),tox,friend_number,data,length,error)
end

function tox_friend_send_lossless_packet(tox::Ptr{Tox},friend_number::Uint32,data::Ptr{Uint8},length::Csize_t,error::Ptr{TOX_ERR_FRIEND_CUSTOM_PACKET})
    ccall((:tox_friend_send_lossless_packet,libtoxcore),Bool,(Ptr{Tox},Uint32,Ptr{Uint8},Csize_t,Ptr{TOX_ERR_FRIEND_CUSTOM_PACKET}),tox,friend_number,data,length,error)
end

function tox_callback_friend_lossy_packet(tox::Ptr{Tox},callback::Ptr{tox_friend_lossy_packet_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_lossy_packet,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_lossy_packet_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_callback_friend_lossless_packet(tox::Ptr{Tox},callback::Ptr{tox_friend_lossless_packet_cb},user_data::Ptr{Void})
    ccall((:tox_callback_friend_lossless_packet,libtoxcore),Void,(Ptr{Tox},Ptr{tox_friend_lossless_packet_cb},Ptr{Void}),tox,callback,user_data)
end

function tox_self_get_dht_id(tox::Ptr{Tox},dht_id::Ptr{Uint8})
    ccall((:tox_self_get_dht_id,libtoxcore),Void,(Ptr{Tox},Ptr{Uint8}),tox,dht_id)
end

function tox_self_get_udp_port(tox::Ptr{Tox},error::Ptr{TOX_ERR_GET_PORT})
    ccall((:tox_self_get_udp_port,libtoxcore),Uint16,(Ptr{Tox},Ptr{TOX_ERR_GET_PORT}),tox,error)
end

function tox_self_get_tcp_port(tox::Ptr{Tox},error::Ptr{TOX_ERR_GET_PORT})
    ccall((:tox_self_get_tcp_port,libtoxcore),Uint16,(Ptr{Tox},Ptr{TOX_ERR_GET_PORT}),tox,error)
end

function tox_callback_group_invite(tox::Ptr{Tox},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_invite,libtoxcore),Void,(Ptr{Tox},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_callback_group_message(tox::Ptr{Tox},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_message,libtoxcore),Void,(Ptr{Tox},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_callback_group_action(tox::Ptr{Tox},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_action,libtoxcore),Void,(Ptr{Tox},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_callback_group_title(tox::Ptr{Tox},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_title,libtoxcore),Void,(Ptr{Tox},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_callback_group_namelist_change(tox::Ptr{Tox},_function::Ptr{Void},userdata::Ptr{Void})
    ccall((:tox_callback_group_namelist_change,libtoxcore),Void,(Ptr{Tox},Ptr{Void},Ptr{Void}),tox,_function,userdata)
end

function tox_add_groupchat(tox::Ptr{Tox})
    ccall((:tox_add_groupchat,libtoxcore),Cint,(Ptr{Tox},),tox)
end

function tox_del_groupchat(tox::Ptr{Tox},groupnumber::Cint)
    ccall((:tox_del_groupchat,libtoxcore),Cint,(Ptr{Tox},Cint),tox,groupnumber)
end

function tox_group_peername(tox::Ptr{Tox},groupnumber::Cint,peernumber::Cint,name::Ptr{Uint8})
    ccall((:tox_group_peername,libtoxcore),Cint,(Ptr{Tox},Cint,Cint,Ptr{Uint8}),tox,groupnumber,peernumber,name)
end

function tox_group_peer_pubkey(tox::Ptr{Tox},groupnumber::Cint,peernumber::Cint,public_key::Ptr{Uint8})
    ccall((:tox_group_peer_pubkey,libtoxcore),Cint,(Ptr{Tox},Cint,Cint,Ptr{Uint8}),tox,groupnumber,peernumber,public_key)
end

function tox_invite_friend(tox::Ptr{Tox},friendnumber::Int32,groupnumber::Cint)
    ccall((:tox_invite_friend,libtoxcore),Cint,(Ptr{Tox},Int32,Cint),tox,friendnumber,groupnumber)
end

function tox_join_groupchat(tox::Ptr{Tox},friendnumber::Int32,data::Ptr{Uint8},length::Uint16)
    ccall((:tox_join_groupchat,libtoxcore),Cint,(Ptr{Tox},Int32,Ptr{Uint8},Uint16),tox,friendnumber,data,length)
end

function tox_group_message_send(tox::Ptr{Tox},groupnumber::Cint,message::Ptr{Uint8},length::Uint16)
    ccall((:tox_group_message_send,libtoxcore),Cint,(Ptr{Tox},Cint,Ptr{Uint8},Uint16),tox,groupnumber,message,length)
end

function tox_group_action_send(tox::Ptr{Tox},groupnumber::Cint,action::Ptr{Uint8},length::Uint16)
    ccall((:tox_group_action_send,libtoxcore),Cint,(Ptr{Tox},Cint,Ptr{Uint8},Uint16),tox,groupnumber,action,length)
end

function tox_group_set_title(tox::Ptr{Tox},groupnumber::Cint,title::Ptr{Uint8},length::Uint8)
    ccall((:tox_group_set_title,libtoxcore),Cint,(Ptr{Tox},Cint,Ptr{Uint8},Uint8),tox,groupnumber,title,length)
end

function tox_group_get_title(tox::Ptr{Tox},groupnumber::Cint,title::Ptr{Uint8},max_length::Uint32)
    ccall((:tox_group_get_title,libtoxcore),Cint,(Ptr{Tox},Cint,Ptr{Uint8},Uint32),tox,groupnumber,title,max_length)
end

function tox_group_peernumber_is_ours(tox::Ptr{Tox},groupnumber::Cint,peernumber::Cint)
    ccall((:tox_group_peernumber_is_ours,libtoxcore),Uint32,(Ptr{Tox},Cint,Cint),tox,groupnumber,peernumber)
end

function tox_group_number_peers(tox::Ptr{Tox},groupnumber::Cint)
    ccall((:tox_group_number_peers,libtoxcore),Cint,(Ptr{Tox},Cint),tox,groupnumber)
end

function tox_group_get_names(tox::Ptr{Tox},groupnumber::Cint,names::Ptr{Array_128_Uint8},lengths::Ptr{Uint16},length::Uint16)
    ccall((:tox_group_get_names,libtoxcore),Cint,(Ptr{Tox},Cint,Ptr{Array_128_Uint8},Ptr{Uint16},Uint16),tox,groupnumber,names,lengths,length)
end

function tox_count_chatlist(tox::Ptr{Tox})
    ccall((:tox_count_chatlist,libtoxcore),Uint32,(Ptr{Tox},),tox)
end

function tox_get_chatlist(tox::Ptr{Tox},out_list::Ptr{Int32},list_size::Uint32)
    ccall((:tox_get_chatlist,libtoxcore),Uint32,(Ptr{Tox},Ptr{Int32},Uint32),tox,out_list,list_size)
end

function tox_group_get_type(tox::Ptr{Tox},groupnumber::Cint)
    ccall((:tox_group_get_type,libtoxcore),Cint,(Ptr{Tox},Cint),tox,groupnumber)
end


text:
# some tests for ALSPGrad

# data

srand(5678)

p = 5
n = 8
k = 3

# Matrices
for T in (Float64, Float32)
    Wg = max.(rand(T, p, k) .- T(0.3), zero(T))
    Hg = max.(rand(T, k, n) .- T(0.3), zero(T))
    X = Wg * Hg

    # test update of H

    H = rand(T, k, n)
    NMF.alspgrad_updateh!(X, Wg, H; maxiter=200)
    @test all(H .>= 0.0)
    @test H ≈ Hg atol=eps(T)^(1/4)

    # test update of W

    W = rand(T, p, k)
    NMF.alspgrad_updatew!(X, W, Hg; maxiter=200)
    @test all(W .>= 0.0)
    @test W ≈ Wg atol=eps(T)^(1/4)

    NMF.solve!(NMF.ALSPGrad{T}(), X, W, H)
end


text:
"""
    MicrophysicsParameters

Module containing 1-moment bulk microphysics parameters.
"""
module MicrophysicsParameters
using CLIMA.ParametersType

@exportparameter MP_n_0 8e6 * 2 "Marshall-Palmer distribution n_0 coeff [1/m4]"
@exportparameter C_drag 0.55 "drag coefficient for rain drops [-]"

@exportparameter τ_cond_evap 10 "condensation/evaporation timescale [s]"
@exportparameter τ_sub_resub 10 "sublimation/resublimation timescale [s]"

@exportparameter q_liq_threshold 5e-4 "autoconversion threshold [-]  ∈(0.5, 1) * 1e-3 "
@exportparameter τ_acnv 1e3 "autoconversion timescale [s]  ∈(1e3, 1e4) "

@exportparameter E_col 0.8 "collision efficiency [-]"

@exportparameter a_vent 1.5 "ventilation factor coefficient [-]"
@exportparameter b_vent 0.53 "ventilation factor coefficient [-]"
@exportparameter K_therm 2.4e-2 "thermal conductivity of air [J/m/s/K] "
@exportparameter D_vapor 2.26e-5 "diffusivity of water vapor [m2/s]"
@exportparameter ν_air 1.6e-5 "kinematic viscosity of air [m2/s]"
@exportparameter N_Sc ν_air / D_vapor "Schmidt number [-]"
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
@doc doc"""
Minimise une fonction en utilisant l'algorithme des régions de confiance avec
    - le pas de Cauchy
ou
    - le pas issu de l'algorithme du gradient conjugue tronqué

# Syntaxe
```julia
xk, nb_iters, f(xk), flag = Regions_De_Confiance(algo,f,gradf,hessf,x0,option)
```

# Entrées :

   * **algo**        : (String) string indicant la méthode à utiliser pour calculer le pas
        - **"gct"**   : pour l'algorithme du gradient conjugué tronqué
        - **"cauchy"**: pour le pas de Cauchy
   * **f**           : (Function) la fonction à minimiser
   * **gradf**       : (Function) le gradient de la fonction f
   * **hessf**       : (Function) la hessiene de la fonction à minimiser
   * **x0**          : (Array{Float,1}) point de départ
   * **options**     : (Array{Float,1})
     * **deltaMax**      : utile pour les m-à-j de la région de confiance
                      ``R_{k}=\left\{x_{k}+s ;\|s\| \leq \Delta_{k}\right\}``
     * **gamma1,gamma2** : ``0 < \gamma_{1} < 1 < \gamma_{2}`` pour les m-à-j de ``R_{k}``
     * **eta1,eta2**     : ``0 < \eta_{1} < \eta_{2} < 1`` pour les m-à-j de ``R_{k}``
     * **delta0**        : le rayon de départ de la région de confiance
     * **max_iter**      : le nombre maximale d'iterations
     * **Tol_abs**       : la tolérence absolue
     * **Tol_rel**       : la tolérence relative

# Sorties:

   * **xmin**    : (Array{Float,1}) une approximation de la solution du problème : ``min_{x \in \mathbb{R}^{n}} f(x)``
   * **fxmin**   : (Float) ``f(x_{min})``
   * **flag**    : (Integer) un entier indiquant le critère sur lequel le programme à arrêter
      - **0**    : Convergence
      - **1**    : stagnation du ``x``
      - **2**    : stagnation du ``f``
      - **3**    : nombre maximal d'itération dépassé
   * **nb_iters** : (Integer)le nombre d'iteration qu'à fait le programme

# Exemple d'appel
```julia
algo="gct"
f(x)=100*(x[2]-x[1]^2)^2+(1-x[1])^2
gradf(x)=[-400*x[1]*(x[2]-x[1]^2)-2*(1-x[1]) ; 200*(x[2]-x[1]^2)]
hessf(x)=[-400*(x[2]-3*x[1]^2)+2  -400*x[1];-400*x[1]  200]
x0 = [1; 0]
options = []
xmin, fxmin, flag,nb_iters = Regions_De_Confiance(algo,f,gradf,hessf,x0,options)
```
"""
function Regions_De_Confiance(algo, f::Function, gradf::Function, hessf::Function, x0, options)

    if options == []
        deltaMax = 10
        gamma1 = 0.5
        gamma2 = 2.00
        eta1 = 0.25
        eta2 = 0.75
        delta0 = 2
        max_iter = 30000
        Tol_abs = sqrt(eps())
        Tol_rel = 1e-15
    else
        deltaMax = options[1]
        gamma1 = options[2]
        gamma2 = options[3]
        eta1 = options[4]
        eta2 = options[5]
        delta0 = options[6]
        max_iter = options[7]
        Tol_abs = options[8]
        Tol_rel = options[9]
    end

    n = length(x0)
    nb_iters = 0
    xk = x0
    xk_1 = x0
    delta = delta0

    for i in 1 : max_iter
        nb_iters = i
        xk = xk_1

        if algo == "gct"
            sk = Gradient_Conjugue_Tronque(gradf(xk), hessf(xk), [delta , 3000, Tol_rel])
        elseif algo == "cauchy"
            sk, ~ = Pas_De_Cauchy(gradf(xk), hessf(xk), delta)
        end

        q(s) = f(xk) + gradf(xk)' * (s-xk) + (1 / 2) * (s - xk)' * hessf(xk) * (s - xk)
        f_xk_sk = f(xk + sk)
        rho_k = (f(xk) - f(xk + sk)) / (q(xk) - q(xk + sk))
        if rho_k >= eta1
            xk_1 = xk + sk
        else
            xk_1 = xk
        end

        if rho_k >= eta2
            delta = min(gamma2 * delta, deltaMax)
        elseif rho_k >= eta1
            delta = delta
        else
            delta = gamma1 * delta
        end

        if norm(gradf(xk_1)) < max(Tol_abs, Tol_rel * norm(gradf(x0)))
            flag = 0
            break
        elseif norm(sk) < max(Tol_abs, Tol_rel * norm(xk))
            flag = 1
            break
        elseif abs(f(xk + sk) - f(xk)) < max(Tol_abs, Tol_rel * abs(f(xk)))
            flag = 2
            break
        end
    end

    if nb_iters == max_iter
        flag = 3
    end

    xmin = xk_1
    f_min = f(xmin)

    return xmin, f_min, flag, nb_iters
end

max_stars_repo_name: MOUDDENEHamza/ENSEEIHT

max_issues_repo_name: MOUDDENEHamza/ENSEEIHT

max_forks_repo_name: MOUDDENEHamza/ENSEEIHT


text:
struct rfd_parametrization
    p_raw::Polyhedron   # original polytope
    p_fd::Polyhedron    # full dimensional polytope
    p_opt::Polyhedron  # rounded, full dimensional polytope
    at_opt2raw::AffineMap
    at_opt2fd::AffineMap
    at_fd2raw::AffineMap

    cheby_p_fd::Vector{Float64}
    cheby_p_opt::Vector{Float64}
end

function eval_opt2raw( rfd_parametrization , x_opt::Array{Float64,2} )
    return rfd_parametrization.at_opt2raw(x_opt)
end

function eval_opt2raw( rfd_parametrization , x_opt::Array{Float64,1} )
    return rfd_parametrization.at_opt2raw(x_opt)
end

function get_rounded_fulldimensional_parametrization( fp::Polyhedron ; remove_red_constr=0)

    ( fp_fd , at_fd ) = MetabolicFluxAnalysis.find_minimal_affine_subspace(fp)
    if(remove_red_constr==1)
        fp_fd_a = MetabolicFluxAnalysis.remove_redundant_constraints( fp_fd )
        fp_fd = fp_fd_a
    end
    ( fp_rounded , at_rounded ) = round_polytope( fp_fd )
    at_opt2fd   = at_rounded
    at_opt2flux = compose(at_fd , at_rounded)

    cheby_p_opt  = zeros(size(fp_rounded.G,2))
    cheby_p_fd   = at_opt2fd(cheby_p_opt)

    return rfd_parametrization( copy(fp) , fp_fd , fp_rounded , at_opt2flux , at_opt2fd , at_fd , cheby_p_fd , cheby_p_opt  )
end


text:
import Printf
import Random
using LinearAlgebra

import DifferentialEquations

import Distributions

import PyPlot

import StatsBase
import StatsFuns

import Printf

import Utilities # https://gitlab.com/RoyCCWang/utilities
import Calculus
import AdaptiveRKHS
import Statistics

include("../src/misc/declarations.jl")

include("../src/RKHS/RKHS.jl")
include("../src/RKHS/kernel.jl")
include("../src/RKHS/interpolators.jl")

include("../src/warp_map/Rieszanalysis.jl")
include("../src/misc/utilities.jl")

include("../src/RKHS/querying.jl")

PyPlot.close("all")
fig_num = 1

Random.seed!(25)

D = 3

σ_oracle = 10.0
ρ_oracle = 28.0
β_oracle = 8/3

params_oracle = [σ_oracle; ρ_oracle; β_oracle]

# p is state vector.
function lorenz!(du::Vector{T}, u, p, t) where T
    σ = p[1]
    ρ = p[2]
    β = p[3]

    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]

    return nothing
end

tspan = (0.0, 100.0)

N_obs = 15 #30 #15
time_stamp_range = LinRange(0.1, tspan[end], N_obs)
time_stamp = collect(time_stamp_range)

u0 = zeros(Float64, D)
u0[1] = 1.0

prob = DifferentialEquations.ODEProblem(lorenz!, u0, tspan, params_oracle, dense = true)
sol = DifferentialEquations.solve(prob,
            DifferentialEquations.Tsit5(),
            reltol = 1e-8, abstol = 1e-8)
y_clean = sol.(time_stamp)

import Plots
Plots.plot(sol, vars = (1,2,3), show = true)

@assert 1==2

N_display = 1000
t_display = LinRange(tspan[1], tspan[2], N_display)

sol_x1_display = collect( sol.u[i][1] for i = 1:length(sol.u) )
sol_x2_display = collect( sol.u[i][2] for i = 1:length(sol.u) )
#y1_display = collect( y[i][d_select] for i = 1:length(y) )

t_dummy = 1.0

y_clean = y_clean #./ 100
σ = 5.0 #/100
y = collect( y_clean[i] + randn(D_state) .* σ for i = 1:length(y_clean) )

y1_clean = collect( y_clean[i][1] for i = 1:length(y_clean) )
y1_itp, dy1_itp, d2y1_itp = Utilities.setupcubicitp(y1_clean, [ time_stamp], 1.0)

dy1_ND = xx->Calculus.gradient(y1_itp, xx)[1]
dy1_itp_no_vec = xx->dy1_itp(xx)[1]
dy2_itp_no_vec = xx->d2y1_itp(xx)[1]

### try RKHS solution since.
θ_oracle = AdaptiveRKHS.GaussianKernel1DType(1.0/100)

σ²_oracle = 1e-7

X_oracle = collect( [sol.t[i]] for i = 1:length(sol.t) )

u_array = Vector{Vector{Float64}}(undef,2)
u_array[1] = collect( sol.u[i][1] for i = 1:length(sol.u) )
u_array[2] = collect( sol.u[i][2] for i = 1:length(sol.u) )

u_GP = fitDEsolutionGP(θ_oracle, σ²_oracle, X_oracle, u_array)
u1_GP = u_GP[1]
u2_GP = u_GP[2]

du_GP = tt->evalpreyderivativeswrttime(u1_GP, u2_GP, p_oracle, tt)

du1_GP = tt->du_GP(tt)[1]
du1_GP_ND = tt->Calculus.gradient(u1_GP, tt)[1]
#d2u1_GP_ND = tt->Calculus.hessian(u1_GP, tt)[1]
d2u1_GP = tt->Calculus.gradient(du1_GP, tt)[1]

du2_GP = tt->du_GP(tt)[1]
du2_GP_ND = tt->Calculus.gradient(u2_GP, tt)[1]

Nq = 500
# xq_range = LinRange(sol.t[1], sol.t[end], Nq)
xq_range = LinRange(time_stamp[1], time_stamp[end], Nq)
xq = collect( [ xq_range[i] ] for i = 1:length(xq_range) )

title_string = Printf.@sprintf("Predator-prey populations")

PyPlot.figure(fig_num)
fig_num += 1
PyPlot.plot(sol.t, sol_x1_display, label = "species 1")
#PyPlot.plot(sol.t, sol_x2_display, label = "species 2")
PyPlot.plot(xq, y1_itp.(xq), label = "y1 itp")
PyPlot.plot(xq, u1_GP.(xq), "x", label = "u1 GP")
PyPlot.title(title_string)
PyPlot.legend()

@assert 1==2

title_string = "derivative"
PyPlot.figure(fig_num)
fig_num += 1
#PyPlot.plot(xq, dy1_ND.(xq), label = "numerical itp")
#PyPlot.plot(xq, dy1_itp_no_vec.(xq), "x", label = "analytical itp")
PyPlot.plot(xq, du1_GP_ND.(xq), "x", label = "numerical GP")
PyPlot.plot(xq, du1_GP.(xq), label = "analytical GP")
PyPlot.title(title_string)
PyPlot.legend()

#@assert 1==2

println("norm(dy1_ND.(xq) - dy1_itp_no_vec.(xq)) = ", norm(dy1_ND.(xq) - dy1_itp_no_vec.(xq)))
println()

# fit RKHS.

θ = AdaptiveRKHS.GaussianKernel1DType(1.0/100)

σ²_RKHS_initial = sqrt(σ) #σ^2 #1e-3

X = collect( [time_stamp_vec[i]] for i = 1:length(time_stamp_vec) )
y1 = collect( y[i][1] for i = 1:length(y) )
η = AdaptiveRKHS.RKHSProblemType( zeros(Float64,length(X)),
                     X,
                     θ,
                     σ²_RKHS_initial)
offset_c = Statistics.mean(y1)
y1_c = y1 .- offset_c
AdaptiveRKHS.fitRKHS!(η, y1_c)

fq = xx->sum( η.c[i]*AdaptiveRKHS.evalkernel(xx, η.X[i], η.θ) for i = 1:length(η.X) )+offset_c

# adaptive kernel.
θ_canonical = AdaptiveRKHS.GaussianKernel1DType(1.0/1000)
#warp map.
# d2f_normed = xx->norm(Calculus.hessian(f_joint,xx))
# d2f_x1 = xx->Calculus.hessian(f_joint,xx)[1,1]
# d2f_x2 = xx->Calculus.hessian(f_joint,xx)[2,2]
# d2f_sum = xx->sum(Calculus.hessian(f_joint,xx))

#d2y1 = xx->Calculus.hessian(fq,xx)[1]
#d2y1 = xx->Calculus.hessian(y1_itp,xx)[1]
d2y1 = xx->Calculus.hessian(u1_GP,xx)[1]
#ϕ = d2y1 #getwarpmap(ϕ)
ϕ = d2u1_GP

title_string = "second derivative"
PyPlot.figure(fig_num)
fig_num += 1
#PyPlot.plot(xq, d2u1_GP_ND.(xq), label = "d2u1_GP_ND")
PyPlot.plot(xq, d2u1_GP.(xq), "^", label = "d2u1_GP")
PyPlot.title(title_string)
PyPlot.legend()

# adaptive RKHS.
M = 30
#σ²_RKHS_array = sqrt(σ) .* ones(Float64, M) # #1.0 #1e-1 #1e-3
σ²_RKHS_array = collect( LinRange(sqrt(σ), σ^2, M) )
amp_factor_array = 100.0 .* ones(Float64, M)

fq_array, c_array, θ_array = iterateGP(fq, ϕ, θ_canonical, σ²_RKHS_array, amp_factor_array, y1, X)

#### visualize.

fq_xq = fq.(xq)

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(time_stamp, y1, "o", label = "observed")
PyPlot.plot(xq, fq_xq, label = "canonical RKHS")
PyPlot.plot(xq, fq_array[1].(xq), label = "adaptive RKHS")
PyPlot.plot(xq, fq_array[2].(xq), label = "adaptive RKHS 2")
PyPlot.plot(xq, fq_array[3].(xq), label = "adaptive RKHS 3")
PyPlot.plot(xq, fq_array[end].(xq), label = "adaptive RKHS end")
PyPlot.plot(sol.t, sol_x1_display, "--", label = "species 1")
#PyPlot.plot(sol.t, sol_x2_display, label = "species 2")

title_string = "1-D RKHS demo"
PyPlot.title(title_string)
PyPlot.legend()

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
const Counter = let
    init() = 0
    op(n, _) = n + 1
    value(n) = n
    (; init, op, value)
end

struct FrequencyAnalysis
    options::Dict{Symbol, Any}
end

function (f::FrequencyAnalysis)(entry::Entry)
    plottype = entry.plottype
    mappings = copy(entry.mappings)
    N = length(getvalue(mappings[1]))
    push!(mappings.positional, Labeled("count", fill(nothing, N)))
    attributes = entry.attributes
    return groupreduce(Counter, Entry(plottype, mappings, attributes))
end

"""
    frequency()

Compute a frequency table of the arguments.
"""
frequency() = Layer((FrequencyAnalysis(Dict{Symbol, Any}()),))


text:
@testset "predict" begin
    Random.seed!(100)

    @model function linear_reg(x, y, σ = 0.1)
        β ~ Normal(0, 1)

        for i ∈ eachindex(y)
            y[i] ~ Normal(β * x[i], σ)
        end
    end

    @model function linear_reg_vec(x, y, σ = 0.1)
        β ~ Normal(0, 1)
        y ~ MvNormal(β .* x, σ^2 * I)
    end

    f(x) = 2 * x + 0.1 * randn()

    Δ = 0.1
    xs_train = 0:Δ:10; ys_train = f.(xs_train);
    xs_test = [10 + Δ, 10 + 2 * Δ]; ys_test = f.(xs_test);

    # Infer
    m_lin_reg = linear_reg(xs_train, ys_train);
    chain_lin_reg = sample(m_lin_reg, NUTS(100, 0.65), 200);

    # Predict on two last indices
    m_lin_reg_test = linear_reg(xs_test, fill(missing, length(ys_test)));
    predictions = Turing.Inference.predict(m_lin_reg_test, chain_lin_reg)

    ys_pred = vec(mean(Array(group(predictions, :y)); dims = 1))

    @test sum(abs2, ys_test - ys_pred) ≤ 0.1

    # Ensure that `rng` is respected
    predictions1 = let rng = MersenneTwister(42)
        predict(rng, m_lin_reg_test, chain_lin_reg[1:2])
    end
    predictions2 = let rng = MersenneTwister(42)
        predict(rng, m_lin_reg_test, chain_lin_reg[1:2])
    end
    @test all(Array(predictions1) .== Array(predictions2))

    # Predict on two last indices for vectorized
    m_lin_reg_test = linear_reg_vec(xs_test, missing);
    predictions_vec = Turing.Inference.predict(m_lin_reg_test, chain_lin_reg)
    ys_pred_vec = vec(mean(Array(group(predictions_vec, :y)); dims = 1))

    @test sum(abs2, ys_test - ys_pred_vec) ≤ 0.1

    # Multiple chains
    chain_lin_reg = sample(m_lin_reg, NUTS(100, 0.65), MCMCThreads(), 200, 2);
    m_lin_reg_test = linear_reg(xs_test, fill(missing, length(ys_test)));
    predictions = Turing.Inference.predict(m_lin_reg_test, chain_lin_reg)

    @test size(chain_lin_reg, 3) == size(predictions, 3)

    for chain_idx in MCMCChains.chains(chain_lin_reg)
        ys_pred = vec(mean(Array(group(predictions[:, :, chain_idx], :y)); dims = 1))
        @test sum(abs2, ys_test - ys_pred) ≤ 0.1
    end

    # Predict on two last indices for vectorized
    m_lin_reg_test = linear_reg_vec(xs_test, missing);
    predictions_vec = Turing.Inference.predict(m_lin_reg_test, chain_lin_reg)

    for chain_idx in MCMCChains.chains(chain_lin_reg)
        ys_pred_vec = vec(mean(
            Array(group(predictions_vec[:, :, chain_idx], :y));
            dims = 1
        ))
        @test sum(abs2, ys_test - ys_pred_vec) ≤ 0.1
    end

    # https://github.com/TuringLang/Turing.jl/issues/1352
    @model function simple_linear1(x, y)
        intercept ~ Normal(0,1)
        coef ~ MvNormal(zeros(2), I)
        coef = reshape(coef, 1, size(x,1))

        mu = intercept .+ coef * x |> vec
        error ~ truncated(Normal(0,1), 0, Inf)
        y ~ MvNormal(mu, error^2 * I)
    end;

    @model function simple_linear2(x, y)
        intercept ~ Normal(0,1)
        coef ~ filldist(Normal(0,1), 2)
        coef = reshape(coef, 1, size(x,1))

        mu = intercept .+ coef * x |> vec
        error ~ truncated(Normal(0,1), 0, Inf)
        y ~ MvNormal(mu, error^2 * I)
    end;

    @model function simple_linear3(x, y)
        intercept ~ Normal(0,1)
        coef = Vector(undef, 2)
        for i in axes(coef, 1)
            coef[i] ~ Normal(0,1)
        end
        coef = reshape(coef, 1, size(x,1))

        mu = intercept .+ coef * x |> vec
        error ~ truncated(Normal(0,1), 0, Inf)
        y ~ MvNormal(mu, error^2 * I)
    end;

    @model function simple_linear4(x, y)
        intercept ~ Normal(0,1)
        coef1 ~ Normal(0,1)
        coef2 ~ Normal(0,1)
        coef = [coef1, coef2]
        coef = reshape(coef, 1, size(x,1))

        mu = intercept .+ coef * x |> vec
        error ~ truncated(Normal(0,1), 0, Inf)
        y ~ MvNormal(mu, error^2 * I)
    end;

    # Some data
    x = randn(2, 100);
    y = [1 + 2 * a + 3 * b for (a,b) in eachcol(x)];

    for model in [simple_linear1, simple_linear2, simple_linear3, simple_linear4]
        m = model(x, y);
        chain = sample(m, NUTS(), 100);
        chain_predict = predict(model(x, missing), chain);
        mean_prediction = [chain_predict["y[$i]"].data |> mean for i = 1:length(y)]
        @test mean(abs2, mean_prediction - y) ≤ 1e-3
    end
end

@testset "Timer" begin
    chain = sample(gdemo_default, MH(), 1000)

    @test chain.info.start_time isa Float64
    @test chain.info.stop_time isa Float64
    @test chain.info.start_time < chain.info.stop_time
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to IAU-76/FK5 nutation algorithm.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/transformations/fk5/nutation.jl
# ===========================================

# Function nutation_fk5
# ---------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# In this example, using JD_TT = 2453101.828154745, one gets the following data
# related to the nutation:
#
#   mϵ_1980 =  23.4387368˚
#   Δϵ_1980 =   0.0020316˚
#   Δψ_1980 =  -0.0034108˚
#
# SatelliteToolbox provides the following results:
#
#   julia> (mϵ_1980, Δϵ_1980, Δψ_1980) = nutation_fk5(2453101.828154745)
#   (0.40908312815877673, 3.5458763448549555e-5, -5.953027070867465e-5)
#
#   julia> mϵ_1980*180/pi
#   23.438736713507268
#
#   julia> Δϵ_1980*180/pi
#   0.0020316374923546377
#
#   julia> Δψ_1980*180/pi
#   -0.0034108332648783257
#
# Scenario 02
# ===========
#
# Example in [1, pp. 233 and 234].
#
# In the day 182.78495062 of year 2000, if one uses only 4 terms in nutation
# computation (IAU-76/FK5 model), then one gets:
#
#   Δψ_1980 = -0.004250260˚
#   Δϵ_1980 = -0.001260854˚
#   mϵ_1980 = +23.43922657˚
#
# If the same algorithm is applied to the day 179.78495062 of year 2000, the one
# gets:
#
#   Δψ_1980 = -0.004337544˚
#   Δϵ_1980 = -0.001247061˚
#   mϵ_1980 = +23.43922657˚
#
################################################################################

@testset "Nutation" begin
    # Scenario 01
    # ===========

    mϵ_1980, Δϵ_1980, Δψ_1980 = nutation_fk5(2453101.828154745)

    @test mϵ_1980 * 180 / π ≈ 23.4387368 atol=1e-7
    @test Δϵ_1980 * 180 / π ≈  0.0020316 atol=1e-7
    @test Δψ_1980 * 180 / π ≈ -0.0034108 atol=1e-7

    # Scenario 02
    # ===========

    JD_UTC = date_to_jd(2000, 1, 1, 0, 0, 0) + 182.78495062 - 1
    JD_TT  = jd_utc_to_tt(JD_UTC)

    (mϵ_1980, Δϵ_1980, Δψ_1980) = nutation_fk5(JD_TT, 4)

    @test mϵ_1980 * 180 / π ≈ 23.43922657  atol=1e-6
    @test Δϵ_1980 * 180 / π ≈ -0.001260854 atol=1e-8
    @test Δψ_1980 * 180 / π ≈ -0.004250260 atol=1e-8

    JD_UTC = date_to_jd(2000, 1, 1, 0, 0, 0) + 179.78495062 - 1
    JD_TT  = jd_utc_to_tt(JD_UTC)

    mϵ_1980, Δϵ_1980, Δψ_1980 = nutation_fk5(JD_TT, 4)

    @test mϵ_1980 * 180 / π ≈ 23.43922657  atol=1e-6
    @test Δϵ_1980 * 180 / π ≈ -0.001247061 atol=1e-8
    @test Δψ_1980 * 180 / π ≈ -0.004337544 atol=1e-8

    # Warnings
    # ========
    JD_TT = 2453101.828154745
    mϵ_1980, Δϵ_1980, Δψ_1980 = nutation_fk5(JD_TT)

    @test (@test_logs(
        (
            :warn,
            "The maximum number of coefficients to compute nutation using IAU-76/FK5 theory is 106."
        ),
        nutation_fk5(JD_TT, 107)
    )) == (mϵ_1980, Δϵ_1980, Δψ_1980)

    @test (@test_logs(
        (
            :warn,
            "n_max must greater than 0. The default value will be used (106)."
        ),
        nutation_fk5(JD_TT, 0)
    )) == (mϵ_1980, Δϵ_1980, Δψ_1980)
end


text:
module test_bukdu_plugs_parsers

using Test
using Bukdu
import .Bukdu: Deps
import .Bukdu.Plug: Parsers

struct FC <: ApplicationController; conn::Conn; end

function urlencoded(c::FC)
    c.params.q
end

function multipart(c::FC)
    c.params
end

routes() do
    post("/urlencoded", FC, urlencoded)
    post("/multipart", FC, multipart)
end

req = Deps.HTTP.Request(
    "POST",
    "/urlencoded",
    ["Content-Type"=>"application/x-www-form-urlencoded", "Content-Length"=>"3"],
    codeunits("q=5")
)
@test Router.call(req).got == "5"

req = Deps.HTTP.Request(
    "POST",
    "/multipart",
    ["Content-Type"=>"multipart/form-data; boundary=boundary"],
    codeunits("""
--boundary
Content-Disposition: form-data; name="field1"

value1
--boundary
Content-Disposition: form-data; name="field2"; filename="example.txt"

value2
--boundary--""")
)
@test Router.call(req).got == Assoc("field1"=>"value1", "field2"=>"value2")

req = Deps.HTTP.Request(
    "POST",
    "/multipart",
    ["Content-Type"=>"multipart/form-data; boundary=boundary"],
    codeunits("""
--boundary
Content-Disposition: form-data; name="field1"

value1
--boundary
Content-Disposition: form-data; name="field2"; filename="example.txt"
Content-Type: text/plain

value2
--boundary--""")
)
@test Router.call(req).got == Assoc("field1"=>"value1", "field2"=>"value2")

req = Deps.HTTP.Request(
    "POST",
    "/multipart",
    ["Content-Type"=>"multipart/form-data; boundary=----WebKitFormBoundaryLeFS8sRvPzeszrWy"],
    codeunits("""
------WebKitFormBoundaryLeFS8sRvPzeszrWy\r\nContent-Disposition: form-data; name=\"user_name\"\r\n\r\n\r\n------WebKitFormBoundaryLeFS8sRvPzeszrWy\r\nContent-Disposition: form-data; name=\"user_famous\"\r\n\r\nfalse\r\n------WebKitFormBoundaryLeFS8sRvPzeszrWy--\r\n""")
)
@test Router.call(req).got == Assoc("user_name"=>"", "user_famous"=>"false")

Routing.empty!()

@test  Parsers.rstripcrlf([0x0d, 0x0a]) == []      # CR LF
@test  Parsers.rstripcrlf([0x0d])       == [0x0d]  # CR

end # module test_bukdu_plugs_parsers

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
function show_Hello(io::IO, x)
    print(io, "Hello, ")
    show(io, x)
    print(io, '!')
end
show_Hello(x) = show_Hello(stdout, x)

show_Hello("Julia")

# %%
str = sprint(show_Hello, "Julia")

# %%
print(str)

# %%
?sprint

# %%
methods(sprint)

# %%

max_stars_repo_path: 0015/sprint(f, x) example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/sprint(f, x) example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/sprint(f, x) example.jl

max_forks_repo_name: genkuroki/public


text:
# # Extending Register: an echo register

# We show how to extend a register that prints the instructions to standard output
# instead of actually exexcuting the instructions.

using Yao

# First you will need to define your own subtype that contains
# your desired data structure

mutable struct EchoReg{B} <: AbstractRegister{B}
    nactive::Int
    nqubits::Int
end

# Then we need to define the most basic API: how many qubits are there

Yao.nactive(reg::EchoReg) = reg.nactive
Yao.nqubits(reg::EchoReg) = reg.nqubits

# And define some instructions, as an echo register, we will just keep printing what we are asked to exexcute

function Yao.instruct!(::EchoReg{B}, args...) where {B, G}
    str = join(string.(args), ", ")
    println("calls: instruct!(reg, $str)")
end

function Yao.focus!(reg::EchoReg{B}, locs) where {B}
    println("focus -> $locs")
    reg.nactive = length(locs)
    return true
end

function Yao.relax!(reg::EchoReg{B}, locs; to_nactive=nqubits(reg)) where {B}
    reg.nactive = to_nactive
    println("relax -> $locs\\$to_nactive")
    return true
end

function Yao.measure!(rng, ::ComputationalBasis, reg::EchoReg{B}, locs) where {B}
    println("measure -> $locs")
    return true
end

# Now we can check what will happen!
r = EchoReg{10}(3, 2)
r |> put(3, 2=>X) |> control(3, 3, 2=>X) |> concentrate(3, put(1, 1=>X), 2:2) |> measure!

# But when we still require some more information from you

r |> cache(X)

# An `NotImplementedError` will be thrown, no worries!