text:
"""
    ArrayAllocators.Windows

Defines array allocators on Windows.

# Examples

```julia
using ArrayAllocators.Windows

Array{UInt8}(WinMemAlign(2^16), 1024)
```
"""
module Windows

import ..ArrayAllocators: AbstractArrayAllocator, nbytes, allocate
import ..ArrayAllocators: AbstractMemAlign, min_alignment, alignment
import ..ArrayAllocators: iszeroinit, lineage_finalizer
import Base: Array

export WinMemAlign

@static if Sys.iswindows()
    const hCurrentProcess = ccall((:GetCurrentProcess, "kernel32"), Ptr{Nothing}, ())
else
    # The value below is the typical default value on Windows when executing the above
    const hCurrentProcess = Ptr{Nothing}(0xffffffffffffffff)
end

const kernel32 = "kernel32"
const kernelbase = "kernelbase"
const MEM_COMMIT      = 0x00001000
const MEM_RESERVE     = 0x00002000
const MEM_RESET       = 0x00080000
const MEM_RESET_UNDO  = 0x10000000
const MEM_LARGE_PAGES = 0x20000000
const MEM_PHYSICAL    = 0x00400000
const MEM_TOP_DOWN    = 0x00100000

const MEM_COMMIT_RESERVE = MEM_COMMIT | MEM_RESERVE

const PAGE_READWRITE = 0x04

const MEM_DECOMMIT = 0x00004000
const MEM_RELEASE  = 0x00008000

const DWORD = Culong

abstract type MemExtendedParameterType end

"""
    MemAddressRequirements

See https://docs.microsoft.com/en-us/windows/win32/api/winnt/ns-winnt-mem_address_requirements
"""
struct MemAddressRequirements
    lowestStartingAddress::Ptr{Nothing}
    highestStartingAddress::Ptr{Nothing}
    alignment::Csize_t
end
MemAddressRequirements(alignment) = MemAddressRequirements(C_NULL, C_NULL, alignment)

"""
    MemExtendedParameterAddressRequirements

This is a Julian structure where the `requirements` field is `Base.RefValue{MemAddressRequirements}`

See https://docs.microsoft.com/en-us/windows/win32/api/winnt/ns-winnt-mem_extended_parameter

See also `_MemExtendedParameterAddressRequirements`
"""
struct MemExtendedParameterAddressRequirements <: MemExtendedParameterType
    type::UInt64
    requirements::Base.RefValue{MemAddressRequirements}
    function MemExtendedParameterAddressRequirements(requirements)
        new(1, Ref(requirements))
    end
end

#=
    _MemExtendedParameterAddressRequirements

Internal structure compatible with C definition

See https://docs.microsoft.com/en-us/windows/win32/api/winnt/ns-winnt-mem_extended_parameter
=#
struct _MemExtendedParameterAddressRequirements
    type::UInt64
    pointer::Ptr{MemAddressRequirements}
end

function Base.convert(::Type{_MemExtendedParameterAddressRequirements}, p::MemExtendedParameterAddressRequirements)
    _MemExtendedParameterAddressRequirements(p.type, pointer_from_objref(p.requirements))
end

#=
LPVOID VirtualAllocEx(
  [in]           HANDLE hProcess,
  [in, optional] LPVOID lpAddress,
  [in]           SIZE_T dwSize,
  [in]           DWORD  flAllocationType,
  [in]           DWORD  flProtect
);
=#
function VirtualAllocEx(hProcess, lpAddress, dwSize, flAllocationType, flProtect)
    ccall((:VirtualAllocEx, kernel32), Ptr{Nothing},
        (Ptr{Nothing}, Ptr{Nothing}, Csize_t, DWORD, DWORD),
        hProcess, lpAddress, dwSize, flAllocationType, flProtect
        )
end
function VirtualAllocEx(dwSize)
    VirtualAllocEx(hCurrentProcess, C_NULL, dwSize, MEM_COMMIT_RESERVE, PAGE_READWRITE)
end

#=
PVOID VirtualAlloc2(
  [in, optional]      HANDLE                 Process,
  [in, optional]      PVOID                  BaseAddress,
  [in]                SIZE_T                 Size,
  [in]                ULONG                  AllocationType,
  [in]                ULONG                  PageProtection,
  [in, out, optional] MEM_EXTENDED_PARAMETER *ExtendedParameters,
  [in]                ULONG                  ParameterCount
);
=#
# https://docs.microsoft.com/en-us/windows/win32/api/memoryapi/nf-memoryapi-virtualalloc2
function VirtualAlloc2(Process, BaseAddress, Size, AllocationType, PageProtection, ExtendedParameters, ParameterCount)
    # Docs say
    ccall((:VirtualAlloc2, kernelbase), Ptr{Nothing},
          (Ptr{Nothing}, Ptr{Nothing}, Csize_t, Culong, Culong, Ptr{Nothing}, Culong),
          Process, BaseAddress, Size, AllocationType, PageProtection, ExtendedParameters, ParameterCount
         )
end
function VirtualAlloc2(Process, BaseAddress, Size, AllocationType, PageProtection, ParameterCount)
    VirtualAlloc2(Process, BaseAddress, Size, AllocationType, PageProtection, C_NULL, ParameterCount)
end

function win_memalign(alignment, num_bytes; lowestStartingAddress::Ptr{Nothing} = C_NULL, highestStartingAddress::Ptr{Nothing} = C_NULL)
    reqs = MemAddressRequirements(lowestStartingAddress, highestStartingAddress, alignment)
    pr = MemExtendedParameterAddressRequirements(reqs)
    p = Ref{_MemExtendedParameterAddressRequirements}(pr)
    GC.@preserve reqs pr p begin
        ptr = VirtualAlloc2(hCurrentProcess, C_NULL, num_bytes, MEM_COMMIT_RESERVE, PAGE_READWRITE, p, 1)
    end
    return ptr
end

const WIN_MEMALIGN_MIN_ALIGNMENT = 2^16
function check_alignment(alignment)
    ispow2(alignment) || throw(ArgumentError("Alignment must be a power of 2"))
    alignment ≥ WIN_MEMALIGN_MIN_ALIGNMENT || throw(ArgumentError("Alignment must be a multiple of $(WIN_MEMALIGN_MIN_ALIGNMENT)"))
    return nothing
end

#=
BOOL VirtualFreeEx(
  [in] HANDLE hProcess,
  [in] LPVOID lpAddress,
  [in] SIZE_T dwSize,
  [in] DWORD  dwFreeType
);
=#

const BOOL = Cint

function VirtualFreeEx(hProcess, lpAddress, dwSize, dwFreeType)
    ccall((:VirtualFreeEx, kernel32),
        BOOL, (Ptr{Nothing}, Ptr{Nothing}, Csize_t, DWORD),
        hProcess, lpAddress, dwSize, dwFreeType
    )
end
function VirtualFreeEx(lpAddress)
    VirtualFreeEx(hCurrentProcess, lpAddress, 0, MEM_RELEASE)
end

function virtual_free(array::Array{T}) where T
    VirtualFreeEx(array)
end

function wrap_virtual(::Type{A}, ptr::Ptr{T}, dims) where {T, A <: AbstractArray{T}}
    if ptr == C_NULL
        throw(OutOfMemoryError())
    end
    arr = unsafe_wrap(A, ptr, dims; own = false)
    lineage_finalizer(virtual_free, arr)
    return arr
end
wrap_virtual(ptr::Ptr{T}, dims) where T = wrap_virtual(Array{T}, ptr, dims)

# == AbstractWinVirtualAllocator == #

abstract type AbstractWinVirtualAllocator{B} <: AbstractArrayAllocator{B} end
allocate(::AbstractWinVirtualAllocator, num_bytes) =  VirtualAllocEx(num_bytes)
Base.unsafe_wrap(::AbstractWinVirtualAllocator, args...) = wrap_virtual(args...)
iszeroinit(::Type{A}) where A <: AbstractWinVirtualAllocator = true

# == WinVirtualAllocator == #

struct WinVirtualAllocator{B} <: AbstractWinVirtualAllocator{B}
end

const virtual = WinVirtualAllocator()

# == WindowsMemAlign == #

"""
    WinMemAlign([alignment, lowestStartingAddress, highestStartingAddress])

Uses `VirtualAlloc2` to allocate aligned memory. `alignment` must be a power of 2 and larger than $(WIN_MEMALIGN_MIN_ALIGNMENT).
"""
struct WinMemAlign{B} <: AbstractMemAlign{B}
    alignment::Int
    lowestStartingAddress::Ptr{Nothing}
    highestStartingAddress::Ptr{Nothing}
    function WinMemAlign{B}(alignment, lowestStartingAddress = C_NULL, highestStartingAddress = C_NULL) where B
        check_alignment(alignment)
        return new{B}(
            alignment,
            reinterpret(Ptr{Nothing}, lowestStartingAddress),
            reinterpret(Ptr{Nothing}, highestStartingAddress)
        )
    end
end
WinMemAlign() = WinMemAlign(WIN_MEMALIGN_MIN_ALIGNMENT)
allocate(alloc::WinMemAlign, num_bytes) = win_memalign(
    alloc.alignment,
    num_bytes;
    lowestStartingAddress = alloc.lowestStartingAddress,
    highestStartingAddress = alloc.highestStartingAddress
)
Base.unsafe_wrap(::WinMemAlign, args...) =  wrap_virtual(args...)
min_alignment(::Type{WinMemAlign}) = WIN_MEMALIGN_MIN_ALIGNMENT
iszeroinit(::Type{A}) where A <: WinMemAlign = true

end # module Windows


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

adjoint(a::AbstractArray) = error("adjoint not defined for $(typeof(a)). Consider using `permutedims` for higher-dimensional arrays.")
transpose(a::AbstractArray) = error("transpose not defined for $(typeof(a)). Consider using `permutedims` for higher-dimensional arrays.")

## Matrix transposition ##

"""
    transpose!(dest,src)

Transpose array `src` and store the result in the preallocated array `dest`, which should
have a size corresponding to `(size(src,2),size(src,1))`. No in-place transposition is
supported and unexpected results will happen if `src` and `dest` have overlapping memory
regions.

# Examples
```jldoctest
julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> B = zeros(Complex{Int64}, 2, 2)
2×2 Array{Complex{Int64},2}:
 0+0im  0+0im
 0+0im  0+0im

julia> transpose!(B, A);

julia> B
2×2 Array{Complex{Int64},2}:
 3+2im  8+7im
 9+2im  4+6im

julia> A
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im
```
"""
transpose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(transpose, B, A)

"""
    adjoint!(dest,src)

Conjugate transpose array `src` and store the result in the preallocated array `dest`, which
should have a size corresponding to `(size(src,2),size(src,1))`. No in-place transposition
is supported and unexpected results will happen if `src` and `dest` have overlapping memory
regions.

# Examples
```jldoctest
julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> B = zeros(Complex{Int64}, 2, 2)
2×2 Array{Complex{Int64},2}:
 0+0im  0+0im
 0+0im  0+0im

julia> adjoint!(B, A);

julia> B
2×2 Array{Complex{Int64},2}:
 3-2im  8-7im
 9-2im  4-6im

julia> A
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im
```
"""
adjoint!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(adjoint, B, A)
function transpose!(B::AbstractVector, A::AbstractMatrix)
    axes(B,1) == axes(A,2) && axes(A,1) == 1:1 || throw(DimensionMismatch("transpose"))
    copyto!(B, A)
end
function transpose!(B::AbstractMatrix, A::AbstractVector)
    axes(B,2) == axes(A,1) && axes(B,1) == 1:1 || throw(DimensionMismatch("transpose"))
    copyto!(B, A)
end
function adjoint!(B::AbstractVector, A::AbstractMatrix)
    axes(B,1) == axes(A,2) && axes(A,1) == 1:1 || throw(DimensionMismatch("transpose"))
    ccopy!(B, A)
end
function adjoint!(B::AbstractMatrix, A::AbstractVector)
    axes(B,2) == axes(A,1) && axes(B,1) == 1:1 || throw(DimensionMismatch("transpose"))
    ccopy!(B, A)
end

const transposebaselength=64
function transpose_f!(f, B::AbstractMatrix, A::AbstractMatrix)
    inds = axes(A)
    axes(B,1) == inds[2] && axes(B,2) == inds[1] || throw(DimensionMismatch(string(f)))

    m, n = length(inds[1]), length(inds[2])
    if m*n<=4*transposebaselength
        @inbounds begin
            for j = inds[2]
                for i = inds[1]
                    B[j,i] = f(A[i,j])
                end
            end
        end
    else
        transposeblock!(f,B,A,m,n,first(inds[1])-1,first(inds[2])-1)
    end
    return B
end
function transposeblock!(f, B::AbstractMatrix, A::AbstractMatrix, m::Int, n::Int, offseti::Int, offsetj::Int)
    if m*n<=transposebaselength
        @inbounds begin
            for j = offsetj .+ (1:n)
                for i = offseti .+ (1:m)
                    B[j,i] = f(A[i,j])
                end
            end
        end
    elseif m>n
        newm=m>>1
        transposeblock!(f,B,A,newm,n,offseti,offsetj)
        transposeblock!(f,B,A,m-newm,n,offseti+newm,offsetj)
    else
        newn=n>>1
        transposeblock!(f,B,A,m,newn,offseti,offsetj)
        transposeblock!(f,B,A,m,n-newn,offseti,offsetj+newn)
    end
    return B
end

function ccopy!(B, A)
    RB, RA = eachindex(B), eachindex(A)
    if RB == RA
        for i = RB
            B[i] = adjoint(A[i])
        end
    else
        for (i,j) = zip(RB, RA)
            B[i] = adjoint(A[j])
        end
    end
    return B
end

"""
    copy(A::Transpose)
    copy(A::Adjoint)

Eagerly evaluate the lazy matrix transpose/adjoint.
Note that the transposition is applied recursively to elements.

This operation is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims), which is non-recursive.

# Examples
```jldoctest
julia> A = [1 2im; -3im 4]
2×2 Array{Complex{Int64},2}:
 1+0im  0+2im
 0-3im  4+0im

julia> T = transpose(A)
2×2 Transpose{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0-3im
 0+2im  4+0im

julia> copy(T)
2×2 Array{Complex{Int64},2}:
 1+0im  0-3im
 0+2im  4+0im
```
"""
copy(::Union{Transpose,Adjoint})

Base.copy(A::Transpose{<:Any,<:AbstractMatrix}) = transpose!(similar(A.parent, reverse(axes(A.parent))), A.parent)
Base.copy(A::Adjoint{<:Any,<:AbstractMatrix}) = adjoint!(similar(A.parent, reverse(axes(A.parent))), A.parent)

function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int},
                         A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int})
    if length(ir_dest) != length(jr_src)
        throw(ArgumentError(string("source and destination must have same size (got ",
                                   length(jr_src)," and ",length(ir_dest),")")))
    end
    if length(jr_dest) != length(ir_src)
        throw(ArgumentError(string("source and destination must have same size (got ",
                                   length(ir_src)," and ",length(jr_dest),")")))
    end
    @boundscheck checkbounds(B, ir_dest, jr_dest)
    @boundscheck checkbounds(A, ir_src, jr_src)
    idest = first(ir_dest)
    for jsrc in jr_src
        jdest = first(jr_dest)
        for isrc in ir_src
            B[idest,jdest] = A[isrc,jsrc]
            jdest += step(jr_dest)
        end
        idest += step(ir_dest)
    end
    return B
end

ext: jl

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


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

mutable struct IoK8sApiCoreV1WeightedPodAffinityTerm <: SwaggerModel
    podAffinityTerm::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodAffinityTerm } # spec name: podAffinityTerm
    weight::Any # spec type: Union{ Nothing, Int32 } # spec name: weight

    function IoK8sApiCoreV1WeightedPodAffinityTerm(;podAffinityTerm=nothing, weight=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1WeightedPodAffinityTerm, Symbol("podAffinityTerm"), podAffinityTerm)
        setfield!(o, Symbol("podAffinityTerm"), podAffinityTerm)
        validate_property(IoK8sApiCoreV1WeightedPodAffinityTerm, Symbol("weight"), weight)
        setfield!(o, Symbol("weight"), weight)
        o
    end
end # type IoK8sApiCoreV1WeightedPodAffinityTerm

const _property_map_IoK8sApiCoreV1WeightedPodAffinityTerm = Dict{Symbol,Symbol}(Symbol("podAffinityTerm")=>Symbol("podAffinityTerm"), Symbol("weight")=>Symbol("weight"))
const _property_types_IoK8sApiCoreV1WeightedPodAffinityTerm = Dict{Symbol,String}(Symbol("podAffinityTerm")=>"IoK8sApiCoreV1PodAffinityTerm", Symbol("weight")=>"Int32")
Base.propertynames(::Type{ IoK8sApiCoreV1WeightedPodAffinityTerm }) = collect(keys(_property_map_IoK8sApiCoreV1WeightedPodAffinityTerm))
Swagger.property_type(::Type{ IoK8sApiCoreV1WeightedPodAffinityTerm }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1WeightedPodAffinityTerm[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1WeightedPodAffinityTerm }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1WeightedPodAffinityTerm[property_name]

function check_required(o::IoK8sApiCoreV1WeightedPodAffinityTerm)
    (getproperty(o, Symbol("podAffinityTerm")) === nothing) && (return false)
    (getproperty(o, Symbol("weight")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1WeightedPodAffinityTerm }, name::Symbol, val)
end


text:
using PyPlot

# This script relies on data (emu_out) that is generated by running smc-abc-script.jl

function plot_distances_training_accepted(emu_out::ABCSMCOutput)
    n_iter = length(emu_out.threshold_schedule)
    for smc_iteration in 1:n_iter
        for parameter_idx in 1:emu_out.n_params
            subplot_idx = (smc_iteration - 1) * emu_out.n_params + parameter_idx
            subplot(n_iter, emu_out.n_params, subplot_idx)
            title("SMC iteration $(smc_iteration)")
            scatter(emu_out.emulators[smc_iteration].gp_training_x[:, parameter_idx],
                emu_out.emulators[smc_iteration].gp_training_y, marker=".",
                label="Emulation - training")
            scatter(emu_out.population[smc_iteration][:, parameter_idx],
                emu_out.distances[smc_iteration], marker=".",
                label="Emulation - accepted")
            xlabel("Theta_$(parameter_idx)")
            ylabel("SMC iteration $(smc_iteration) - L2 distance")
            legend()
        end
    end
end

ion()
f = figure()
ioff()
plot_distances_training_accepted(emu_out)
subplots_adjust(
left    =  0.08,
bottom  =  0.06,
right   =  0.96,
top     =  0.97,
wspace  =  0.26,
hspace  =  0.3
)
show(f)

savefig("/Users/tanhevg/Desktop/l2-norm.pdf")


text:
using ContMechTensors
using BenchmarkTools
using JLD

include("generate_report.jl")

const SUITE = BenchmarkGroup()

function create_tensors()
    tensor_dict = Dict{Tuple{Int, Int, DataType}, AbstractTensor}()
    symtensor_dict = Dict{Tuple{Int, Int, DataType}, AbstractTensor}()
    for dim in 1:3
        for order in (1,2,4)
            for T in (Float32, Float64)
                tensor_dict[(dim, order, T)] = rand(Tensor{order, dim, T})
                if order != 1
                    symtensor_dict[(dim, order, T)] = rand(SymmetricTensor{order, dim, T})
                else
                    symtensor_dict[(dim, order, T)] = rand(Tensor{order, dim, T})
                end
            end
        end
    end
    return tensor_dict, symtensor_dict
end

tensor_dict, symtensor_dict = create_tensors()

include("benchmark_functions.jl")
include("benchmark_ad.jl")

function run_benchmarks(name, tagfilter = @tagged ALL)
    const paramspath = joinpath(dirname(@__FILE__), "params.jld")
    if !isfile(paramspath)
        println("Tuning benchmarks...")
        tune!(SUITE, verbose=true)
        JLD.save(paramspath, "SUITE", params(SUITE))
    end
    loadparams!(SUITE, JLD.load(paramspath, "SUITE"), :evals, :samples)
    results = run(SUITE[tagfilter], verbose = true, seconds = 2)
    JLD.save(joinpath(dirname(@__FILE__), name * ".jld"), "results", results)
end

function generate_report(v1, v2)
    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
    v2_res = load(joinpath(dirname(@__FILE__), v2 * ".jld"), "results")
    open(joinpath(dirname(@__FILE__), "results_$(v1)_$(v2).md"), "w") do f
        printreport(f, judge(minimum(v1_res), minimum(v2_res)); iscomparisonjob = true)
    end
end

function generate_report(v1)
    v1_res = load(joinpath(dirname(@__FILE__), v1 * ".jld"), "results")
    open(joinpath(dirname(@__FILE__), "results_$(v1).md"), "w") do f
        printreport(f, minimum(v1_res); iscomparisonjob = false)
    end
end


text:
function test_InvariantDistribution(; parallel = false, print = false)
  ts3::Vector{Float64} = [1, 5, 4, 2, 3.2, -5, -6, -7, -4, 4, 5, -1, -2, 3, 2, 1, 5, 4, 2, 3.2,
  -5, -6, -7, -4, 4, 5, -1, -2, 3, 2, 1, 5, 4, 2, 3.2, -5, -6, -7, -4, 4, 5,
  -1, -2, 3, 2, 1, 5, 4, 2, 3.2, -5, -6, -7, -4, 4, 5, -1, -2, 3, 2]

  if parallel
    SS = markovmatrix_parallel(ts3, 3, 2, print = print)
  elseif !parallel
    SS = markovmatrix(ts3, 3, 2, print = print)
  end

  invmeasure = invariantdist(SS.P)

  # Invariant distribution
  invdist = invmeasure[1]

  # Reference to the simplices with stritly positive measure (column vector of
  # indices referring to the rows of the triangulation, each corresponding to
  # a simplex).
  positive_measure_simplex_inds = invmeasure[2]

  # Get the simplices with strictly positive measure
  simplex_inds = SS.triangulation_indices[positive_measure_simplex_inds, :]

  return(invdist)
end

@test sum(test_InvariantDistribution()) ≈ 1.0 atol = 1/10^7
@time test_InvariantDistribution()


text:
#######
#additional functions to compute a kernel on a graph.
#####

#####
#Define the edge rules, resulting in the weight on the graph between edges
#####
"""
    EdgeRule

Determines the value for the edge.
"""
abstract type EdgeRule end

# define the type of the edge, default is Boolean
Base.eltype(::EdgeRule) = Float64

"""
    NCommon(n::Int)
Edgerule,evaluteds if there are at least n tiles are in common.

"""

struct NCommon <: EdgeRule
    n::Int
end

OneCommon() = NCommon(1)
TwoCommon() = NCommon(2)

"""
    edgevalue(rule::NCommon, sᵢ::Set, sⱼ::Set)

Return 1 if the two sets; sᵢ, sⱼ, have minimum `n` equal modules, 0 if not.
`n` obtained form `Ncommon` struct
"""
edgevalue(rule::NCommon, sᵢ::Set, sⱼ::Set) = length(sᵢ ∩ sⱼ) ≥ rule.n |> Int

"""
    edgevalue(rule::NCommon, xᵢ, xⱼ)

Return 1 if the two vectors; xᵢ, xⱼ,  have minimum `n` equal modules, 0 if not.
`n` obtained form `Ncommon` struct
"""

edgevalue(rule::NCommon, xᵢ, xⱼ) = edgevalue(rule, Set(xᵢ), Set(xⱼ))

"""
    LevRule
Edgerule, base on the Levenshtein distance between strings
"""
struct LevRule <: EdgeRule end

"""
    edgevalue(rule::LevRule, xᵢ, xⱼ)

Returns ``1-\\dfrac{d_\\text{lev}(x_i,x_j)}{\\text{max}}``.
Where ``d_\\text{lev}(x_i,x_j)`` is levenstein distance between xᵢ, xⱼ.
and max is the maximum levenstein distance between xᵢ, xⱼ possible.
"""

function edgevalue(rule::LevRule, xᵢ, xⱼ)
    max = maximum([length(xᵢ),length(xⱼ)])
    return (1 - evaluate(Levenshtein(), xᵢ, xⱼ)/max)
end

"""
    CosRule(q::Int)
Edgerule, base on the Cosine distance between strings

"""
struct CosRule <: EdgeRule  end

"""
    edgevalue(rule::CosRule, xᵢ, xⱼ)
Returns the Cosine similarity between xᵢ, xⱼ.
Calulation base on the Cosine function of
[StringDistances.jl](https://github.com/matthieugomez/StringDistances.jl)
current implementation, q = 1 if fixed.
Extention to use the full capacity of this function can be made
"""

edgevalue(rule::CosRule, xᵢ, xⱼ) = (1 - evaluate(StringDistances.Cosine(1), xᵢ, xⱼ))

#######
#Function of the graph
#######
"""
    Base.exp(E::Eigen)
Efficient exponential of a matrix based on the Eigenvalue decomposition
"""
Base.exp(E::Eigen) = eigvecs(E) * Diagonal(exp.(eigvals(E))) * eigvecs(E)'
Base.:*(a::Real, E::Eigen) = Eigen(a * eigvals(E), eigvecs(E))

"""
    degree(A::AbstractMatrix)
Returns a vector with the degree of the graph
"""

degree(A::AbstractMatrix) = Diagonal(sum(A, dims=1)[:])

"""
laplacian(A::AbstractMatrix)
Constructs a Laplacian based on the adjacency matrix `A`.
"""
laplacian(A::AbstractMatrix) = degree(A) - A

"""
    diffusion(L::AbstractMatrix, β::Real=1.0)

Compute the diffusion kernel directly on the Laplacian.
"""
diffusion(L::AbstractMatrix, β::Real=1.0) = exp(β * L)

"""
    diffusion(E::Eigen, β::Real=1.0)

Compute the diffusion kernel based on the Eigenvalue decomposition
of the Laplacian.
"""
diffusion(E::Eigen, β::Real=1.0) = exp(β * E)

"""
    norm_laplacian(L)

Normalizes a Laplacian matrix such that the diagonal values equal -1.
"""
norm_laplacian(A) = I - (degree(A)^(-0.5))*A*(degree(A)^(-0.5))

"""
adjacency(x::AbstractVector, edgerule::EdgeRule)

Constructs an adjacency matrix of the graph based on a list of nodes `x` and and the `edgerule`.
The Diagnal is set to zero, no self-loops allowed in the graph
"""

function adjacency(X::AbstractVector, d::EdgeRule)
    n = length(X)
    A = Matrix{Float64}(undef, n, n)
    for (i, xᵢ) in enumerate(X)
        for (j, xⱼ) in enumerate(X)
            A[i,j] = i==j ?  0 : edgevalue(d,xᵢ, xⱼ)
        end
    end
    return A
end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
# Broyden
## NOTE: For the time being it is better to use `LimitedMemoryBroydenSolver`
struct BroydenCache{J,F,X}
    Jinv::J
    fx::F
    Δfx::F
    fx_old::F
    x::X
    Δx::X
    x_old::X
end

function BroydenCache(x)
    fx, Δfx, fx_old = copy(x), copy(x), copy(x)
    x, Δx, x_old = copy(x), copy(x), copy(x)
    Jinv = _init_identity_matrix(x)
    return BroydenCache(Jinv, fx, Δfx, fx_old, x, Δx, x_old)
end

BroydenCache(vec_length::Int, device) = BroydenCache(device(zeros(vec_length)))

"""
    BroydenSolver(; T=Float32, device, original_dims::Tuple{Int,Int}, batch_size, maxiters::Int=50, ϵ::Real=1e-6,
                  abstol::Union{Real,Nothing}=nothing, reltol::Union{Real,Nothing}=nothing)

Broyden Solver ([broyden1965class](@cite)) for solving Discrete DEQs. It is recommended to use [`LimitedMemoryBroydenSolver`](@ref) for better performance.

## Arguments

* `T`: The type of the elements of the vectors. (Default: `Float32`)
* `device`: The device to use. Pass `gpu` to use the GPU else pass `cpu`.
* `original_dims`: Dimensions to reshape the arrays into (excluding the batch dimension).
* `batch_size`: The batch size of the problem. Your inputs can have a different batch size, but having
                them match allows us to efficiently cache internal statistics without reallocation.
* `maxiters`: Maximum number of iterations to run.
* `ϵ`: Tolerance for convergence.
* `abstol`: Absolute tolerance.
* `reltol`: Relative tolerance. (This value is ignored by `BroydenSolver` at the moment)

See also: [`LimitedMemoryBroydenSolver`](@ref)
"""
struct BroydenSolver{C<:BroydenCache,T<:Real}
    cache::C
    maxiters::Int
    batch_size::Int
    ϵ::T
end

function BroydenSolver(; T=Float32, device, original_dims::Tuple{Int,Int}, batch_size, maxiters::Int=50, ϵ::Real=1e-6,
                       abstol::Union{Real,Nothing}=nothing, reltol::Union{Real,Nothing}=nothing)
    ϵ = abstol !== nothing ? abstol : ϵ

    if reltol !== nothing
        @warn "reltol is set to $reltol, but `BroydenSolver` ignores this value" maxlog=1
    end

    x = device(zeros(T, prod(original_dims) * batch_size))
    cache = BroydenCache(x)

    return BroydenSolver(cache, maxiters, batch_size, T(ϵ))
end

function (broyden::BroydenSolver{C,T})(f!, x_::AbstractVector{T}) where {C,T}
    @unpack Jinv, fx, Δfx, fx_old, x, Δx, x_old = broyden.cache
    if size(x) != size(x_)
        # This might happen when the last batch with insufficient batch_size
        # is passed.
        @unpack Jinv, fx, Δfx, fx_old, x, Δx, x_old = BroydenCache(x_)
    end
    x .= x_

    f!(fx, x)
    _init_identity_matrix!(Jinv)

    maybe_stuck = false
    max_resets = 3
    resets = 0

    for i in 1:(broyden.maxiters)
        x_old .= x
        fx_old .= fx

        p = -Jinv * fx_old

        ρ, σ₂ = T(0.9), T(0.001)

        x .= x_old .+ p
        f!(fx, x)

        if norm(fx, 2) ≤ ρ * norm(fx_old, 2) - σ₂ * norm(p, 2)^2
            α = T(1)
        else
            α = _approximate_norm_descent(f!, fx, x, p)
            x .= x_old .+ α * p
            f!(fx, x)
        end

        Δx .= x .- x_old
        Δfx .= fx .- fx_old

        maybe_stuck = all(abs.(Δx) .<= eps(T)) || all(abs.(Δfx) .<= eps(T))
        if maybe_stuck
            Jinv = _init_identity_matrix(x)
            resets += 1
            maybe_stuck = (resets ≤ max_resets) && maybe_stuck
        else
            ΔxJinv = Δx' * Jinv
            Jinv .+= ((Δx .- Jinv * Δfx) ./ (ΔxJinv * Δfx)) * ΔxJinv
        end

        maybe_stuck = false

        # Convergence Check
        norm(Δfx, 2) ≤ broyden.ϵ && return x
    end

    return x
end

# https://doi.org/10.1080/10556780008805782
# FIXME: We are dropping some robustness tests for now.
function _approximate_norm_descent(f!, fx::AbstractArray{T,N}, x::AbstractArray{T,N}, p; λ₀=T(1), β=T(0.5), σ₁=T(0.001),
                                   η=T(0.1), max_iter=50) where {T,N}
    λ₂, λ₁ = λ₀, λ₀

    f!(fx, x)
    fx_norm = norm(fx, 2)

    # TODO: Test NaN/Finite
    # f!(fx, x .- λ₂ .* p)
    # fxλp_norm = norm(fx, 2)
    # TODO: nan backtrack

    j = 0

    f!(fx, x .+ λ₂ .* p)
    converged = _test_approximate_norm_descent_convergence(f!, fx, x, fx_norm, p, σ₁, λ₂, η)

    while j < max_iter && !converged
        j += 1
        λ₁, λ₂ = λ₂, β * λ₂
        converged = _test_approximate_norm_descent_convergence(f!, fx, x, fx_norm, p, σ₁, λ₂, η)
    end

    return λ₂
end

function _test_approximate_norm_descent_convergence(f!, fx, x, fx_norm, p, σ₁, λ₂, η)
    f!(fx, x .+ λ₂ .* p)
    n1 = norm(fx, 2)

    f!(fx, x)
    n2 = norm(fx, 2)

    return n1 ≤ fx_norm - σ₁ * norm(λ₂ .* p, 2) .^ 2 + η * n2
end


text:
struct Machine{M<:Model, T<:Number}
    model::M
    preference::Preference{T}
    users::Int
    items::Int
end

model(machine::Machine) = machine.model
users(machine::Machine) = machine.users
items(machine::Machine) = machine.items
preference(machine::Machine) = machine.preference

checkuser(machine::Machine, user::Int) = (user >= 0 && user <= users(machine)) ? user : error("invalid user id ($user)")
checkitem(machine::Machine, item::Int) = (item >= 0 && item <= items(machine)) ? item : error("invalid item id ($item)")

function train!(machine::Machine, dataset::Dataset; max_epochs::Int = 1000)
    m = model(machine)
    error = Inf
    for epoch in 1:max_epochs
        update!(m, dataset)
        e = objective(m, dataset)

        if e > error
            break
        end

        error = e
    end
end

function machine(model::Model, dataset::Dataset)
    u = users(dataset)
    v = items(dataset)
    pref = preference(dataset)
    return Machine(model, pref, u, v)
end

function Base.getindex(machine::Machine{M, T}, user::Int, item::Int) where {M<:Model, T<:Number}
    checkuser(machine, user)
    checkitem(machine, item)

    m = model(machine)
    pref = preference(machine)
    value = predict(m, pref, user, item)

    if !ismissing(value)
        return PredictRating(correct(value, pref), pref)
    end

    return MissingRating{T}()
end

function Base.getindex(machine::Machine{M, T}, dataset::Dataset{T}) where {M<:Model, T<:Number}
    predicts = Vector{UserPreference{T}}(undef, length(dataset))

    for i in 1:length(dataset)
        u, v, r = dataset[i]
        predicts[i] = UserPreference{T}(u, v, machine[u, v])
    end

    return predicts
end


text:
struct QueryReport
    query::NamedTuple
    influence_curve::Vector{Float64}
    estimate::Float64
    initial_estimate::Float64
end

influencecurve(covariate, y, observed_fluct, ct_fluct, estimate) =
    covariate .* (float(y) .- observed_fluct) .+ ct_fluct .- estimate

"""
    briefreport(qr::QueryReport; tail=:both, alpha=0.05)

For a given QueryReport, provides a summary of useful statistics.

# Arguments:
    - qr: A query report, for instance extracted via `getqueryreport`
    - tail: controls weither the test is single or two sided: eg :left, :right or :both
    - alpha: level of the test
"""
function briefreport(qr::QueryReport; tail=:both, level=0.95)
    testresult = ztest(qr)
    return (query=qr.query,
            pvalue=pvalue(testresult, tail=tail),
            confint=confint(testresult, level=level, tail=tail),
            estimate=qr.estimate,
            initial_estimate=qr.initial_estimate,
            stderror=testresult.stderr,
            mean_inf_curve=mean(qr.influence_curve))
end

"""
    briefreport(mach::Machine{TMLEstimator}; tail=:both, alpha=0.05)

For a given Machine{<:TMLEstimator}, provides a summary of useful statistics for each query.

# Arguments:
    - mach: The fitted machine
    - tail: controls weither the test is single or two sided: eg :left, :right or :both
    - alpha: level of the test
"""
briefreport(mach::Machine{TMLEstimator}; tail=:both, level=0.95) =
    Tuple(briefreport(getqueryreport(mach, i), tail=tail, level=level)
            for (i, q) in enumerate(mach.model.queries))

queryreportname(i::Int) = Symbol("queryreport_$i")

getqueryreport(mach::Machine{<:TMLEstimator}, i::Int) =
    getfield(mach.report, queryreportname(i))

# Simple Test
ztest(qr::QueryReport) =
    OneSampleZTest(qr.estimate .+ qr.influence_curve)
ztest(mach::Machine{<:TMLEstimator}, arg::Int) = ztest(getqueryreport(mach, arg))
ztest(mach::Machine{<:TMLEstimator}, args::Vararg{Int}) =
    Tuple(ztest(getqueryreport(mach, arg)) for arg in args)

# Paired Test
ztest(qr₁::QueryReport, qr₂::QueryReport) =
    OneSampleZTest(qr₁.influence_curve .+ qr₁.estimate, qr₂.influence_curve .+ qr₂.estimate)
ztest(mach::Machine{<:TMLEstimator}, pair::Pair{Int, Int}) =
    ztest(getqueryreport(mach, pair[1]), getqueryreport(mach, pair[2]))

max_stars_repo_name: olivierlabayle/TMLE

max_issues_repo_name: olivierlabayle/TMLE

max_forks_repo_name: olivierlabayle/TMLE


text:
using Base: SHA1
using JSON: JSON
using HTTP: HTTP
using Pkg.GitTools: tree_hash
using Pkg.Artifacts
using SHA: sha256
using TimeZones.TZData: extract, tzdata_url, tzdata_versions

# Code loosely based upon: https://julialang.github.io/Pkg.jl/dev/artifacts/#Using-Artifacts-1
function bind_artifact_url!(artifacts_toml::String, name::String, url::String)
    # Query the `Artifacts.toml` file for the hash associated with artifact name. If no such
    # binding exists within the file then `nothing` will be returned.
    artifact_hash = Artifacts.artifact_hash(name, artifacts_toml)

    if artifact_hash === nothing
        @info "Processing new artifact: $name"

        archive = download(url)
        try
            # Compute the `sha256` of the archive.
            archive_sha = bytes2hex(open(sha256, archive))

            # Compute the `git-tree-sha1`. Usually this is computed via
            # `create_artifact` but we want to avoid actually storing this data as an
            # artifact.
            # Alternatively: `Tar.tree_hash(IOBuffer(inflate_gzip(archive)))`
            artifact_hash = mktempdir() do temp_dir
                extract(archive, temp_dir)
                SHA1(tree_hash(temp_dir))
            end

            bind_artifact!(
                artifacts_toml,
                name,
                artifact_hash,
                lazy=true,
                download_info=[(url, archive_sha)],
            )
        finally
            rm(archive)
        end
    end
end

function update_artifacts!(artifacts_toml::String)
    @info "Checking for missing tzdata versions..."
    versions = tzdata_versions()
    for version in versions
        artifact_name = "tzdata$version"
        bind_artifact_url!(artifacts_toml, artifact_name, tzdata_url(version))
    end
    latest_tzdata = last(versions)

    @info "Checking for latest Unicode CLDR release..."
    response = HTTP.get("https://api.github.com/repos/unicode-org/cldr/releases/latest")
    json = JSON.parse(String(response.body))
    latest_unicode_cldr = json["tag_name"]  # latest release

    @info "Latest Unicode CLDR release: $latest_unicode_cldr"
    url = "https://github.com/unicode-org/cldr/archive/$latest_unicode_cldr.tar.gz"
    bind_artifact_url!(artifacts_toml, "unicode-cldr-$latest_unicode_cldr", url)

    # Display the information on the latest releases which is useful for updating the
    # TimeZones package defaults.
    println("""

        Latest releases:
           tzdata:       $latest_tzdata
           unicode-cldr: $latest_unicode_cldr
        """
    )
end

# Execute with: `julia --project=dev/ dev/generate_artifacts.jl`
if abspath(PROGRAM_FILE) == @__FILE__
    artifact_toml = if length(ARGS) >= 1
        ARGS[1]
    else
        joinpath(@__DIR__, "..", "Artifacts.toml")
    end

    @time update_artifacts!(artifact_toml)
end


text:
import YaoBlocks.Optimise.to_basictypes
export simplify_blocktypes

function to_basictypes_ex(pb::PutBlock{N,C,<:ChainBlock}) where {N,C}
    chain(N, [put(N, pb.locs=>b) for b in subblocks(content(pb))])
end

function to_basictypes_ex(pb::PutBlock{N,C,<:ControlBlock}) where {N,C}
    c = content(pb)
    ControlBlock{N}(map(loc->pb.locs[loc], c.ctrl_locs), c.ctrl_config,
            content(c), map(loc->pb.locs[loc], c.locs))
end

function to_basictypes_ex(pb::PutBlock{N,C,<:PutBlock}) where {N,C}
    c = content(pb)
    PutBlock{N}(content(c), map(loc->pb.locs[loc], c.locs))
end
to_basictypes_ex(pb::AbstractBlock) = pb

simplify_blocktypes(c) = simplify(c, rules=[to_basictypes, to_basictypes_ex])


text:
module LogitTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: VBLogitOut, VBLogitIn1, VBLogitXi

#-------------
# Update rules
#-------------

@testset "VBLogitOut" begin
    @test VBLogitOut <: NaiveVariationalRule{Logit}
    @test outboundType(VBLogitOut) == Message{Bernoulli}
    @test isApplicable(VBLogitOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleVBLogitOut(nothing, ProbabilityDistribution(Univariate, GaussianMeanVariance, m=2.0, v=1.0), ProbabilityDistribution(Univariate, PointMass, m=3.0)) == Message(Univariate, Bernoulli, p=1/(1+exp(-2.0)))
end

@testset "VBLogitIn1" begin
    @test VBLogitIn1 <: NaiveVariationalRule{Logit}
    @test outboundType(VBLogitIn1) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(VBLogitIn1, [ProbabilityDistribution, Nothing, ProbabilityDistribution])

    @test ruleVBLogitIn1(ProbabilityDistribution(Univariate, Bernoulli, p=0.8), nothing, ProbabilityDistribution(Univariate, PointMass, m=3.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=0.3, w=(1/(1+exp(-3.0)) - 0.5)/3.0)
end

@testset "VBLogitXi" begin
    @test VBLogitXi <: NaiveVariationalRule{Logit}
    @test outboundType(VBLogitXi) == Message{Function}
    @test isApplicable(VBLogitXi, [ProbabilityDistribution, ProbabilityDistribution, Nothing])

    @test ruleVBLogitXi(ProbabilityDistribution(Univariate, Bernoulli, p=0.8), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=2.0, v=1.0), nothing) == Message(Univariate, Function, mode=sqrt(5.0))
end

@testset "averageEnergy" begin
    @test averageEnergy(Logit, ProbabilityDistribution(Univariate, Bernoulli, p=0.8), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=2.0, v=1.0), ProbabilityDistribution(Univariate, PointMass, m=3.0)) == (1/(1+exp(-3.0)) - 0.5)/6.0*(5.0 - 9.0) + 2.5 + log(1+exp(-3.0)) - 1.6
end

end # module


text:
using Random
traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
urdf_folder = joinpath(traj_folder, "dynamics","urdf")
urdf_kuka_orig = joinpath(urdf_folder, "kuka_iiwa.urdf")
urdf_kuka = joinpath(urdf_folder, "temp","kuka.urdf")

function write_kuka_urdf()
    kuka_mesh_dir = joinpath(TrajectoryOptimization.root_dir(),"dynamics","urdf","kuka_iiwa_mesh")
    temp_dir = joinpath(TrajectoryOptimization.root_dir(),"dynamics","urdf","temp")
    if !isdir(temp_dir)
        mkdir(temp_dir)
    end
    open(urdf_kuka_orig,"r") do f
        open(urdf_kuka, "w") do fnew
            for ln in eachline(f)
                pre = findfirst("<mesh filename=",ln)
                post = findlast("/>",ln)
                if !(pre isa Nothing) && !(post isa Nothing)
                    inds = pre[end]+2:post[1]-2
                    pathstr = ln[inds]
                    file = splitdir(pathstr)[2]
                    ln = ln[1:pre[end]+1] * joinpath(kuka_mesh_dir,file) * ln[post[1]-1:end]
                end
                println(fnew,ln)
            end
        end
    end
end

function get_kuka_ee(kuka)
    ee_body = findbody(kuka, "iiwa_link_ee")
    ee_point = Point3D(default_frame(ee_body),0.,0.,0.)
    return ee_body, ee_point
end

function get_kuka_ee_postition_fun(kuka::Mechanism,statecache=StateCache(kuka)) where {O}
    ee_body, ee_point = Dynamics.get_kuka_ee(kuka)
    world = root_frame(kuka)
    nn = num_positions(kuka)

    function ee_position(x::AbstractVector{T}) where T
        state = statecache[T]
        set_configuration!(state, x[1:nn])
        RigidBodyDynamics.transform(state, ee_point, world).v
    end
end

function calc_ee_position(kuka::Mechanism,X::Trajectory)
    ee = zero.(X)
    N = length(X)
    state = MechanismState(kuka)
    world = root_frame(kuka)
    ee_point = get_kuka_ee(kuka)[2]
    nn = num_positions(kuka)
    for k = 1:N
        set_configuration!(state, X[k][1:nn])
        ee[k] = RigidBodyDynamics.transform(state, ee_point, world).v
    end
    return ee
end

function kuka_ee_ik(kuka::Mechanism,point::Vector,ik_iterations=1000,attempts=20,tol=1e-2)
    state = MechanismState(kuka)
    world = root_frame(kuka)

    # Get end-effector
    ee_body, ee_point = get_kuka_ee(kuka)

    # Run IK
    err = Inf
    iter = 1
    while err > tol
        rand!(state)
        goal = Point3D(world,point)
        ik_res = jacobian_transpose_ik!(state,ee_body,ee_point,goal,iterations=ik_iterations)
        point_res = RigidBodyDynamics.transform(ik_res,ee_point,world).v
        err = norm(point-point_res)
        if iter > attempts
            error("IK cannot get sufficiently close to the goal")
        end
        return ik_res
    end
end

function jacobian_transpose_ik!(state::MechanismState,
                               body::RigidBody,
                               point::Point3D,
                               desired::Point3D;
                               α=0.1,
                               iterations=100)
    mechanism = state.mechanism
    world = root_frame(mechanism)

    # Compute the joint path from world to our target body
    p = path(mechanism, root_body(mechanism), body)
    # Allocate the point jacobian (we'll update this in-place later)
    Jp = point_jacobian(state, p, RigidBodyDynamics.transform(state, point, world))

    q = copy(configuration(state))

    for i in 1:iterations
        # Update the position of the point
        point_in_world = RigidBodyDynamics.transform(state, point, world)
        # Update the point's jacobian
        point_jacobian!(Jp, state, p, point_in_world)
        # Compute an update in joint coordinates using the jacobian transpose
        Δq = α * Array(Jp)' * (RigidBodyDynamics.transform(state, desired, world) - point_in_world).v
        # Apply the update
        q .= configuration(state) .+ Δq
        set_configuration!(state, q)
    end
    state
end

function hold_trajectory(n,m,N, mech::Mechanism, q)
    state = MechanismState(mech)
    nn = num_positions(state)
    set_configuration!(state, q[1:nn])
    vd = zero(state.q)
    u0 = dynamics_bias(state)

    if length(q) > m
        throw(ArgumentError("system must be fully actuated to hold an arbitrary position ($(length(q)) should be > $m)"))
    end
    U0 = zeros(m,N)
    for k = 1:N
        U0[:,k] = u0
    end
    return U0
end

# Write new urdf file with correct absolute paths
write_kuka_urdf()

kuka = Model(urdf_kuka)
end_effector_function = Dynamics.get_kuka_ee_postition_fun(parse_urdf(urdf_kuka,remove_fixed_tree_joints=false))

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
function forward(backend::GPUBackend, state::HingeLossLayerState, inputs::Vector{Blob})
  pred = inputs[1]
  label = inputs[2]

  data_type = eltype(pred)
  n = length(pred)

  x_block = div(n + CUDA.THREADS_PER_BLOCK_X-1, CUDA.THREADS_PER_BLOCK_X)

  if length(state.loss_blob) < x_block
    destroy(state.loss_blob)
    state.loss_blob = make_blob(backend, data_type, x_block)
  end

  if data_type == Float32
    kernel = backend.mocha.hinge_loss_forward_float
  elseif data_type == Float64
    kernel = backend.mocha.hinge_loss_forward_double
  else
    error("Unsupported data type $data_type")
  end

  CUDA.launch(kernel, (x_block,1), (CUDA.THREADS_PER_BLOCK_X, 1),
              (pred.ptr.p, label.ptr.p, n, state.loss_blob.ptr.p))

  losses = Array{data_type}(size(state.loss_blob)...)
  copy!(losses, state.loss_blob)
  state.loss = state.layer.weight * sum(losses[1:x_block]) / get_num(pred)

  # accumulate statistics
  state.loss_accum = (state.loss_accum*state.n_accum + state.loss*get_num(pred)) / (state.n_accum+get_num(pred))
  state.n_accum += get_num(pred)
end

function backward(backend::GPUBackend, state::HingeLossLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  const pred = inputs[1]
  const label = inputs[2]

  const data_type = eltype(pred)
  const n = length(pred)
  const neg_weight :: data_type = -state.layer.weight / get_num(pred)

  const x_block = div(n + CUDA.THREADS_PER_BLOCK_X-1, CUDA.THREADS_PER_BLOCK_X)

  if data_type == Float32
    kernel = backend.mocha.hinge_loss_backward_float
  elseif data_type == Float64
    kernel = backend.mocha.hinge_loss_backward_double
  else
    error("Unsupported data type $data_type")
  end

  const gradient1 :: CuPtr = isa(diffs[1], CuTensorBlob) ? diffs[1].ptr : CuPtr()
  const gradient2 :: CuPtr = isa(diffs[2], CuTensorBlob) ? diffs[2].ptr : CuPtr()

  CUDA.launch(kernel, (x_block,1), (CUDA.THREADS_PER_BLOCK_X,1),
              (pred.ptr.p, label.ptr.p, gradient1.p, gradient2.p, n, neg_weight))
end


text:
abstract Model # this is either an LSTM or RNN

type NNMatrix # Neural net layer's weights & gradients
    n::Int
    d::Int
    w::Matrix{Float64} # matrix of weights
    dw::Matrix{Float64} # matrix of gradients
    NNMatrix(n::Int) = new(n, 1, zeros(n,1), zeros(n,1))
    NNMatrix(n::Int, d::Int) = new(n, d, zeros(n,d), zeros(n,d))
    NNMatrix(n::Int, d::Int, w::Array, dw::Array) = new(n, d, w, dw)
end

randNNMat(n::Int, d::Int, std::FloatingPoint=1.) = NNMatrix(n, d, randn(n,d)*std, zeros(n,d))
onesNNMat(n::Int, d::Int) = NNMatrix(n, d, ones(n, d), zeros(n, d))

function softmax(m::NNMatrix)
    out = NNMatrix(m.n,m.d)
    maxval = maximum(m.w)
    out.w[:] = exp(m.w - maxval)
    out.w /= sum(out.w)
    return out
end

max_issues_repo_name: Andy-P/RecurrentNN

max_forks_repo_name: Andy-P/RecurrentNN


text:
using BinaryBuilder, Pkg

name = "Ipopt"
version = v"3.13.2"

sources = [
    ArchiveSource("https://github.com/coin-or/Ipopt/archive/releases/$(version).tar.gz",
                  "891ab9e9c7db29fc8ac5c779ccec6313301098de7bbf735ca230cd5544c49496"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd Ipopt-releases-*

# Remove misleading libtool files
rm -f ${prefix}/lib/*.la
rm -f /opt/${target}/${target}/lib*/*.la

LIBASL=(-L${libdir} -lasl)
if [[ "${target}" == *-linux-* ]]; then
  LIBASL+=(-lrt)
fi

./configure --enable-shared \
            --prefix=${prefix} \
            --with-lapack-lflags=-lopenblas \
            --with-mumps-cflags="-I$prefix/include/mumps_seq" \
            --with-mumps-lflags="-L${libdir} -ldmumps -lmpiseq -lmumps_common -lopenblas -lpord" \
            --with-asl-cflags="-I${prefix}/include" \
            --with-asl-lflags="${LIBASL[*]}" \
            --host=${target}

# parallel build fails
make
make install
"""

platforms = expand_gfortran_versions(expand_cxxstring_abis(supported_platforms()))

# The products that we will ensure are always built
products = [
    LibraryProduct("libipopt", :libipopt),
    ExecutableProduct("ipopt", :amplexe),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("ASL_jll"),
    Dependency("OpenBLAS32_jll", v"0.3.9"),  # Ipopt uses 32-bit ints
    Dependency(PackageSpec(; name="MUMPS_seq_jll", version=v"5.2.1")),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Justly
using Justly: Song, Chord, Note, Interval
using Documenter: doctest
using Test: @test, @test_throws, @testset
using AudioSchedules: AudioSchedule, duration
using Base.Meta: ParseError
using Unitful: s

cd(joinpath(pkgdir(Justly), "test")) do
    edit_song("song.justly"; test = true)
    song = read_justly("song.justly")
    @test length(song.chords[1].notes) == 4
    @test duration(AudioSchedule(song)) == 1.55s

    @testset "parsing" begin
        @test_throws ParseError("Can't parse ? as beats on line 3") read_justly(
            "bad_beats.justly",
        )
        @test_throws ParseError("Can't parse ? as chord on line 3") read_justly(
            "bad_chord.justly",
        )
        @test_throws ParseError("Can't parse ? as denominator on line 3") read_justly(
            "bad_denominator.justly",
        )
        @test_throws ParseError("Can't parse o as interval on line 3") read_justly(
            "bad_interval.justly",
        )
        @test_throws ParseError("Can't parse ? as note on line 3") read_justly(
            "bad_note.justly",
        )
        @test_throws ParseError("Can't parse ? as numerator on line 3") read_justly(
            "bad_numerator.justly",
        )
        @test_throws ParseError("Can't parse ? as octave on line 3") read_justly(
            "bad_octave.justly",
        )
        @test_throws ParseError("Can't parse ? as tempo on line 1") read_justly(
            "bad_tempo.justly",
        )
        @test_throws ParseError("Can't parse ? as initial key on line 1") read_justly(
            "bad_initial_key.justly",
        )
        @test_throws ParseError("Can't parse ? as initial key and tempo on line 1") read_justly(
            "bad_initial_key_and_tempo.justly",
        )
    end
end

if v"1.6" <= VERSION < v"1.7"
    doctest(Justly)
end

max_stars_repo_name: bramtayl/Justly

max_issues_repo_name: bramtayl/Justly

max_forks_repo_name: bramtayl/Justly


text:
# These functions (variants of sa_numerical_method)
# return an instance of a numerical method to solve
# saturation adjustment, for different combinations
# of thermodynamic variable inputs.

# @print only accepts literal strings, so we must
# branch to print which method is being used.
function print_numerical_method(
    ::Type{sat_adjust_method},
) where {sat_adjust_method}
    if sat_adjust_method <: NewtonsMethod
        @print("    Method=NewtonsMethod")
    elseif sat_adjust_method <: NewtonsMethodAD
        @print("    Method=NewtonsMethodAD")
    elseif sat_adjust_method <: SecantMethod
        @print("    Method=SecantMethod")
    elseif sat_adjust_method <: RegulaFalsiMethod
        @print("    Method=RegulaFalsiMethod")
    else
        error("Unsupported numerical method")
    end
end

#####
##### Thermodynamic variable inputs: ρ, e_int, q_tot
#####
function sa_numerical_method(
    ::Type{NM},
    param_set::APS,
    ρ::FT,
    e_int::FT,
    q_tot::FT,
    phase_type::Type{<:PhaseEquil},
) where {FT, NM <: NewtonsMethod}
    _T_min::FT = T_min(param_set)
    T_init =
        max(_T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
    return NewtonsMethod(
        T_init,
        T_ -> ∂e_int_∂T(param_set, heavisided(T_), e_int, ρ, q_tot, phase_type),
    )
end

function sa_numerical_method(
    ::Type{NM},
    param_set::APS,
    ρ::FT,
    e_int::FT,
    q_tot::FT,
    phase_type::Type{<:PhaseEquil},
) where {FT, NM <: NewtonsMethodAD}
    _T_min::FT = T_min(param_set)
    T_init =
        max(_T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
    return NewtonsMethodAD(T_init)
end

function sa_numerical_method(
    ::Type{NM},
    param_set::APS,
    ρ::FT,
    e_int::FT,
    q_tot::FT,
    phase_type::Type{<:PhaseEquil},
) where {FT, NM <: SecantMethod}
    _T_min::FT = T_min(param_set)
    q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice
    T_2 = air_temperature(param_set, e_int, q_pt)
    T_1 = max(_T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
    T_2 = bound_upper_temperature(T_1, T_2)
    return SecantMethod(T_1, T_2)
end

function sa_numerical_method(
    ::Type{NM},
    param_set::APS,
    ρ::FT,
    e_int::FT,
    q_tot::FT,
    phase_type::Type{<:PhaseEquil},
) where {FT, NM <: RegulaFalsiMethod}
    _T_min::FT = T_min(param_set)
    q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice
    T_2 = air_temperature(param_set, e_int, q_pt)
    T_1 = max(_T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
    T_2 = bound_upper_temperature(T_1, T_2)
    return RegulaFalsiMethod(T_1, T_2)
end

#####
##### Thermodynamic variable inputs: ρ, p, q_tot
#####

function sa_numerical_method_ρpq(
    ::Type{NM},
    param_set::APS,
    ρ::FT,
    p::FT,
    q_tot::FT,
    phase_type::Type{<:PhaseEquil},
) where {FT, NM <: NewtonsMethodAD}
    q_pt = PhasePartition(q_tot)
    T_init = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt)
    return NewtonsMethodAD(T_init)
end

function sa_numerical_method_ρpq(
    ::Type{NM},
    param_set::APS,
    ρ::FT,
    p::FT,
    q_tot::FT,
    phase_type::Type{<:PhaseEquil},
) where {FT, NM <: RegulaFalsiMethod}
    q_pt = PhasePartition(q_tot)
    T_1 = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt) - 5
    T_2 = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt) + 5
    return RegulaFalsiMethod(T_1, T_2)
end

#####
##### Thermodynamic variable inputs: p, e_int, q_tot
#####

function sa_numerical_method_peq(
    ::Type{NM},
    param_set::APS,
    p::FT,
    e_int::FT,
    q_tot::FT,
    phase_type::Type{<:PhaseEquil},
) where {FT, NM <: NewtonsMethodAD}
    _T_min::FT = T_min(param_set)
    T_init =
        max(_T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
    return NewtonsMethodAD(T_init)
end

function sa_numerical_method_peq(
    ::Type{NM},
    param_set::APS,
    p::FT,
    e_int::FT,
    q_tot::FT,
    phase_type::Type{<:PhaseEquil},
) where {FT, NM <: SecantMethod}
    _T_min::FT = T_min(param_set)
    q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice
    T_2 = air_temperature(param_set, e_int, q_pt)
    T_1 = max(_T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
    T_2 = bound_upper_temperature(T_1, T_2)
    return SecantMethod(T_1, T_2)
end


text:
function build_quad_diff_cache!(model)
    problem_data, index_map = get_problem_data(model.optimizer)
    (
        Q, q, G, h, A, b,
        nz, var_list,
        nineq_le, le_con_idx,
        nineq_ge, ge_con_idx,
        nineq_sv_le, le_con_sv_idx,
        nineq_sv_ge, ge_con_sv_idx,
        neq, eq_con_idx,
        neq_sv, eq_con_sv_idx,
    ) = problem_data

    z = MOI.get(model.optimizer, MOI.VariablePrimal(), var_list)

    # separate λ, ν

    _λ = -MOI.get.(model.optimizer, MOI.ConstraintDual(), le_con_idx)
    append!(
        _λ,
        MOI.get.(model.optimizer, MOI.ConstraintDual(), ge_con_idx),
    )
    append!(
        _λ,
        -MOI.get.(model.optimizer, MOI.ConstraintDual(), le_con_sv_idx),
    )
    append!(
        _λ,
        MOI.get.(model.optimizer, MOI.ConstraintDual(), ge_con_sv_idx),
    )
    λ = convert(Vector{Float64}, _λ)
    # We want to stay consistent with the variable `ν` defined in (3) of
    # Left hand side of eq. (6) in https://arxiv.org/pdf/1703.00443.pdf
    # However, in eq. (6), they put it in the lagrangian as
    # `+ ν ⋅ (Az - b)`
    # while in MOI, we put it as
    # `- ν ⋅ (Az - b)`
    # so the we should reverse the sign if we want to use the same equations
    # as in the paper.
    _ν = -MOI.get.(model.optimizer, MOI.ConstraintDual(), eq_con_idx)
    append!(
        _ν,
        -MOI.get.(model.optimizer, MOI.ConstraintDual(), eq_con_sv_idx),
    )
    ν = convert(Vector{Float64}, _ν)

    LHS = create_LHS_matrix(z, λ, Q, G, h, A)
    model.gradient_cache = QPCache(
        problem_data,
        λ,
        ν,
        z,
        LHS,
        index_map,
    )
    return nothing
end

# TODO: create test functions for the methods

# """
#     Left hand side of eqn(6) in https://arxiv.org/pdf/1703.00443.pdf
# """
# function create_LHS_matrix(z, λ, Q, G, h, A=nothing)
#     if A == nothing || size(A)[1] == 0
#         return [Q                G';
#                 Diagonal(λ) * G    Diagonal(G * z - h)]
#     else
#         @assert size(A)[2] == size(G)[2]
#         p, n = size(A)
#         m    = size(G)[1]
#         return [Q                  G'                    A';
#                 Diagonal(λ) * G    Diagonal(G * z - h)   zeros(m, p);
#                 A                  zeros(p, m)           zeros(p, p)]
#     end
# end

"""
    create_LHS_matrix(z, λ, Q, G, h, A=nothing)

Inverse matrix specified on RHS of eqn(7) in https://arxiv.org/pdf/1703.00443.pdf

Helper method while calling `_backward_quad`
"""
function create_LHS_matrix(z, λ, Q, G, h, A=nothing)::AbstractMatrix{Float64}
    if (A === nothing || size(A)[1] == 0) && (G === nothing || size(G)[1] == 0)
        return Q
    elseif A === nothing || size(A)[1] == 0
        return [Q         G' * Diagonal(λ);
                G         Diagonal(G * z - h)]
    elseif G === nothing || size(G)[1] == 0
        p, n = size(A)
        return [Q         A';
                A         spzeros(p, p)]
    else
        p, n  = size(A)
        m, n2 = size(G)
        if n != n2
            throw(DimensionError("Sizes of $A and $G do not match"))
        end
        return [Q         G' * Diagonal(λ)       A';
                G         Diagonal(G * z - h)    spzeros(m, p);
                A         spzeros(p, m)          spzeros(p, p)]
    end
end
# TODO: this is the transpose, check back for usage

# """
#     Right hand side of eqn(6) in https://arxiv.org/pdf/1703.00443.pdf
# """
# function create_RHS_matrix(z, dQ, dq, λ, dG, dh, ν=nothing, dA=nothing, db=nothing)
#     if dA == nothing || size(dA)[1] == 0
#         return -[dQ * z + dq + dG' * λ      ;
#                  Diagonal(λ) * (dG * z - dh)]
#     else
#         return -[dQ * z + dq + dG' * λ + dA' * ν;
#                  Diagonal(λ) * (dG * z - dh)    ;
#                  dA * z - db                    ]
#     end
# end

const _QP_SET_TYPES = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    # MOI.Interval{Float64},
}

const _QP_FUNCTION_TYPES = Union{
    MOI.VariableIndex,
    MOI.ScalarAffineFunction{Float64},
}

_qp_supported(::Type{F}, ::Type{S}) where {F <: _QP_FUNCTION_TYPES, S <: _QP_SET_TYPES} = true
_qp_supported(::Type{F}, ::Type{S}) where {F, S} = false
function _qp_supported(model::MOI.AbstractOptimizer)
    return all(FS -> _qp_supported(FS...), MOI.get(model, MOI.ListOfConstraintTypesPresent()))
end

"""
    get_problem_data(model::MOI.AbstractOptimizer)

Return problem parameters as matrices along with other program info such as number of constraints, variables, etc
"""
function get_problem_data(model::MOI.AbstractOptimizer)
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        if !_qp_supported(F, S)
            throw(MOI.UnsupportedConstraint{F,S}("DiffOpt does not support this constraint type for its Quadratic Programming differentiation. Maybe try the Conic Programming differentiation ? For this, do `MOI.set(model, DiffOpt.ProgramClass(), DiffOpt.CONIC)`."))
        end
    end
    var_list = MOI.get(model, MOI.ListOfVariableIndices())
    nz = length(var_list)

    index_map = MOIU.IndexMap()
    for (i,vi) in enumerate(var_list)
        index_map[vi] = VI(i)
    end

    # handle inequality constraints
    le_con_idx = MOI.get(
                        model,
                        MOI.ListOfConstraintIndices{
                            MOI.ScalarAffineFunction{Float64},
                            MOI.LessThan{Float64},
                        }())
    ge_con_idx = MOI.get(
                        model,
                        MOI.ListOfConstraintIndices{
                            MOI.ScalarAffineFunction{Float64},
                            MOI.GreaterThan{Float64},
                        }())
    nineq_le = length(le_con_idx)
    nineq_ge = length(ge_con_idx)
    le_con_sv_idx = MOI.get(
                        model,
                        MOI.ListOfConstraintIndices{
                            MOI.VariableIndex,
                            MOI.LessThan{Float64},
                        }())
    ge_con_sv_idx = MOI.get(
                        model,
                        MOI.ListOfConstraintIndices{
                            MOI.VariableIndex,
                            MOI.GreaterThan{Float64},
                        }())
    nineq_sv_le = length(le_con_sv_idx)
    nineq_sv_ge = length(ge_con_sv_idx)

    G = spzeros(nineq_le + nineq_ge + nineq_sv_le + nineq_sv_ge, nz)
    h = spzeros(nineq_le + nineq_ge + nineq_sv_le + nineq_sv_ge)

    ineq_cont = 0
    eq_cont = 0

    for i in 1:nineq_le
        con = le_con_idx[i]

        func = MOI.get(model, MOI.ConstraintFunction(), con)
        set = MOI.get(model, MOI.ConstraintSet(), con)

        for (j, var_idx) in enumerate(var_list)
            for term in func.terms
                if term.variable == var_idx
                    G[i,j] = MOI.coefficient(term)
                end
            end
        end
        h[i] = set.upper - func.constant

        ineq_cont += 1
        index_map[con] =
            CI{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}(ineq_cont)
    end
    for i in 1:nineq_ge
        # note: ax >= b needs to be converted in Gx <= h form
        con = ge_con_idx[i]

        func = MOI.get(model, MOI.ConstraintFunction(), con)
        set = MOI.get(model, MOI.ConstraintSet(), con)

        for (j, var_idx) in enumerate(var_list)
            for term in func.terms
                if term.variable == var_idx
                    G[i+nineq_le,j] = -MOI.coefficient(term)
                end
            end
        end
        h[i+nineq_le] = func.constant - set.lower
        ineq_cont += 1
        index_map[con] =
            CI{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}(ineq_cont)
    end
    for i in eachindex(le_con_sv_idx)
        con = le_con_sv_idx[i]
        func = MOI.get(model, MOI.ConstraintFunction(), con)
        set = MOI.get(model, MOI.ConstraintSet(), con)
        vidx = findfirst(v -> v == func, var_list)
        G[i+nineq_le+nineq_ge,vidx] = 1
        h[i+nineq_le+nineq_ge] = MOI.constant(set)
        ineq_cont += 1
        index_map[con] =
            CI{MOI.VariableIndex, MOI.LessThan{Float64}}(ineq_cont)
    end
    for i in eachindex(ge_con_sv_idx)
        # note: x >= b needs to be converted in Gx <= h form
        con = ge_con_sv_idx[i]
        func = MOI.get(model, MOI.ConstraintFunction(), con)
        set = MOI.get(model, MOI.ConstraintSet(), con)
        vidx = findfirst(isequal(func), var_list)
        G[i+nineq_le+nineq_ge+nineq_sv_le,vidx] = -1
        h[i+nineq_le+nineq_ge+nineq_sv_le] = -MOI.constant(set)
        ineq_cont += 1
        index_map[con] =
            CI{MOI.VariableIndex, MOI.GreaterThan{Float64}}(ineq_cont)
    end

    # handle equality constraints
    eq_con_idx = MOI.get(
                        model,
                        MOI.ListOfConstraintIndices{
                            MOI.ScalarAffineFunction{Float64},
                            MOI.EqualTo{Float64}
                        }())
    neq = length(eq_con_idx)

    eq_con_sv_idx = MOI.get(
        model,
        MOI.ListOfConstraintIndices{
            MOI.VariableIndex,
            MOI.EqualTo{Float64}
        }())
    neq_sv = length(eq_con_sv_idx)

    A = spzeros(neq + neq_sv, nz)
    b = spzeros(neq + neq_sv)

    for i in 1:neq
        con = eq_con_idx[i]

        func = MOI.get(model, MOI.ConstraintFunction(), con)
        set = MOI.get(model, MOI.ConstraintSet(), con)

        for x in func.terms
            # never nothing, variable is present
            vidx = findfirst(isequal(x.variable), var_list)
            A[i, vidx] = x.coefficient
        end
        b[i] = set.value - func.constant

        eq_cont += 1
        index_map[con] =
            CI{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}(eq_cont)
    end
    for i in 1:neq_sv
        con = eq_con_sv_idx[i]
        func = MOI.get(model, MOI.ConstraintFunction(), con)
        set = MOI.get(model, MOI.ConstraintSet(), con)
        vidx = findfirst(isequal(func), var_list)
        A[i+neq,vidx] = 1
        b[i+neq] = set.value
        eq_cont += 1
        index_map[con] =
            CI{MOI.VariableIndex, MOI.EqualTo{Float64}}(eq_cont)
    end

    # handle objective
    # works both for any objective function convertible to a ScalarQuadraticFunction.
    # So in particular VariableIndex, ScalarAffineFunction and ScalarQuadraticFunction should work.
    objective_function = MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    sparse_array_obj = sparse_array_representation(objective_function, nz, index_map)

    return (
        sparse_array_obj.quadratic_terms, sparse_array_obj.affine_terms,
        G, h, A, b,
        nz, var_list,
        nineq_le, le_con_idx,
        nineq_ge, ge_con_idx,
        nineq_sv_le, le_con_sv_idx,
        nineq_sv_ge, ge_con_sv_idx,
        neq, eq_con_idx,
        neq_sv, eq_con_sv_idx,
    ), index_map
end


text:
export AbstractInput,
    StartTargetInput, TargetSystemInput, HomotopyInput, ParameterSystemInput

const Inputs = Union{AbstractSystem,MPPolys,Composition,ModelKit.System}
const MPPolyInputs = Union{MPPolys,Composition}

const input_supported_keywords = [
    :parameters,
    :generic_parameters,
    :start_parameters,
    :target_parameters,
    :target_gamma,
    :start_gamma,
    :p₁,
    :p₀,
    :γ₁,
    :γ₀,
    :variable_ordering,
    # deprecated
    :startparameters,
    :targetparameters,
    :targetgamma,
    :startgamma,
]

"""
    AbstractInput

An abstract type to model different input types.
"""
abstract type AbstractInput end

"""
    StartTargetInput(start::MPPolyInputs, target::MPPolyInputs, startsolutions)

Construct a `StartTargetInput` out of two systems of polynomials.
"""
struct StartTargetInput{P1<:Inputs,P2<:Inputs} <: AbstractInput
    start::P1
    target::P2

    function StartTargetInput{P1,P2}(start::P1, target::P2) where {P1<:Inputs,P2<:Inputs}
        if length(start) ≠ length(target)
            error("Cannot construct `StartTargetInput` since the lengths of `start` and `target` don't match.")
        end
        new(start, target)
    end
end
function StartTargetInput(start::P1, target::P2) where {P1<:Inputs,P2<:Inputs}
    StartTargetInput{P1,P2}(start, target)
end

"""
    HomotopyInput(H::Homotopy, startsolutions)

Construct a `HomotopyInput` for a homotopy `H` with given `startsolutions`.
"""
struct HomotopyInput{Hom<:AbstractHomotopy} <: AbstractInput
    H::Hom
end

"""
    TargetSystemInput(system::Inputs)

Construct a `TargetSystemInput`. This indicates that the system `system`
is the target system and a start system should be assembled.
"""
struct TargetSystemInput{S<:Inputs} <: AbstractInput
    system::S
end

"""
    ParameterSystemInput(F, parameters, p₁, p₀, startsolutions, γ₁=nothing, γ₂=nothing)

Construct a `ParameterSystemInput`.
"""
struct ParameterSystemInput{S<:Inputs} <: AbstractInput
    system::S
    parameters::Union{Nothing,Vector{<:MP.AbstractVariable}}
    p₁::AbstractVector
    p₀::AbstractVector
    γ₁
    γ₀
end

"""
    input_startsolutions(F::MPPolyInputs)
    input_startsolutions(F::AbstractSystem)
    input_startsolutions(G::MPPolyInputs, F::MPPolyInputs, startsolutions)
    input_startsolutions(F::MPPolyInputs, parameters, startsolutions; kwargs...)
    input_startsolutions(H::AbstractHomotopy, startsolutions)

Construct an `AbstractInput` and pass through startsolutions if provided.
Returns a named tuple `(input=..., startsolutions=...)`.
"""
function input_startsolutions(
    F::Union{MPPolyInputs,ModelKit.System};
    variable_ordering = nothing,
    # parameter homotopy indication,
    generic_parameters = nothing,
    start_parameters = generic_parameters,
    p₁ = start_parameters,
    target_parameters = generic_parameters,
    p₀ = target_parameters,
    kwargs...,
)
    # if p₀ or p₁ is provided this is actually the input constructor
    # for a parameter homotopy, but no startsolutions are provided
    if !isnothing(p₁) || !isnothing(p₀)
        return parameter_homotopy(
            F,
            nothing;
            variable_ordering = variable_ordering,
            p₁ = p₁,
            p₀ = p₀,
            kwargs...,
        )
    end

    if isa(F, MPPolyInputs)
        if variable_ordering !== nothing && nvariables(F) != length(variable_ordering)
            throw(ArgumentError("Number of assigned variables is too small."))
        end
        remove_zeros!(F)
        if has_constant_polynomial(F)
            throw(ArgumentError("System contains a non-zero constant polynomial."))
        end
    end

    (input = TargetSystemInput(F), startsolutions = nothing)
end

function input_startsolutions(F::AbstractSystem; variable_ordering = nothing)
    (input = TargetSystemInput(F), startsolutions = nothing)
end

# create ModelKit.System and use this case
function input_startsolutions(
    F::Vector{<:ModelKit.Expression},
    args...;
    variable_ordering::Union{Nothing,Vector{ModelKit.Variable}} = nothing,
    parameters::Union{Nothing,Vector{ModelKit.Variable}} = nothing,
    kwargs...,
)
    if variable_ordering === nothing
        throw(ArgumentError("`variable_ordering = ...` needs to be passed as a keyword argument to indicate the order of the variables."))
    end
    if isnothing(parameters)
        system = ModelKit.System(F, variable_ordering)
    else
        system = ModelKit.System(F, variable_ordering, parameters)
    end
    input_startsolutions(system, args...; kwargs...)
end

function input_startsolutions(
    G::MPPolyInputs,
    F::MPPolyInputs,
    startsolutions = nothing;
    variable_ordering = nothing,
)
    if length(G) ≠ length(F)
        throw(ArgumentError("Start and target system don't have the same length"))
    end
    if variable_ordering !== nothing && (
        nvariables(F) != length(variable_ordering) ||
        nvariables(G) != length(variable_ordering)
    )
        throw(ArgumentError("Number of assigned variables is too small."))
    end

    check_zero_dimensional(F)
    if startsolutions === nothing
        startsolutions = [randn(ComplexF64, nvariables(F))]
    elseif isa(startsolutions, AbstractVector{<:Number})
        startsolutions = [startsolutions]
    end
    (input = StartTargetInput(G, F), startsolutions = startsolutions)
end

# need
input_startsolutions(F::Inputs, starts; kwargs...) =
    parameter_homotopy(F, starts; kwargs...)

function parameter_homotopy(
    F::Inputs,
    startsolutions;
    parameters = F isa MPPolyInputs ?
                 throw(ArgumentError("Argument `parameters = ...` needs to be provided.")) :
                 nothing,
    variable_ordering = nothing,
    generic_parameters = nothing,
    start_parameters = generic_parameters,
    p₁ = start_parameters,
    target_parameters = generic_parameters,
    p₀ = target_parameters,
    start_gamma = nothing,
    γ₁ = start_gamma,
    target_gamma = nothing,
    γ₀ = target_gamma,
)
    if isnothing(p₁)
        throw(ArgumentError(
            "You need to pass `generic_parameters = ...`, " *
            "`start_parameters = ...` or `p₁ = ...` as a keyword argument",
        ))
    end
    if isnothing(p₀)
        throw(ArgumentError(
            "`target_parameters = ...` or `p₀ = ...` need " *
            "to be passed as a keyword argument.",
        ))
    end
    if length(p₁) != length(p₀)
        throw(ArgumentError("Number of parameters doesn't match!"))
    end
    if F isa MPPolyInputs && length(parameters) != length(p₀)
        throw(ArgumentError("Number of parameters doesn't match!"))
    end
    if F isa MPPolyInputs && !isnothing(variable_ordering) &&
       nvariables(F, parameters) != length(variable_ordering)
        throw(ArgumentError("Number of assigned variables is too small."))
    end

    if F isa MPPolyInputs && variable_ordering !== nothing &&
       nvariables(F, parameters) != length(variable_ordering)
        throw(ArgumentError("Number of assigned variables is too small."))
    end

    if isnothing(startsolutions)
        startsolutions = [randn(ComplexF64, nvariables(F, parameters))]
    elseif isa(startsolutions, AbstractVector{<:Number})
        startsolutions = [startsolutions]
    end

    (
        input = ParameterSystemInput(F, parameters, p₁, p₀, γ₁, γ₀),
        startsolutions = startsolutions,
    )
end

function input_startsolutions(
    H::AbstractHomotopy,
    startsolutions;
    variable_ordering = nothing,
)
    if isa(startsolutions, AbstractVector{<:Number})
        startsolutions = [startsolutions]
    end
    (input = HomotopyInput(H), startsolutions = startsolutions)
end


text:
# Structs and Functions
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#chap16
#include("ch15module.jl")

using Printf

# Time
"""
Represents the time of day.

hour:
minute:
second:
"""
mutable struct MyTime
    hour::Int64
    minute::Int64
    second::Int64
    # Inner constructor
    function MyTime(hour::Int64=0, minute::Int64=0,
        second::Int64=0)
        new(hour, minute, second)
    end
end

# Exercise 16-1
function printtime(mytime)
    @printf("%02d:%02d:%02d\n", mytime.hour, mytime.minute, mytime.second)
end

# Exercise 16-2
function isafter(t2, t1)
    tosecond(t2) < tosecond(t1)
end
function tosecond(mytime)
    mytime.hour*60*60 + mytime.minute*60 + mytime.second
end

# Pure Function
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_pure_functions

# A pure function does not modify its parameter
function addtimeDEAD(t1, t2)
    second = t1.second + t2.second
    minute = t1.minute + t2.minute
    hour = t1.hour + t2.hour

    if second >= 60
        second -= 60
        minute += 1
    end
    if minute >= 60
        minute -= 60
        hour += 1
    end

    MyTime(hour, minute, second)
end

# Improved version
function addtime(t1, t2)
    totalsecond = timetosecond(t1) + timetosecond(t2)
    secondtotime(totalsecond)
end

# Modifiers
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#modifiers
# """
# Represents the time of day.
#
# hour:
# minute:
# second:
# """
# mutable struct MMyTime
#     hour::Int64
#     minute::Int64
#     second::Int64
# end

# Modifier version
function increment!DEAD(time, seconds)
    time.second += seconds

    if time.second >= 60
        time.second -= 60
        time.minute += 1
    end
    if time.minute >= 60
        time.minute -= 60
        time.hour += 1
    end
    # Exercise 16-3
    if time.second >= 60
        increment!(time, 0)
    end
end

# Exercise 16-5
function increment!(time, seconds)
    totalsecond = timetosecond(time) + seconds
    newtime = secondtotime(totalsecond)

    time.second = newtime.second
    time.minute = newtime.minute
    time.hour = newtime.hour
end

# Exercise 16-4: Pure version
"""
    increment

Add seconds to the time.

```
time:
addsecond:
```
"""
function increment(time::MyTime, addsecond::Int64)
    second = time.second + addsecond
    minute = time.minute
    hour = time.hour

    if second >= 60
        second -= 60
        minute += 1
    end
    if minute >= 60
        minute -= 60
        hour += 1
    end

    rettime = MyTime(hour, minute, second)

    if rettime.second >= 60
    # Recursive to move second to minute to hour
        rettime = increment(rettime, 0)
    end

    return rettime
end

# Prototyping Versus Planning
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#prototyping_versus_planning
function timetosecond(time)
    seconds = time.hour*3600 + time.minute*60 + time.second
end

function secondtotime(totalsecond)
    (minute, second) = divrem(totalsecond, 60)
    (hour, minute) = divrem(minute, 60)
    return MyTime(hour, minute, second)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2b2f0b1a-3265-11ec-1237-95caba1fa9ab
md"""

# Managing external packages for your Genie app

Genie fully takes advantage of Julia's excellent package manager, `Pkg` -- while allowing Genie developers to use any 3rd party package available in Julia's ecosystem. This is achieved by taking a common sense approach: Genie apps are just plain Julia projects. This can be easily confirmed by looking at the familiar `Project.toml` and `Manifest.toml` files in the root of your Genie app.

In order to add extra packages to your Genie app, thus, we need to use Julia's `Pkg` mode:

1. start a Genie REPL with your app: `$ bin/repl`. This will automatically load the package environment of the app.
2. switch to `Pkg` mode: `julia> ]`
3. add the package you want, for example `OhMyREPL`: `(MyGenieApp) pkg> add OhMyREPL`

That's all! Now you can use the packages at the Genie REPL or anywhere in your files via `using` or `import`.

Use the same approach to update all the packages, via: `pkg> up` to apply all available updates, or `pkg> up OhMyREPL` to update a single package.

Finally, you can use `pkg> remove OhMyREPL` to uninstall a dependency.

"""

# ╔═╡ Cell order:
# ╠═2b2f0b1a-3265-11ec-1237-95caba1fa9ab


text:
using PyPlot
using LaTeXStrings

#data = readdlm("data",'\t')
#data_13dpg = readdlm("combined_data_combo_13dpg")
#data = vcat(data,data_13dpg)
#data = readdlm("data_corrected",'\t')
#combined_data_combo_new = readdlm("combined_data_combo_new")
#combined_data_combo_13dpg = readdlm("13dpg")
#data = vcat(data,combined_data_combo_new,combined_data_combo_13dpg)

data = readdlm("FINAL",'\t')
data = readdlm("combined_data_SA_238ens",'\t')

labels = data[:,1]
error = data[:,2]
FVA_norm = data[:,3]

species_indices = Dict()
for i in 1:length(labels)
	species_indices[labels[i]] = i
end

error /= error[1]
FVA_norm /= FVA_norm[1]

m = 8
a = .3

no_glc = Array{Int64}(readdlm("no_glc"))
no_g6p = Array{Int64}(readdlm("no_g6p"))

#figure(figsize=(22.5,11.5))
figure(figsize=(12,8))

#semilogx(error[3:end],FVA_norm[3:end],"go",markersize=m,markeredgecolor="none",alpha=a)
#for i in 3:length(error)
#	if in(i,no_glc) && !in(i,no_g6p)
#		semilogx(error[i],FVA_norm[i],color="blue","o",markersize=m,markeredgecolor="none",alpha=a)
#	elseif in(i,no_g6p) && !in(i,no_glc)
#		semilogx(error[i],FVA_norm[i],"ro",markersize=m,markeredgecolor="none",alpha=a)
#	elseif in(i,no_glc) && in(i,no_g6p)
#		semilogx(error[i],FVA_norm[i],color="k","o",markersize=m,markeredgecolor="none",alpha=a)
#	else
#		semilogx(error[i],FVA_norm[i],"go",markersize=m,markeredgecolor="none",alpha=a)
#	end
#end
for i in 3:length(error)
	if in(i,no_glc)
		semilogx(error[i],FVA_norm[i],color="orange","o",markersize=m,markeredgecolor="none",alpha=a)
	else
		semilogx(error[i],FVA_norm[i],"go",markersize=m,markeredgecolor="none",alpha=a)
	end
end

#semilogx(error[2],FVA_norm[2],"ro",markersize=m,markeredgecolor="none")
semilogx(error[1],FVA_norm[1],"k*",markersize=m*2,markeredgecolor="none")
xlabel("error / base case error")
ylabel("flux uncertainty / base case flux uncertainty")
#xlabel("error")
#ylabel("flux uncertainty")

ax = collect(axis())
#ax[3] = 15000
#ax[4] = 25000
axis(ax)

#for i in 1:length(labels)
#	key = labels[i]
#	if in(key,["–GLC"; "–PYR"; "–ADP"; "–ATP"; "+CIT"; "–TYR"; "–UMP"; "+R5P"; "–MAL"])
#		annotate(key,xy=(error[species_indices[key]]+.005,FVA_norm[species_indices[key]]))
#	elseif in(key,["–AC"])
#		annotate(key,xy=(error[species_indices[key]],FVA_norm[species_indices[key]]-.003))
#	elseif in(key,["+AKG"])
#		annotate(key,xy=(error[species_indices[key]]-.082,FVA_norm[species_indices[key]]-.0015))
#	elseif in(key,["–PHE"])
#		annotate(key,xy=(error[species_indices[key]]-.035,FVA_norm[species_indices[key]]-.005))
#	elseif in(key,["+2PG"])
#		annotate(key,xy=(error[species_indices[key]]-.04,FVA_norm[species_indices[key]]-.005))
#	elseif in(key,["+F6P"])
#		annotate(key,xy=(error[species_indices[key]]+.007,FVA_norm[species_indices[key]]-.002))
#	elseif in(key,["+G6P"])
#		annotate(key,xy=(error[species_indices[key]]-.04,FVA_norm[species_indices[key]]-.005))
#	end
#end

savefig("New.pdf")


text:
using GeneralizedSampling
using Base.Test

println("Testing NFFT...")

#=
Test if the (N)FFT part of a Freq2Wave type behaves like a DFT matrix by
removing the parts related to wavelets.

For integer samples [0:N-1] the DFT matrix is Hermitian.
=#

# ------------------------------------------------------------
# 1D

begin
	N = 16.0
	samples = collect( 0:N-1 )

	J = Int(log2(N))
	T = Freq2Wave(samples, "haar", J)
	for idx in eachindex(T.internal)
		T.internal[idx] = 1.0 + 0.0*im
	end

	F = collect(T)
	B = F'*F / N

	@test_approx_eq B eye(B)
end

# ------------------------------------------------------------
# 2D

begin
	M = 4
	samples = grid(M)

	J = Int(log2(M))
	T = Freq2Wave(samples, "haar", J)
	for idx in eachindex(T.internal)
		T.internal[idx] = 1.0 + 0.0*im
	end

	F = collect(T)
	N = size(T,2)
	B = F'*F / N

	@test_approx_eq B eye(B)
end


text:
using GeomDAE.Fields: Field, Q, P, V, F

"Special Ordinary Differential Equation"
immutable SODE{T} <: Equation{T}
    d::Int
    n::Int

    q₀::Array{T, 2}
    p₀::Array{T, 2}
    t₀::T

    q::Function
    p::Function
    v::Function
    f::Function

    fields::Tuple{Vararg{Field}}
    functs::Tuple{Vararg{Field}}

    function SODE(d, n, q₀, p₀, t₀, q, p, v, f)
        @assert d == size(q₀,1) == size(p₀,1)
        @assert n == size(q₀,2) == size(p₀,2)
        @assert T == eltype(q₀) == eltype(p₀)

        if ndims(q₀) == 1
            q₀ = reshape(q₀, d, n)
        end

        if ndims(p₀) == 1
            p₀ = reshape(p₀, d, n)
        end

        fields = []
        functs = []

        for (field, symb) in ((q, Q()),
                              (p, P()),
                              (v, V()),
                              (f, F()))
            if field == do_nothing
                push!(fields, symb)
            else
                push!(functs, symb)
            end
        end

        @assert length(fields) == length(functs) == 2

        new(d, n, q₀, p₀, t₀, q, p, v, f, tuple(fields...), tuple(functs...))
    end
end

function SODE{T}(q₀::DenseArray{T}, p₀::DenseArray{T}, t₀::Real=0;
                 q::Function=do_nothing, p::Function=do_nothing,
                 v::Function=do_nothing, f::Function=do_nothing)
    @assert size(q₀) == size(p₀)
    SODE{T}(size(q₀, 1), size(q₀, 2), q₀, p₀, t₀, q, p, v, f)
end

Base.hash(ode::SODE, h::UInt) = hash(ode.d, hash(ode.n,
                                hash(ode.t₀, hash(ode.q₀, hash(ode.p₀,
                                hash(ode.q, hash(ode.p, hash(ode.v, hash(ode.f,
                                hash(ode.fields, hash(ode.functs, h)))))))))))

Base.:(==){T1, T2}(ode1::SODE{T1}, ode2::SODE{T2}) = (ode1.d == ode2.d
                                                   && ode1.n == ode2.n
                                                   && ode1.t₀ == ode2.t₀
                                                   && ode1.q₀ == ode2.q₀
                                                   && ode1.p₀ == ode2.p₀
                                                   && ode1.q == ode2.q
                                                   && ode1.p == ode2.p
                                                   && ode1.v == ode2.v
                                                   && ode1.f == ode2.f
                                                   && ode1.fields == ode2.fields
                                                   && ode1.functs == ode2.functs)


text:
procs = addprocs(12)

using OPCSPs
using MCTS

N = 1000
par = true

problems = [gen_informative(noise_2norm=100.0,
                                 n_nodes=10,
                                 n_connected=4,
                                 rbar=10.0,
                                 rng=MersenneTwister(i))
            for i in 1:N]

iss = initial_states(problems, rng_offset=100)

naive = evaluate_performance(problems, iss, GurobiExactSolver(multithreaded=false), rng_offset=1000, parallel=par)
@show mean(naive)

cheating = evaluate_performance(problems, iss, Cheater(), rng_offset=1000, parallel=par)
@show mean(cheating)

upper = evaluate_performance(problems, iss, OnlyUnobservableUncertainty(), rng_offset=1000, parallel=par)
@show mean(upper)

mean_feedback = evaluate_performance(problems, iss, FeedbackSolver(GurobiExactSolver(multithreaded=false)), rng_offset=1000, parallel=par)
@show mean(mean_feedback)

s = GurobiExactSolver(multithreaded=false)

srng = MersenneTwister(1947)
solver = AgUCTSolver(
    aggregator=VoronoiOPCSPAg(50.0),
    rollout_solver=FeedbackSolver(s),
    exploration_constant=100.0,
    n_iterations=10000,
    rng=srng
)

@time ag = evaluate_performance(problems, iss, solver, rng_offset=1000, parallel=par)
@show mean(ag)

rmprocs(procs)


text:
abstract type Shooting end
struct shooting end
abstract type Collocation end
struct collocation end

function optimization_template(model, u0, xf, x0, ::Shooting;
    quadrature = rk4(),
    nt=10,
    ts=range(0.0, stop=1.0, length=nt),
    lu = zeros(length(ts) * model.nu),  # lower bounds on u
    uu = ones(length(ts) * model.nu) .* 10000 * pi/30,  # upper bounds on u
    objective = (x) -> x[end],
    solver = SNOW.IPOPT()
)
    xs = zeros(length(ts) * model.nx)
    function objective!(g, u)
        integrate_dyanmics!(xs, model, x0, u, ts, quadrature)
        # x0 constraint
        g[1:model.nx] .= xs[1:model.nx] .- x0
        # xf constraint
        g[model.nx+1:2*model.nx] .= xs[end-model.nx+1:end] .- xf

        return objective(xs)
    end

    ng = 5*model.nx  # number of constraints
    lg = -Inf*ones(ng)  # lower bounds on g
    lg[1:2*model.nx] .= 0.0 # equality constraints on start and end states
    ug = zeros(ng)  # upper bounds on g
    options = SNOW.Options(solver=solver)  # choosing IPOPT solver

    return ts, xs, objective!, u0, ng, lu, uu, lg, ug, options
end

function optimization_template(model, u0, xf, ::Shooting;
    quadrature = rk4(),
    nt=10,
    ts=range(0.0, stop=1.0, length=nt),
    lu = zeros(length(ts) * model.nu),  # lower bounds on u
    uu = ones(length(ts) * model.nu) .* 10000 * pi/30,  # upper bounds on u
    objective = (x) -> x[end],
    solver = SNOW.IPOPT()
)
    x0 = trim(model)
end

# uopt, fopt, info = SNOW.minimize(objective!, u0, ng, lu, uu, lg, ug, options)
# println("ustar = ", uopt)
# println("fstar = ", fopt)
# println("info = ", info)

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
seed!(42)
L = 3
D = 10
N = 20
nSamples = 10
b = 5
x = rand(N, D)
y = rand(N)
@testset "Inference" begin
    i = AnalyticVI()
    @test length(i) == 1
end


text:
mutable struct ReferenceGraph
    graph
    adjacency_graph::WeightedGraph
    get_vertex_property::Function
    get_edge_property::Function
    edge_property_limits::Dict{Symbol,Any} # Any to allow lists, ranges, tuples (first and second must be defined on it)
    vertex_property_limits::Dict{Symbol,Any}
end

function ReferenceGraph(
    g::AbstractGraph;
    get_vertex_property=(args...)->nothing,
    get_edge_property=(args...)->nothing,
    edge_property_limits=Dict{Symbol, Any}(),
    vertex_property_limits=Dict{Symbol, Any}()
)
    adjacency_graph = is_directed(g) ? SimpleDiGraph(g) : SimpleGraph(g)
    return ReferenceGraph(g, deepcopy(adjacency_graph), get_vertex_property, get_edge_property, vertex_property_limits, edge_property_limits)
end

function ReferenceGraph(
    g,
    adjacency_matrix;
    get_vertex_property=(args...)->nothing,
    get_edge_property=(args...)->nothing,
    edge_property_limits=Dict{Symbol, Any}(),
    vertex_property_limits=Dict{Symbol, Any}()
)
    adjacency_graph = graph_from_matrix(adjacency_matrix)
    return ReferenceGraph(g, adjacency_graph, get_vertex_property, get_edge_property, vertex_property_limits, edge_property_limits)
end


text:
# include -
include("Include.jl")

# This script requires Flux.jl (we manually include here, its takes a while to load)
using Flux
using Flux: @epochs

# how many training data sets do we have?
number_of_synthetic_patients = 500

trajectory_range = range(1,stop=2001,step=100) |> collect
N = length(trajectory_range)

# initialize storage for the training data -
training_data = Vector{Tuple{Vector{Float32},Vector{Float32}}}()
for i ∈ 1:number_of_synthetic_patients

    # load the output -
    output_df = CSV.read(joinpath(_PATH_TO_TMP,"OUT-Synthetic-1k-P$(i).csv"), DataFrame)
    output_vector = convert.(Float32, output_df[!,:simulated])

    # load the FIIa -
    simulation_df = CSV.read(joinpath(_PATH_TO_TMP,"SIM-Synthetic-1k-P$(i).csv"), DataFrame)
    FIIa_vector = convert.(Float32, simulation_df[trajectory_range,:FIIa])

    # make a training tuple -
    training_tuple = (
        output_vector,FIIa_vector
    );

    # insert -
    push!(training_data, training_tuple)
end

# build the model -
FIIa_trajectory_model = Chain(Dense(5, 10, σ), Dense(10, N));

# setup a loss function -
loss(x, y) = Flux.Losses.mae(FIIa_trajectory_model(x), y; agg = mean)

# pointer to params -
ps = Flux.params(FIIa_trajectory_model)

# # use old school gradient descent -
opt = Momentum(0.1, 0.95)

@epochs 1000 Flux.train!(loss, ps, training_data, opt)

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
#!/bin/bash
export JULIA_NUM_THREADS=4

#good_orders = 55 .+ collect(1:14);
 #good_orders = 55 .+ collect(17:24)
 #good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41))
# julia -e 'idx_day_to_use=1; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=2; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=3; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=4; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=5; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=6; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=7; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=8; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=9; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=10; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=11; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
# julia -e 'idx_day_to_use=12; good_orders = 55 .+ collect(1:14); include("examples/calc_ccfs_subset_orders.jl")';
#julia -e 'idx_day_to_use=13; include("examples/calc_ccfs.jl")';
#julia -e 'idx_day_to_use=14; include("examples/calc_ccfs.jl")';

#good_orders = 55 .+ collect(17:24);
 julia -e 'idx_day_to_use=1; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=1; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';

 julia -e 'idx_day_to_use=2; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=2; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';

 julia -e 'idx_day_to_use=3; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=3; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=4; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=4; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';

 julia -e 'idx_day_to_use=5; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=5; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';

 julia -e 'idx_day_to_use=6; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=6; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';

  julia -e 'idx_day_to_use=7; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=8; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=9; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=10; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=11; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=12; good_orders = 55 .+ collect(17:24); include("examples/calc_ccfs_subset_orders.jl")';
#julia -e 'idx_day_to_use=13; include("examples/calc_ccfs.jl")';
#julia -e 'idx_day_to_use=14; include("examples/calc_ccfs.jl")';

#good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));
 julia -e 'idx_day_to_use=7; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=8; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=9; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=10; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=11; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';
 julia -e 'idx_day_to_use=12; good_orders = 55 .+ vcat(collect(27:29), collect(37:38),collect(40:41));  include("examples/calc_ccfs_subset_orders.jl")';
#julia -e 'idx_day_to_use=13; include("examples/calc_ccfs.jl")';
#julia -e 'idx_day_to_use=14; include("examples/calc_ccfs.jl")';


text:
VERSION < v"0.5" && warn("This example uses dot broadcast syntax "*
                         "and may not work on Julia versions lower"*
                         " than v0.5!")

using BenchmarkTools

# Define a parent BenchmarkGroup to contain our suite
const suite = BenchmarkGroup()

# Add some child groups to our benchmark suite.
suite["string"] = BenchmarkGroup(["unicode"])
suite["trig"] = BenchmarkGroup(["math", "triangles"])
suite["dot"] = BenchmarkGroup(["broadcast", "elementwise"])

# This string will be the same every time because we're seeding the RNG
teststr = join(rand(MersenneTwister(1), 'a':'d', 10^4))

# Add some benchmarks to the "utf8" group
suite["string"]["replace"] = @benchmarkable replace($teststr, "a", "b") seconds=Float64(π)
suite["string"]["join"] = @benchmarkable join($teststr, $teststr) samples=42

# Add some benchmarks to the "trig"/"dot" group
for f in (sin, cos, tan)
    for x in (0.0, pi)
        suite["trig"][string(f), x] = @benchmarkable $(f)($x)
        suite["dot"][string(f), x] = @benchmarkable $(f).([$x, $x, $x])
    end
end

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `suite` every time the file is included.
paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
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 IoK8sApiPolicyV1beta1IDRange <: SwaggerModel
    max::Any # spec type: Union{ Nothing, Int64 } # spec name: max
    min::Any # spec type: Union{ Nothing, Int64 } # spec name: min

    function IoK8sApiPolicyV1beta1IDRange(;max=nothing, min=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1IDRange, Symbol("max"), max)
        setfield!(o, Symbol("max"), max)
        validate_property(IoK8sApiPolicyV1beta1IDRange, Symbol("min"), min)
        setfield!(o, Symbol("min"), min)
        o
    end
end # type IoK8sApiPolicyV1beta1IDRange

const _property_map_IoK8sApiPolicyV1beta1IDRange = Dict{Symbol,Symbol}(Symbol("max")=>Symbol("max"), Symbol("min")=>Symbol("min"))
const _property_types_IoK8sApiPolicyV1beta1IDRange = Dict{Symbol,String}(Symbol("max")=>"Int64", Symbol("min")=>"Int64")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1IDRange }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1IDRange))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1IDRange }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiPolicyV1beta1IDRange[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1IDRange }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1IDRange[property_name]

function check_required(o::IoK8sApiPolicyV1beta1IDRange)
    (getproperty(o, Symbol("max")) === nothing) && (return false)
    (getproperty(o, Symbol("min")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiPolicyV1beta1IDRange }, name::Symbol, val)
end


text:
using Interpolations,Plots

t = 0:.1:1
x = sin.(2π*t)
y = cos.(2π*t)
A = hcat(x,y)

itp = scale(interpolate(A, (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), t, 1:2)

tfine = 0:.01:1
xs, ys = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]

scatter(x, y, label="knots")
plot!(xs, ys, label="spline")

max_stars_repo_name: MarsBighead/mustang

max_issues_repo_name: MarsBighead/mustang

max_forks_repo_name: MarsBighead/mustang