text:
######################################################################
# Pia Link - Pia's standard transport protocol.
# -----
# Licensed under Apache 2.0
module LinkProtocol
using ..Futures
using ..Lockables
using ..Pia
using ..Packets
using ..PacketFactories
using ..Applications
using Base.Threads
using ExtraFun
using Sockets

# ===== UDPSocket Wrapper =====
export Socket
struct Socket
    raw::UDPSocket
    endpoint::Endpoint
end
function Socket(endpoint::Endpoint)
    socket = UDPSocket()
    bind(socket, endpoint.host, endpoint.port)
    Socket(socket, endpoint)
end
Socket(host::IPAddr, port::Integer) = Socket(Endpoint(host, port))

Base.write(socket::Socket, remote::Endpoint, data) = send(socket.raw, remote.host, remote.port, data)
Base.read(socket::Socket) = recvfrom(socket.raw)

Base.close(socket::Socket) = close(socket.raw)
Base.isopen(socket::Socket) = isopen(socket.raw)

# ===== Protocol =====
export LinkPeer
mutable struct LinkPeer <: Peer
    endpoint::Endpoint
    deadline::Float64
end

export PiaLink
mutable struct PiaLink <: TransportProtocol
    num_desired_connections::UInt16
    num_max_connections::UInt16
    peer_lifetime::Float64

    socket::Socket
    wakeup_socket::Socket
    bootstrap::Set{Endpoint}
    connections::Lockable{Dict{Endpoint, Peer}}
    expectance::Lockable{Dict{Endpoint, Future{Vector{UInt8}}}}
    terminate::Atomic{Bool}
end
function PiaLink(
        host::IPAddr = IPv6(0),
        port::Integer = DEFAULT_PORT;
        num_desired_connections::Integer = Pia.DEFAULT_DESIRED_CONNECTIONS,
        num_max_connections::Integer = Pia.DEFAULT_CONNECTION_HARDCAP,
        peer_lifetime::Real = Pia.DEFAULT_PEER_LIFETIME,
        wakeup_port::Integer = randport(50000)
    )
    PiaLink(
        num_desired_connections,
        num_max_connections,
        peer_lifetime,

        Socket(host, port),
        Socket(loopback(host), wakeup_port),
        Set(),
        Lockable(Dict{Endpoint, Peer}()),
        Lockable(Dict{Endpoint, Future{Vector{UInt8}}}()),
        Atomic{Bool}(false)
    )
end

Pia.endpoint(link::PiaLink) = link.socket.endpoint
Pia.endpoint(peer::LinkPeer) = peer.endpoint

Pia.connect(link::PiaLink, app::AbstractApplication, host::IPAddr, port::Integer = Pia.DEFAULT_PORT) = Pia.connect(link, app, Endpoint(host, port))
function Pia.connect(link::PiaLink, app::AbstractApplication, endpoint::Endpoint; timeout::Real = Pia.DEFAULT_TIMEOUT)::Peer
    if haskey(link.connections[], endpoint)
        lock(link.connections) do
            if haskey(link.connections[], endpoint)
                return link.connections[][endpoint]
            end
        end
    end

    @with ch = PeerChannel(link, endpoint) begin
        handshake_id = rand(UInt8)
        write(ch, packetbytes(Packets.HandshakeStart, appname(app), appversion(app), handshake_id))

        (pktcls, pktid), response_id = packetcontent(read(ch, timeout=timeout))
        if pktcls !== Packets.ConnectivityClass || pktid !== Packets.HandshakeResponse || response_id != handshake_id+0x1
            throw(ConnectionError("invalid handshake response"))
        end

        write(ch, packetbytes(Packets.HandshakeConfirm))
    end

    peer = LinkPeer(endpoint, time() + Pia.DEFAULT_PEER_LIFETIME)
    lock(link.connections) do
        link.connections[][endpoint] = peer
    end
    return peer
end
function Pia.connect(link::PiaLink, app::AbstractApplication)::Vector{Peer}
    if isempty(link.bootstrap)
        throw(StateError("no known bootstrap peers"))
    end

    peers = Peer[]
    errors = []
    for remote ∈ link.bootstrap
        try
            push!(peers, Pia.connect(link, app, remote))
        catch ex
            push!(errors, ex) # TODO: Capture stack trace
        end
    end
    if isempty(peers)
        throw(ConnectionError("failed to connect to bootstrap peers", errors))
    end
    return peers
end

Pia.isconnected(link::PiaLink) = !isempty(link.connections[])

Pia.bootstrap(link::PiaLink, ip::IPAddr, port::Integer = Pia.DEFAULT_PORT) = Pia.bootstrap(link, Endpoint(ip, port))
function Pia.bootstrap(link::PiaLink, endpoint::Endpoint)
    push!(link.bootstrap, endpoint)
    link
end

Base.write(link::PiaLink, ::AbstractApplication, endpoint::Endpoint, data) = write(link, endpoint, packetbytes(Packets.AppData, data))
Base.write(link::PiaLink, endpoint::Endpoint, data) = write(link.socket, endpoint, data)

function Base.read(link::PiaLink, ::AbstractApplication, endpoint::Endpoint; kwargs...)
    (clsid, pktid), bytes = packetcontent(read(link, endpoint; kwargs...))
    @assert clsid === Packets.DataClass
    @assert pktid === Packets.AppData
    return bytes
end
function Base.read(link::PiaLink, endpoint::Endpoint; timeout::Real = Pia.DEFAULT_TIMEOUT)
    future = Future{Vector{UInt8}}()
    lock(link.expectance) do
        link.expectance[][endpoint] = future
    end
    return wait(future, timeout=timeout)
end

function Base.close(link::PiaLink)
    if !atomic_cas!(link.terminate, false, true)
        write(link.wakeup_socket, Endpoint(loopback(link.socket.endpoint.host), link.socket.endpoint.port), packetbytes(Packets.Wakeup))

        lock(link.connections) do
            bytes = packetbytes(Packets.Goodbye)
            for peer ∈ values(link.connections[])
                write(link.socket, peer.endpoint, bytes)
            end
        end

        close(link.wakeup_socket)
        close(link.socket)
    end
end

function Base.show(io::IO, link::PiaLink)
    endpoint    = link.socket.endpoint
    num_conns   = length(link.connections[])
    num_max     = link.num_max_connections
    num_desired = link.num_desired_connections

    print(io, "Pia Link @ $endpoint($num_conns/$num_max peers ($num_desired desired))")
end

# ===== PeerInterface Implementation =====
struct LinkPeerInterface <: PeerInterface
    link::PiaLink
end

Pia.peers(link::PiaLink) = LinkPeerInterface(link)

function Base.getindex(iface::LinkPeerInterface, endpoint::Endpoint)
    lock(iface.link.connections) do
        return iface.link.connections[][endpoint]
    end
end

Base.haskey(iface::LinkPeerInterface, endpoint::Endpoint) = haskey(iface.link.connections[], endpoint)
Base.lock(cb, iface::LinkPeerInterface) = lock(cb, iface.link.connections)
Base.lock(iface::LinkPeerInterface) = lock(iface.link.connections)
Base.unlock(iface::LinkPeerInterface) = unlock(iface.link.connections)

function refresh(iface::LinkPeerInterface, endpoint::Endpoint)
    link = iface.link
    lock(link.connections) do
        if haskey(link.connections[], endpoint)
            link.connections[][endpoint].deadline = time()+link.peer_lifetime
        end
    end
    return iface
end

# ===== Pia Link Loop =====
# --- Connection Multiplexing Solution ---
function Pia.loop(link::PiaLink, app::AbstractApplication)
    while !link.terminate[]
        remote, packet = @await yield() recvfrom(link.socket.raw)
        if isself(link, remote) continue end

        try
            _, pktid = peekpacket(packet)
            Applications.process_packet(Ident{pktid}(), app, PeerChannel(link, remote), packet)
        catch ex
            println(stderr, ex)
        end
    end
end
isself(link::PiaLink, other::Endpoint) = isself(link.socket.endpoint, other) || isself(link.wakeup_socket.endpoint, other)
isself(me::Endpoint, other::Endpoint) = me.host == other.host && me.port == other.port

# ===== Packet Integration =====
# --- Expectancy Anti-DOS Solution ---
function Applications.process_packet(::Ident, ::AbstractApplication, link::PiaLink, endpoint::Endpoint, bytes)
    lock(link.expectance) do
        if haskey(link.expectance[], endpoint)
            link.expectance[][endpoint][] = bytes
        end
    end
    nothing
end

# --- Handshake Start Packet ---
function Applications.process_packet(::Ident{Packets.HandshakeStart}, app::AbstractApplication, link::PiaLink, remote::Endpoint, bytes)
    # Validate packet
    remote_appname, remote_appversion, handshake_id = packetcontent(bytes, without_id=true)
    if remote_appname != appname(app)
        return LogConnectivity("invalid app name $remote_appname from peer $(ch.endpoint)")
    end
    if remote_appversion != appversion(app)
        return LogConnectivity("invalid app version $remote_appversion from peer $(ch.endpoint)")
    end

    # Send HandshakeResponse
    write(link, remote, packetbytes(Packets.HandshakeResponse, handshake_id+0x1))

    # Anticipate HandshakeConfirm (async w/ timeout)
    @async begin
        (rcls, rpkt), _ = packetcontent(read(link, remote; timeout=Pia.DEFAULT_TIMEOUT))
        if rcls === Packets.ConnectivityClass && rpkt === Packets.HandshakeConfirm
            lock(link.connections) do
                link.connections[][remote] = LinkPeer(remote, time()+link.peer_lifetime)
            end
        else
            LogConnectivity("invalid handshake confirm from peer $(ch.endpoint)")
        end
    end
end

# --- Goodbye Packet ---
function Applications.process_packet(::Ident{Packets.Goodbye}, ::AbstractApplication, link::PiaLink, remote::Endpoint, _)
    lock(link.connections) do
        delete!(link.connections[], remote)
    end
end

end # module LinkProtocol


text:
# This file is part of JuliaFMI.
# Licensed under MIT: https://github.com/AnHeuermann/JuliaFMI/blob/master/LICENSE.txt

# This file contains tools to verify results from a csv file.

import DataFrames
import Trajectories
import CSV

"""
    csvFilesEqual(csvFile1Path::String, csvFile1Path::String, epsilon::Real)

Checks if trajectories of variables in two CSV files are equal.
"""
function csvFilesEqual(csvFile1Path::String, csvFile2Path::String,
    epsilon::Real)

    # Load csv files
    csvData1 = CSV.File(csvFile1Path) |> DataFrame
    csvData2 = CSV.File(csvFile2Path) |> DataFrame

    # Compare trajectories for all variables
    allVars = setdiff(names(csvData1), [:time])

    # grep only those variables that are in both csv files
    checkvars = intersect(allVars, names(csvData2))

    return csvCompareVars(csvData1, csvData2, checkvars, epsilon)
end

"""
    csvFilesEqual(csvFile1Path::String, csvFile2Path::String, checkVars::Array{String,1}, epsilon::Real)

Checks if trajectories of variables `checkVars` in two CSV files are equal.
"""
function csvFilesEqual(csvFile1Path::String, csvFile2Path::String,
    checkVars::Array{String,1}, epsilon::Real)

    # Load csv files
    csvData1 = CSV.File(csvFile1Path) |> DataFrame
    csvData2 = CSV.File(csvFile2Path) |> DataFrame

    # Compare trajectories for specified variables
    return csvCompareVars(csvData1, csvData2, checkVars, epsilon)
end

"""
    csvFilesEqual(csvFile1Path::String, csvFile2Path::String, checkVars::Array{String,1})

Checks if trajectories of variables `checkVars` in two CSV files are equal for nonspecific epsilon.
"""
function csvFilesEqual(csvFile1Path::String, csvFile2Path::String,
    checkVars::Array{String,1})

    # Load csv files
    csvData1 = CSV.File(csvFile1Path) |> DataFrame
    csvData2 = CSV.File(csvFile2Path) |> DataFrame
    epsilon = 0.118

    # Compare trajectories for specified variables
    return csvCompareVars(csvData1, csvData2, checkVars, epsilon)
end

"""
    csvFilesEqual(csvFile1Path::String, csvFile2Path::String)

Checks if trajectories in two CSV files are equal for nonspecific epsilon.
"""
function csvFilesEqual(csvFile1Path::String, csvFile2Path::String)

    # Load csv files
    csvData1 = CSV.File(csvFile1Path) |> DataFrame
    csvData2 = CSV.File(csvFile2Path) |> DataFrame
    epsilon = 0.118

    # Compare trajectories for all variables
    allVars = setdiff(names(csvData1), [:time])

    # grep only those variables that are in both csv files
    checkvars = intersect(allVars, names(csvData2))

    # Compare trajectories for specified variables
    return csvCompareVars(csvData1, csvData2, checkvars, epsilon)
end

"""
    csvCompareVars(csvData1:::DataFrames.DataFrame, csvData2:::DataFrames.DataFrame, checkVars::Array{Symbol,1}, epsilon::Real)

Helper function to check if trajectories of variables `checkVars` in two
DataFrames are equal.
"""
function csvCompareVars(csvData1::DataFrames.DataFrame,
    csvData2::DataFrames.DataFrame, checkVars::Array{Symbol,1},
    epsilon::Real)

    isEqual = true

    if  isempty(checkVars)
        @error "Cannot compare Files. They not have any same Variables."
        return false
    end

    # rename variables like der(x) to der_x_
    newNames = [replaceNames!(names(csvData1)), replaceNames!(names(csvData2))]
    for (i, t) in enumerate(newNames[1])
        rename!(csvData1, names(csvData1)[i] => newNames[1][i])
    end
    for (i, t) in enumerate(newNames[2])
        rename!(csvData2, names(csvData2)[i] => newNames[2][i])
    end
    checkVars = replaceNames!(checkVars)

    # Check if comparable
    if unique(names(csvData1)) != names(csvData1)
        @error "Got duplicate names in csv file. Can't compare."
        return false
    end

    # Compare trajectories for each specified variables
    for varName in checkVars
        infoStr = "Checking Trajectories for Variable: $varName  \t"
        trajectory1 = Trajectories.trajectory(csvData1.time, csvData1[!,Symbol(varName)])
        trajectory2 = Trajectories.trajectory(csvData2.time, csvData2[!,Symbol(varName)])

        try
            if trajectoriesEqual(trajectory1::Trajectories.Trajectory, trajectory2::Trajectories.Trajectory, epsilon)
                infoStr *= "true"
            else
                infoStr *= "false"
                isEqual = false
            end

        catch
            infoStr *= "false"
            isEqual = false
        end
        @info infoStr
    end

    @info "Check complete."
    return isEqual
end

"""
    function trajectoriesEqual(trajectory1::Trajectory, trajectory2::Trajectory, epsilon)

Compares if two trajectories are equal by linear interpolation and comparing to
definded error ϵ.
"""
function trajectoriesEqual(trajectory1::Trajectories.Trajectory, trajectory2::Trajectories.Trajectory,
    epsilon::Real)

    time1 = trajectory1.t
    values1 = trajectory1.x
    time2 = trajectory2.t
    values2 = trajectory2.x

    intersectionTime = vcat(max(time1[1], time2[1]), time1[max(time1[1], time2[1]) .< time1 .< min(time1[end], time2[end])], min(time1[end], time2[end]))
    if intersectionTime[1] >= intersectionTime[end]
        error("Can't compare trajectories. Time intervalls not intersecting")
    end

    #find Events and check if equal in both files
    events_1 = []
    for (i, t) in enumerate(time1)
        if i != 1
            if time1[i - 1] == time1[i]
                append!(events_1, time1[i])
            end
        end
    end
    events_1 = unique(events_1)
    events_2 = []
    for (i, t) in enumerate(time2)
        if i != 1
            if time2[i - 1] == time2[i]
                append!(events_2, time2[i])
            end
        end
    end
    events_2 = unique(events_2)

    if length(events_1) != length(events_2)
        error("Files have not the same number of events")
    elseif length(events_1) != 0
        for i in 1:lastindex(events_1)
            if abs(events_1[i] - events_2[i]) > epsilon
                error("Events have not the same timestamps.")
            end
        end
    end
    # filter values
    intersection_time_values1 = values1[findall(x->x == intersectionTime[1], time1)[1]:findall(x->x == intersectionTime[end], time1)[1]]
    for (i, t) in enumerate(intersectionTime)
        if intersectionTime[i] in events_1
                continue
        end
        if abs(Trajectories.interpolate(Trajectories.Linear(), trajectory2, t) - intersection_time_values1[i]) > epsilon
            return false
        end
    end

    return true
end

"""
Helper function to replace brackets in DataFrames names.
"""
function replaceNames!(data::Array{Symbol,1})

    newNames = Array{String}(undef, length(data))
    for (i, name) in enumerate(data)
        newNames[i] = String(data[i])
        newNames[i] = replace(newNames[i], "(" => "_")
        newNames[i] = replace(newNames[i], ")" => "_")
    end

    return Symbol.(newNames)
end

ext: jl

max_stars_repo_name: AnHeuermann/JuliaFMU

max_issues_repo_name: AnHeuermann/JuliaFMU

max_forks_repo_name: jwscook/JuliaFMI


text:
include("../src/vm_basis.jl")

function test_dim_Sk()
   print("dimensions.Sk...")

   for k in [-2, 0, 2, 3, 4, 6, 8, 10, 14]
      @test dim_Sk(k) == 0
   end

   for k in [12, 16, 18, 20, 22, 26]
      @test dim_Sk(k) == 1
   end

   @test dim_Sk(24) == 2

   println("PASS")
end

function test_dim_Mk()
   print("dimensions.Mk...")

   for k in [-2, 2, 3]
      @test dim_Mk(k) == 0
   end

   for k in [0, 4, 6, 8, 10, 14]
      @test dim_Mk(k) == 1
   end

   for k in [12, 16, 18, 20, 22, 26]
      @test dim_Mk(k) == 2
   end

   @test dim_Mk(24) == 3

   println("PASS")
end

function test_dimensions()
   test_dim_Sk()
   test_dim_Mk()

   println("")
end


text:
x = [0 1 0 0 1;
     0 0 1 0 0;
     0 0 0 1 0;
     0 0 0 0 1;
     0 0 0 0 0]

#=
@inline function _istriu(A::AbstractMatrix, k)
    m, n = size(A)
    for j in 1:min(n, m + k - 1)
        all(iszero, view(A, max(1, j - k + 1):m, j)) || return false
    end
    return true
end

    function AdjacencyMap(n::Integer, cptr::Vector, rval::Vector, m::M, d::StaticBool) where {M}
        T = promote_type(eltype(cptr), eltype(rval))
        _check_type(n, T)
        # String interpolation is a performance bottleneck when it's part of the same function,
        # ensure we only do it once committed to the error.
        throwstart(ckp) = throw(ArgumentError("$ckp == colptr[1] != 1"))
        throwmonotonic(ckp, ck, k) = throw(ArgumentError("$ckp == colptr[$(k-1)] > colptr[$k] == $ck"))
        _check_length("colptr", cptr, n+1, String) # don't check upper bound
        ckp = T(1)
        ckp == cptr[1] || throwstart(ckp)
        @inbounds for k = 2:n+1
            ck = cptr[k]
            ckp <= ck || throwmonotonic(ckp, ck, k)
            ckp = ck
        end
        _check_length("rowval", rval, ckp - 1, T)
        if !(M <: NoMetadata)
            # we allow empty nzval
            _check_length("metadata", m, 0, T)
        end

        # silently shorten rval and nzval to usable index positions.
        maxlen = abs(widemul(m, n))
        isbitstype(T) && (maxlen = min(maxlen, typemax(Ti) - 1))
        length(rval) > maxlen && resize!(rval, maxlen)
        length(nzval) > maxlen && resize!(nzval, maxlen)
        _goodbuffers(n, cptr, rval, m)
        new{T,M,d === static(true)}(DynamicCount(n), cptr, rval, m)
    end
=#


text:
    import CircularList: List, Node
    import Base: getindex, setindex!, size, IndexStyle, copy
    import Base: deleteat!, insert!, findfirst
    function getnode(CL::List, i::Int)
        node = head(CL)
        forward_steps = mod(i-1, length(CL))
        if forward_steps == 0
            return node
        end
        backward_steps = abs(mod(i-1, -length(CL)))
        if forward_steps <= backward_steps
            for i=1:forward_steps
                node = node.next
            end
        else
            for i=1:backward_steps
                node = node.prev
            end
        end
        node
    end
    function getindex(CL::List, i::Int)
        node = getnode(CL, i)
        node.data
    end
    function getindex(CL::List{T}, r::UnitRange) where T
        n = length(r)
        nodes = Vector{Node{T}}( undef, n )
        nodes[1] = getnode( CL, first(r) )
        for i=2:n
            nodes[i] = nodes[i-1].next
        end
        sublist = List{T}(nodes, nodes[1], n, n, n)
    end
    function getindex(CL::List, I::Int...)
        if length(I) == 0
            CL.current
        else
            getindex(CL, I[1])
        end
    end
    function setindex!(CL::List{T}, v, i::Int) where T
        node = getnode(CL, i)
        node.data = v
    end
    function setindex!(CL::List{T}, v, r::UnitRange{S}) where {T, S}
        node = getnode(CL, first(r))
        for i=1:length(r)
            node.data = v[ mod1(i, length(v) ) ]
            node = node.next
        end
    end
    function insert!(CL::List{T}, i::Integer, item) where {T}
        current = CL.current
        jump!( CL, getnode(CL,i-1) )
        insert!(CL, item)
        CL.current = current
        CL
    end
    function insert!(CL::List{T}, i::Integer, items::AbstractVector{S}) where {T,S}
        current = CL.current
        jump!( CL, getnode(CL,i-1) )
        for item in items
            insert!(CL, convert(T,item))
        end
        CL.current = current
        CL
    end
    "Insert items::List{T} such that the items.node[1] is at index i"
    function insert!(CL::List{T}, i::Integer, items::List{T}) where T
        insert!(CL, getnode( CL, i - 1), items)
    end
    "Insert items::List{T} after node"
    function insert!(CL::List{T}, node::Node{T}, items::List{T}) where T
        # No allocations needed!
        #current = CL.current
        #jump!( CL, getnode(CL,i-1) )
        if node in items.nodes
            error("Cannot insert items into themselves")
        end

        # Remove snippet from where it came from first
        source_prev = items.nodes[1].prev
        source_next = items.nodes[end].next

        # Link old ends
        source_prev.next = source_next
        source_next.prev = source_prev

        # Identify nodes at target
        # We must do this at the beginning
        target_prev = node
        target_next = target_prev.next

        # Connect beginning of insert
        target_prev.next = items.nodes[1]
        items.nodes[1].prev = target_prev

        # Connect end of insert
        items.nodes[end].next = target_next
        target_next.prev = items.nodes[end]

        # Reset head node to where it was
        #CL.current = current
        CL
    end
    function copy(CL::List{T}) where {T}
        List{T}(CL.nodes, CL.current, CL.length, CL.last, CL.capacity)
    end
    function deleteat!(CL::List, i::Integer)
        current = CL.current
        jump!( CL, getnode( CL, i ) )
        delete!(CL)
        CL.current = current
        CL
    end
    function deleteat!(CL::List, r::UnitRange)
        current = CL.current
        jump!( CL, getnode( CL, last(r) ) )
        for i=1:length(r)
            delete!(CL)
        end
        # CL.current = current
        CL
    end
    function findfirst( predicate::Function, CL::List )
        current = CL.current
        node = current
        counter = 1
        while !predicate(node.data)
            node = node.next
            counter += 1
        end
        counter
    end
    function findfirstnode( predicate::Function, CL::List )
        current = CL.current
        node = current
        while !predicate(node.data)
            node = node.next
        end
        node
    end

    # Node
    getindex(node::Node) = node.data
    function getindex(node::Node, i::Int)
        i -= 1
        if i >= 0
            for i=1:i
                node = node.next
            end
        else
            for i=1:abs(i)
                node = node.prev
            end
        end
        node
    end
    function getindex(node::Node, r::UnitRange)
        node = node[ first(r) ]
        collect( Iterators.take(node, length(r)) )
    end
    import Base: iterate, eltype, IteratorSize
    function iterate(node::Node, state = node)
        (state, state.next)
    end
    eltype(::Type{Node{T}}) where T = Node{T}
    IteratorSize(::Type{Node{T}}) where T = Base.IsInfinite()


text:
using Pkg
Pkg.activate(".")

using PowerDynamics
using OrdinaryDiffEq
using CSV
using Plots
plotly()

##

dir = @__DIR__

## set per unit MV
const base_power = 1E6 # 1MW
const base_voltage = 20E3 # 20kV
const base_current = base_power / base_voltage # 50A
const base_admittance = base_power / base_voltage^2 # 0.0025Ω^-1
const ω = 2 * π * 50.0 # 314.1593rad/s
# per unit HV
const base_voltage_HV = 110E3 # 110kV
const base_admittance_HV = base_power / base_voltage_HV^2 # 8.264462809917356e-5
const base_current_HV = base_power / base_voltage_HV

## import special functions

# ADN control layer
include("$dir/control.jl")
# custom types
include("$dir/DGUnit.jl")
include("$dir/OLTC.jl")
# load actual data
include("$dir/cigre_static.jl")
# load pre-defined test cases
include("$dir/mwe.jl")
# failure model
include("$dir/short_circuit.jl")

## setup system

# define events
t_step = 1.
t_fault = 3. # onset of node short circuit
t_duration = 0.15 # duration of node short circuit

Ron = 2.5 # Ω
Znsc = complex(Ron, 0.0) * base_admittance #  fixed by Emocosy
nsc_node = 2 # sc on low-voltage side

#S_total = complex(-0.49998009701576795, -0.49758893704171214)
S_total = complex(0.5000200975012079,0.5024116164039469)
P_ref(t) = t > t_step ? 0.9 * real(S_total) : real(S_total)
Q_ref(t) = t > t_step ? 0.9 * imag(S_total) : imag(S_total)
#P_ref(t) = real(S_total)
#Q_ref(t) = imag(S_total)

## load data from powerfactory for comparison

powerfactory =
    CSV.File("$dir/FRT_Vac_id_iq.csv", header = true, normalizenames = true)

## setup system

DG = DGUnit(;
    K_pll = 1632.993, #Hz/pu
    K_PT1 = 1.0,
    T_PT1 = 0.001,
    K_FRT = 2.0,
    I_max = 1.0, # pu
    Vref = 1.0, # pu
    Vdead = 0.1, #pu
    S_pq = -complex(0.5, 0.5),
    Y_n = 0.0,
)

# Minimalbeispiel -> S=complex(0.5, 0.5) positive !!
PQ = PQAlgebraic(
    P = -0.5,
    Q = -0.5,
    Y_n = 0.,
)

## find operating point
pg = testbed(PQ; P_ref = t->P_ref(0), Q_ref = t->Q_ref(0))
ic_guess = initial_guess(pg)
ic_pf = pf_sol(PowerGrid(pg), ic_guess, nothing)

ic_pf[1,:s] # -> S_total

cpg = testbed(DG; P_ref = t->P_ref(0), Q_ref = t->Q_ref(0))

# set all initial voltages to slack value
op = find_operationpoint(cpg, initial_guess(cpg, ic_pf[:, :u]))
check_operationpoint(cpg, op)

# ode = rhs(cpg)
# du = similar(op.vec)
# p = cpg.controller(op.vec, nothing, 0.)
# ode(du, op.vec, nothing, 0.)
# du

@show op[2, :iqp]

@show op[2, :P_g]
@show op[2, :Q_g]

@show op[2, :v] # 1.00241 pu
@show op[2, :φ] .|> rad2deg # 0.14 deg
@show op[1, :s] # -499.98 kW, -497.61 kvar
@show op[2, :s] # 500.0 kW, 500.0 kvar
# verketteter Strom / sqrt(3)
@show op[1, :iabs] * base_current_HV * 1e-3 / sqrt(3) # 0.004 kA
@show op[2, :iabs] * base_current * 1e-3 / sqrt(3) # 0.020 kA

## check dynamics

# perturb = Perturbation(2, :θ, Dec(π))
# u0 = perturb(op)
#
# prob = ODEProblem(ode, u0.vec, (0.0, 0.1))
#
# _sol = solve(prob, Rodas4())
# sol = PowerGridSolution(_sol, cpg)

## node short circuit

cpg = testbed(DG; P_ref = P_ref, Q_ref = Q_ref)
nsc = NodeShortCircuitAlt(;
    node_number = nsc_node,
    Y = 50/Znsc,
    tspan_fault = (t_fault, t_fault + t_duration),
)
@time sol = simulate(nsc, cpg, op.vec, (0., 4.))

## plots

# TODO: bugs in solution object, e.g. (sol, 2, :v), (sol, 2, :φ)
# (sol, 2, :iabs) triggers call to rhs and hangs
# --> compare with indexing sol

plot(sol, :, :v, label = false)#, yscale=:log10)
plot(sol, 2:2, :iqp, label = false)
plot!(2.5 .+ powerfactory.b_tnow_in_s, powerfactory.s_Vac_ist, label = "PF", ms = 2, ls=:dash)
hline!(op[:, :v], c=:gray, label = "op")
#vline!([t_fault, t_fault + t_duration], label = "fault", c = :black, ls=:dot)

Δ = [
    cpg.controller(u, nothing, t) |> first
    for (t, u) in zip(sol.dqsol.t, sol.dqsol.u)
]

plot(sol.dqsol.t, first.(Δ))
plot!(sol.dqsol.t, last.(Δ))

plot(sol, 2, :iqp)

plot(sol, :, :φ, label = false)
plot!(sol, 2, :θ, c = :black, label = "PLL")
hline!(op[:, :φ], c = :gray, label = "op")
#vline!([t_fault, t_fault + t_duration], label = "fault", c = :black)
#ylims!(-0.01, 0.01)

plot(sol, 2, :P_g)
plot!(sol, 2, :Q_g)
hline!(op[2:2, :P_g], label = "op")
hline!(op[2:2, :Q_g], label = "op")
vline!([t_fault, t_fault + t_duration], label = "fault", c = :black, ls=:dot)

plot(sol, 2, :P_int)
plot!(sol, 2, :Q_int)
hline!(op[2:2, :P_int], label = "op")
hline!(op[2:2, :Q_int], label = "op")
vline!([t_fault, t_fault + t_duration], label = "fault", c = :black, ls=:dot)

inc = (last(sol.dqsol.prob.tspan) - first(sol.dqsol.prob.tspan)) / 100.0
t_array = range(
    first(sol.dqsol.prob.tspan);
    stop = last(sol.dqsol.prob.tspan),
    step = inc,
)

# currents
node = nsc_node
# current set point
I_r = [
    complex(sol(t, node, :P_g), -sol(t, node, :Q_g)) / sol(t, node, :v)
    for t in t_array
]
plot(t_array, abs.(I_r), label = "I_r_$node")

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
module APSInstructions

using HDF5, Compat

export APSInstr, wf_entry, flow_call, flow_return, flow_repeat, load_repeat, sync, wait_trig, flow_loadcmp, flow_cmp, flow_goto
# write_file, write_hdf5_file

import Base.show

type APSInstr
	data::Vector{UInt16}
	name::Symbol
	target::Symbol
end

function __init__()
	global PRINTSHORT = false
	global const emptySymbol = symbol("")
	global const OPCODES = @compat Dict(
		0x0 => :WFM,
		0x1 => :MARKER,
		0x2 => :WAIT,
		0x3 => :LOAD_REPEAT,
		0x4 => :REPEAT,
		0x5 => :CMP,
		0x6 => :GOTO,
		0x7 => :CALL,
		0x8 => :RETURN,
		0x9 => :SYNC,
		0xA => :PREFETCH,
		0xB => :LOAD_CMP)
end

APSInstr(data::Vector{UInt16}) = APSInstr(data, emptySymbol, emptySymbol)
APSInstr(data::Vector{UInt16}, name::Symbol) = APSInstr(data, name, emptySymbol)
header(instr::APSInstr) = (instr.data[1] >> 8) % UInt8
opcode(instr::APSInstr) = header(instr) >> 4

function show(io::IO, e::APSInstr)
	PRINTSHORT || print(io, string(OPCODES[opcode(e)], ": "))
	for b in e.data
		print(io, hex(b, 4))
	end
end

function write_file(filename, seq)
	global PRINTSHORT = true
	open(filename, "w") do f
		for s = seq
			println(f, s)
		end
	end
	PRINTSHORT = false
end

function write_hdf5_file(filename, seq, wfA, wfB)

	if isfile(filename)
		rm(filename)
	end

	#Stack the instruction data
	instrs = UInt64[]
	for s in seq
		r = UInt64(0)
		for ct in 1:4
			r += UInt64(s.data[ct]) << 16*(4-ct)
		end
		push!(instrs, r)
	end

	#write it to the h5 file
	h5write(filename, "chan_1/instructions", instrs)

	#write the waveform data
	h5write(filename, "chan_1/waveforms", wfA);
	h5write(filename, "chan_2/waveforms", wfB);
end

function read_hdf5_file(filename)
	raw_instrs = h5read(filename, "chan_1/instructions")
	instrs = Array(APSInstr, length(raw_instrs))
	for ct in 1:length(raw_instrs)
		data = raw_instrs[ct]
		# unstack the data
		instrs[ct] = APSInstr(UInt16[(data >> 16*(4-i)) & 0xffff for i in 1:4])
	end
	return instrs
end

# instruction encodings
WFM = 0x0000
MARKER = 0x0001
WAIT_TRIG = 0x0002
LOAD_REPEAT = 0x0003
DEC_REPEAT = 0x0004
CMP = 0x0005
GOTO = 0x0006
CALL = 0x0007
RET = 0x0008
SYNC = 0x0009
# PREFETCH = 0x000A
LOAD_CMP = 0x000B

#The top bits of the 48 bit instruction payload set the sync or wait for trigger in the playback engines
WAIT_TRIG_INSTR_WORD = 1 << 14;
SYNC_INSTR_WORD = 1 << 15;
PREFETCH_INSTR_WORD = 3 << 14;

# CMP op encodings
EQUAL = 0x0000
NOTEQUAL = 0x0001
GREATERTHAN = 0x0002
LESSTHAN = 0x0003

#Waveform payload slv is 47 downto 0
# Bits 47-46 are op-code (play wf / wait for trig / wait for sync )
# Bit 45 is the TAPair flag
# Bits 44 downto 24 is count
# Bits 23 downto 0 is addr
function wf_entry(addr, count; TAPair=false, writeData=true, label=emptySymbol)
	@assert (addr & 0xffffff) == addr "Oops! We can only handle 24 bit addresses"
	@assert (count & 0xfffff) == count "Oops! We can only handle 20 bit counts"

	APSInstr(UInt16[WFM << 12 | (writeData ? 1 : 0) << 8, (TAPair ? 1 : 0) << 13  | count >> 8 , (count & 0x00ff) << 8 | addr >> 16, addr & 0xffff], label)
end

wf_prefetch(addr) = APSInstr(UInt16[WFM << 12 | 1 << 8, PREFETCH_INSTR_WORD, addr >> 16, addr & 0xffff])

function marker_entry(target, state, count, transitionWord; writeData=true, label=emptySymbol)
	@assert (count & 0xffffffff) == count "Oops! We can only handle 32 bit trigger counts"

	APSInstr(UInt16[MARKER << 12 | (writeData ? 1 : 0) << 8 | (target & 0x3) << 10, (transitionWord & 0xf) << 1 | (state & 0x1), count >> 16, count & 0xffff], label)
end

flow_entry(cmd, jumpAddr::Integer, payload; label=emptySymbol) =
	APSInstr(UInt16[(cmd << 12), payload, jumpAddr >> 16, jumpAddr & 0xffff], label)

flow_entry(cmd, jumpAddr::Symbol, payload; label=emptySymbol) =
	APSInstr(UInt16[(cmd << 12), payload, 0x0000, 0x0000], label, jumpAddr)

function wait_trig(; label=emptySymbol)
	instr = flow_entry(WAIT_TRIG, 0, WAIT_TRIG_INSTR_WORD; label=label)
	instr.data[1] |= 1 << 8; # set broadcast bit
	return instr
end

function sync(; label=emptySymbol)
	instr = flow_entry(SYNC, 0, SYNC_INSTR_WORD; label=label)
	instr.data[1] |= 1 << 8; # set broadcast bit
	return instr
end

function load_repeat(loadValue; label=emptySymbol)
	@assert loadValue < 2^16 "Oops! We can only handle 16 bit repeat counts."
	flow_entry(LOAD_REPEAT, loadValue, 0; label=label)
end

flow_loadcmp(; label=emptySymbol) = flow_entry(LOAD_CMP, 0, 0; label=label)
flow_cmp(cmpOp, mask; label=emptySymbol) = APSInstr(UInt16[CMP << 12, 0, 0, ((cmpOp & 0x03) << 8) | (mask & 0xff)], label)
flow_repeat(target; label=emptySymbol) = flow_entry(DEC_REPEAT, target, 0; label=label)
flow_goto(target; label=emptySymbol) = flow_entry(GOTO, target, 0; label=label)
flow_call(target; label=emptySymbol) = flow_entry(CALL, target, 0; label=label)
flow_return(; label=emptySymbol) = flow_entry(RET, 0, 0; label=label)

# updates the address in a command instruction
function updateAddr!(entry, addr)
	entry.data[3] = addr >> 16
	entry.data[4] = addr & 0xffff
	entry
end

function resolve_symbols!(seq)
	labeledEntries = filter((x) -> x[2].name != emptySymbol, collect(enumerate(seq)))
	symbolDict = [entry.name => idx-1 for (idx, entry) in labeledEntries]
	println(symbolDict)
	for entry in seq
		if entry.target != emptySymbol && haskey(symbolDict, entry.target)
			#println("Updating target of $(entry.target) to $(symbolDict[entry.target]). Before:")
			#println(entry)
			updateAddr!(entry, symbolDict[entry.target])
			#println("After:")
			#println(entry)
		end
	end
	seq
end

# # test sequences

# function simpleSeq()
# 	seq = APSInstr[]
# 	push!(seq, wf_entry(0, 63; label=:A))
# 	push!(seq, wf_entry(0, 7; label=:B))
# 	push!(seq, wf_entry(8, 7; label=:C))
# 	push!(seq, wf_entry(16, 7))
# 	push!(seq, wf_entry(24, 7))
# 	push!(seq, wf_entry(32, 7))
# 	push!(seq, flow_entry(GOTO, 0, 0, 0))
# 	resolve_symbols!(seq)
# end
function simpleSeq()
	seq = APSInstr[
		sync(),
		load_repeat(1),
		wf_entry(0x1A, 3; label=:A),
		flow_loadcmp(),
		flow_cmp(EQUAL, 0x01),
		flow_call(:C),
		wf_entry(0x2A, 3),
		flow_repeat(:A),
		load_repeat(2),
		wf_entry(0x3A, 3; label=:B),
		wf_entry(0x4A, 4; label=:C),
		wf_entry(0x5A, 4),
		wf_entry(0x6A, 3),
		flow_repeat(:B),
		wf_entry(0x00, 3; label=:D),
		flow_loadcmp(),
		flow_cmp(EQUAL, 0x02),
		flow_return(),
		flow_goto(:D)
	]
	resolve_symbols!(seq)
end

function ramsey()
	seq = APSInstr[]
	for delay in 10:10:100
		push!(seq, sync())
		push!(seq, wait_trig())
		push!(seq, wf_entry(0, 9, TAPair=true, writeData=false))
		push!(seq, marker_entry(0, 0, 8, 0xC))
		push!(seq, wf_entry(13, 11))
		push!(seq, marker_entry(0, 1, 11, 0x3))
		push!(seq, wf_entry(0, delay, TAPair=true))
		push!(seq, marker_entry(0, 0, delay, 0xC))
		push!(seq, wf_entry(1, 11))
		push!(seq, marker_entry(0, 1, 11, 0x3))
		push!(seq, marker_entry(0, 0, 1, 0x0))
	end
	push!(seq, flow_goto(0x00))
	resolve_symbols!(seq)
end

function cpmg()
	seq = APSInstr[]
	for (ct,rep) in enumerate(2.^[1:5])
		push!(seq, sync())
		push!(seq, wait_trig())

		#First 90
		push!(seq, wf_entry(0, 9, TAPair=true, writeData=false))
		push!(seq, marker_entry(0, 0, 8, 0xC))
		push!(seq, wf_entry(1, 11))
		push!(seq, marker_entry(0, 1, 11, 0x3))

		#CPMG block repeats
		push!(seq, load_repeat(rep-1))

		push!(seq, wf_entry(0, 0xA, TAPair=true; label=symbol("C$ct")))
		push!(seq, marker_entry(0, 0, 9, 0xC))
		push!(seq, wf_entry(25, 11))
		push!(seq, marker_entry(0, 1, 11, 0x3))
		push!(seq, wf_entry(0, 0x14, TAPair=true))
		push!(seq, marker_entry(0, 0, 0x13, 0xC))
		push!(seq, wf_entry(25, 11))
		push!(seq, marker_entry(0, 1, 11, 0x3))
		push!(seq, wf_entry(0, 0xA, TAPair=true))
		push!(seq, marker_entry(0, 0, 9, 0x0))
		push!(seq, flow_repeat(symbol("C$ct")))

		#Final 90
		push!(seq, wf_entry(1, 11))
		push!(seq, marker_entry(0, 1, 11, 0x3))
	end
	push!(seq, flow_goto(0x00))
	resolve_symbols!(seq)
end

function rabi_amp()
	seq = APSInstr[]
	for ct in 1:100
		push!(seq, sync())
		push!(seq, wait_trig())
		push!(seq, wf_entry(0, 9, TAPair=true, writeData=false))
		push!(seq, marker_entry(0, 0, 8, 0xC))
		push!(seq, wf_entry((ct-1)*256 + 1, 255))
		push!(seq, marker_entry(0, 1, 255, 0x3))
		push!(seq, wf_entry(0, 9, TAPair=true))
		push!(seq, marker_entry(0, 0, 9, 0x0))
	end
	push!(seq, flow_goto(0x00))
end

function test_prefetch()
	seq = APSInstr[]
	offset = 0
	for offset = 0:2^15:(3*2^15)
		push!(seq, wf_prefetch(offset))
		for ct in 1:100
			push!(seq, sync())
			push!(seq, wait_trig())
			push!(seq, wf_entry(offset+0, 9, TAPair=true, writeData=false))
			push!(seq, marker_entry(0, 0, 8, 0xC))
			push!(seq, wf_entry(offset+ (ct-1)*256 + 1, 255))
			push!(seq, marker_entry(0, 1, 255, 0x3))
			push!(seq, wf_entry(offset+0, 9, TAPair=true))
			push!(seq, marker_entry(0, 0, 9, 0x0))
		end
	end
	push!(seq, flow_goto(0x00))
end

function writeTestFiles()
	seq = simpleSeq();
	write_file("simpleSeq.dat", seq)
end

end


text:
x = -1 + 0im
y = sqrt(x)
println(y)

max_stars_repo_path: Basic Part/complex_compute.jl

max_issues_repo_path: Basic Part/complex_compute.jl

max_forks_repo_path: Basic Part/complex_compute.jl


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -12.266285154621722 <= q <= 12.793564807427844)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.46956444239618067 + 0.18451073354487368*$(x[1]) + -0.3056648113863898*$(x[2]) + -0.8733496775879632*$(x[3]) + -0.7482422964857776*$(x[4]) + 0.39599386992219143*$(x[5])))) + 1/(1 + exp(-(0.2831226968799312 + -0.000968797402465249*$(x[1]) + -0.6607077518367617*$(x[2]) + 0.11528579207894962*$(x[3]) + -0.7069435392252328*$(x[4]) + 0.7900850540395843*$(x[5])))) + 1/(1 + exp(-(0.41681338081805386 + 0.17922182588365843*$(x[1]) + 0.5927695896331957*$(x[2]) + -0.635120363897093*$(x[3]) + 0.31189019033409826*$(x[4]) + 0.11999673266347255*$(x[5])))) + 1/(1 + exp(-(-0.4140404108339504 + 0.583322476422016*$(x[1]) + 0.9921442381523282*$(x[2]) + 0.45546141054937284*$(x[3]) + -0.3096177212085385*$(x[4]) + 0.6860556002335301*$(x[5])))) + 1/(1 + exp(-(0.44730860193520705 + -0.9763312564707065*$(x[1]) + 0.297857633230612*$(x[2]) + -0.7735590056574844*$(x[3]) + -0.3545666189247747*$(x[4]) + -0.48025799425370996*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Robotlib
# using Robotlib.Calibration
using DSP # For filtfilt
include("../dynamics.jl")
include("gravityFridaLS.jl")
include("../gravityFridaLS2.jl")
# import JuMP
# using Ipopt
h = 0.004032;

dh = DHYuMi()
# xi = DH2twistsPOE(dh)

pathopen = "/work/fredrikb/extRosetta/frida_gravity_2.txt"
pathsave = "/work/fredrikb/log.mat"

# pathopen = "frida_gravity_3.txt"
# pathsave = "log.mat"

# data    = orcalog2mat(pathopen, pathsave)
data    = readmat(pathsave);
ds      = 1
q       = getData("robot_1.*posRawAbs", data, ds);
q̇       = getData("irb2ext_robot_1.*velRef", data, ds);
τ       = getData("robot_1.*trqRaw", data, ds);
f       = getData("force", data, ds);

abb2logical!(q);
abb2logical!(q̇);
abb2logical!(τ);
q = q*dh.GR';
q̇ = q̇*dh.GR';
τ = τ*inv(dh.GR');

q̈ = filtfilt(ones(50),[50.],smartDiff(q̇));
q̇ = filtfilt(ones(10),[10.],q̇);
# plot(abs([q̇, q̈]))

lowAcc  = all(abs(q̈) .< 3e-4,2);
q       = q[lowAcc,:];
q̇       = q̇[lowAcc,:];
τ       = τ[lowAcc,:];
f       = f[lowAcc,:];
N       = size(q,1)

baseAnglesLeft  = [-0.63 , 0.95 , -0.18]
Rbase           = rpy2R(baseAnglesLeft,"xyz")
Tbase           = eye(4)
Tbase[1:3,1:3]  = Rbase
T  = cat(3,[Tbase*fkinePOE(xi,q[i,:]') for i = 1:N]...);

# include("../dynamics.jl")

N = size(q,1)
np = 28+28
function getRegressor(q,q̇)
    N = size(q,1)
    A = Array(Float64,7N,np)
    ii = 1
    for i = 1:N
        A[ii:ii+6,:] = [gravityFridaLS(q[i,:])  diagm(q̇[i,:][:] .> 0)   diagm(q̇[i,:][:] .< 0) diagm(q̇[i,:][:] .> 0)*diagm(q̇[i,:][:])  diagm(q̇[i,:][:] .< 0)*diagm(q̇[i,:][:])]
        ii+= 7
    end
    return A
end

@time A = getRegressor(q,q̇);

still = vec(abs(q̇)') .< 2e-4
still[1:3:end] = true
still[2:3:end] = true

estA(A) = A[!still,:]
esty(y) = vec(y')[!still]

k = estA(A)\esty(τ)
# k = [estA(A); 0.1eye(np)]\[esty(τ);zeros(np)]
τhat = reshape(A*k,7,N)'
err = τ-τhat
println("Error: ", rms(esty(err)))

# println("Estimating a position dependent friction model for different regularization parameters")
# include("../posDepFric.jl")
# include("/local/home/fredrikb/SystemIdentification/src/least_squares.jl")
# n_basis = [10, 4]
# nppdf = prod(n_basis)
# centers = getCenters(n_basis, q, q̇)
# @time Apdf = frictionRBFN(q, q̇, centers, normalized=true);
# trials = 30
# λ = logspace(-15,6,trials)
# @time kpdf,errpdf,normk = dsvd(estA(Apdf),esty(err),λ)
# for i = 1:trials
#     τpdf = reshape(Apdf*kpdf[:,i],7,N)';
#     errpdf[i] = rms(esty(τ-τhat-τpdf))
# end
# Lcurve(errpdf, normk, λ)
# τpdf = reshape(Apdf*kpdf[:,end-8],7,N)';

plot(τ,layout=7)
plot!(τhat,c=:red,layout=7)
# subplot!(τhat + τpdf,c=:green,nr=7)


text:
#Code to extract genes from NCBI database using their API
#API found here: https://api.ncbi.nlm.nih.gov/variation/v0/

#takes "rs" out of snpids/rsids if they start with "rs###"
function getSNPID(snpid::String)
    id = ""
    try
        id = match(r"(?<=rs)\w+", snpid).match
    catch
        id = snpid
    end
    return id
end

#checks if gene information is present, if not, return no gene found
function genecheck(info::LazyJSON.Object{Nothing,String})
    locus = ""
    try
        locus = info["primary_snapshot_data"]["allele_annotations"][1]["assembly_annotation"][1]["genes"][1]["locus"]
    catch
        locus = "No gene listed"
    end
    return locus
end

#checks to make sure the HTTP has valid information
function httpcheck(snpid::AbstractString)
    http = HTTP.Messages.Response()
    try
        http = HTTP.get("https://api.ncbi.nlm.nih.gov/variation/v0/beta/refsnp/" * snpid)
    catch
        http = HTTP.Messages.Response()
    end
    return http
end

"""
    getgenes(data::DataFrame; idvar::AbstractString)

# Position arguments

- `data::DataFrame`: A DataFrame containing a column with the Ref SNP IDs. By default, assumes that the variable name is "snpid". The variable name can be specified using the `idvar` keyword.

# Keyword arguments

- `idvar::AbstractString`: the variable name in the dataframe that specifies the Ref SNP ID (rsid).

    getgenes(snps::AbstractArray)

# Position arguments

- `snps::AbstractArray`: Ref SNP IDs (rsid) to get loci names for.

# Output

Returns an array of gene loci associated to the Ref SNP IDs.

"""
function getgenes(snpids::AbstractArray)
    snpids = getSNPID.(snpids)
    loci = Vector{String}(undef, 0)
    for snpid in snpids
        http = httpcheck(snpid)
        if isempty(http.body)
            locus = "snpid not in database"
        else
            str = String(http.body)
            info = LazyJSON.parse(str)
            locus = genecheck(info)
        end
        push!(loci, locus)
    end
    return loci
end

function getgenes(df::DataFrame; idvar::AbstractString = "snpid")
    rsidsym = Meta.parse(idvar)
    if !(rsidsym in names(df))
        throw(ArgumentError(idvar * " is not in the dataframe. Please rename
        the column of rsids to `snpid` or specify the correct name using the `idvar` argument"))
    end
    getgenes(df[rsidsym])
end

function getgenes(snpid::Union{String, Int64})
    if typeof(snpid) == Int64
        snpid = string(snpid)
    end
    snpid = getSNPID(snpid)
    locus = ""
    http = httpcheck(snpid)
    if isempty(http.body)
        locus = "snpid not in database"
    else
        str = String(http.body)
        info = LazyJSON.parse(str)
        locus = genecheck(info)
    end
    return locus
end

#gets additional gene information
function getinfo(info)
    snp_id = info["refsnp_id"]
    annotations = info["primary_snapshot_data"]["allele_annotations"]
    geneinfo = Dict()
    for a in annotations
        assembly_annotations = a["assembly_annotation"]
        for aa in assembly_annotations
            seq_id = aa["seq_id"]
            annotation_release = aa["annotation_release"]
            genes = aa["genes"]
                for g in genes
                    gene_name = g["name"]
                    gene_id = string(g["id"])
                    gene_locus = g["locus"]
                    gene_is_pseudo = string(Int(g["is_pseudo"] == "true"))
                    gene_orientation = string(Int(g["orientation"] == "minus"))
                    geneinfo = Dict("seq_id" => seq_id, "annotation_release" => annotation_release,
                    "gene_name" => gene_name, "gene_id" => gene_id, "gene_locus" => gene_locus,
                    "gene_is_pseudo" => gene_is_pseudo, "gene_orientation" => gene_orientation)
                end
        end
    end
    return geneinfo
end

"""
    getgeneinfo(snpid::AbstractString)

# Position arguments

- `snp::AbstractString`: Ref SNP ID (rsid) to get gene information for.

# Output

Returns a dictionary of gene information associated with the Ref SNP ID. Entries are `seq_id`, `annotation_release`, `gene_name`, `gene_id`, `gene_locus`, `gene_is_pseudo`, and `gene_orientation`.

"""
function getgeneinfo(snpid::AbstractString)
    geneinfo = Dict()
    rsid = getSNPID(snpid);
    http = httpcheck(rsid);
    if isempty(http.body)
        throw(ArgumentError("Ref SNP ID (rsid) not in database"))
    end
    str = String(http.body)
    info = LazyJSON.parse(str)
    geneinfo = getinfo(info)
    if isempty(geneinfo)
        geneinfo = Dict("Empty" => "No gene info listed for ref SNP ID")
    end
    return geneinfo
end


text:
module Readme

using LightGraphs
using GraphPlot, GraphPlot.Compose
using StorageGraphs
using Test
using Logging

log = SimpleLogger(stdout, Logging.Debug)

function draw_graph(g, x, y, name::String, args...; ns=1, C=5)
    layout = (x...)->spring_layout(x...; C=C)
    draw(SVG("$(@__DIR__)/../assets/$name.svg", x*cm, y*cm),
        plot_graph(g, layout=layout, nodesize=ns))
end

g = StorageGraph()
@test StorageGraphs.add_vertex!.(Ref(g), [(x=1,),(x=2,),(x=3,)]) |> all
@test get_prop(g, 1) == (x=1,)
@test get_prop(g, 2) == (x=2,)
@test get_prop(g, 3) == (x=3,)

draw_graph(g, 10, 4, "ex1")

# with_logger(log) do
#     add_derived_values!(g, (x=[1,2,3],), (y=[1,4,9],))
# end

g = StorageGraph()
add_nodes!(g, (x=[1,2,3],)=>(y=[1,4,9],))

@test nv(g) == 2
@test ne(g) == 1
@test get_prop(g, 1) == (x=[1,2,3],)
@test get_prop(g, 2) == (y=[1,4,9],)

draw_graph(g, 10, 4, "ex2", C=8)

# Code snippets

using StorageGraphs

g = StorageGraph()
add_nodes!(g, (x=[1,2,3],)=>(y=[1,4,9],))

###

@test isempty(setdiff(g[:x][1], [1,2,3]))

###

using StorageGraphs

g = StorageGraph()

add_nodes!(g, (P=1,)=>(alg="alg1",)=>(x=[10., 20., 30.],))

draw_graph(g, 10, 4, "ic_graph", C=8)

simulation(x; alg) = alg == "alg1" ? x.+2 : x.^2

# retrieve the previously stored initial conditions
x = g[:x, (P=1,)=>(alg="alg1",)][1]
# there is only one node depending on (P=1,)=>(alg="alg1",), so that's why
# we take only the first element
results = simulation(x, alg="alg1")
add_nodes!(g, foldr(=>, ((P=1,), (alg="alg1",), (x=x,), (r=results,))))
# foldr(=>, dep) can be useful for building long dependency chanis form parts

draw_graph(g, 10, 6, "sim_graph", C=8)

add_nodes!(g, foldr(=>,((P=2,), (alg="alg1",), (x=2x,), (r=2results,))))

draw_graph(g, 11, 11, "complicated_graph", C=6)

end  # module Readme

@testset "Readme" begin
    using .Readme
end


text:
# Filters
# =======

# Creation of a boolean mask array
# --------------------------------

# If the input is a matrix, dropdims the singleton dimension

function _column_mask(mask::AbstractMatrix{Bool}, msa)
    @assert size(mask, 1) == 1 "The mask should be a vector or a matrix of size (1,ncol)"
    @assert size(mask, 2) == ncolumns(msa) "One boolean value per column is needed."
    dropdims(mask, dims=1)
end

function _sequence_mask(mask::AbstractMatrix{Bool}, msa)
    @assert size(mask, 2) == 1 "The mask should be a vector or a matrix of size (nseq,1)"
    @assert size(mask, 1) == nsequences(msa) "One boolean value per sequence is needed."
    dropdims(mask, dims=2)
end

# If the input is a function, apply it to the respective slice

function _column_mask(mask::Function, msa)
    out = mapslices(mask, msa, dims=1)
    if size(out) != (1, ncolumns(msa)) || eltype(out) != Bool
        error("The function must return a Bool element per column.")
    end
    dropdims(out, dims=1)
end

function _sequence_mask(mask::Function, msa)
    out = mapslices(mask, msa, dims=2)
    if size(out) != (nsequences(msa), 1) || eltype(out) != Bool
        error("The function must return a Bool element per column.")
    end
    dropdims(out, dims=2)
end

# If the mask is an AbstractVector{Bool}, return it without changes

function _column_mask(mask::AbstractVector{Bool}, msa)
    @assert length(mask) == ncolumns(msa) "One boolean value per column is needed."
    mask
end

function _sequence_mask(mask::AbstractVector{Bool}, msa)
    @assert length(mask) == nsequences(msa) "One boolean value per sequence is needed."
    mask
end

# If the mask is list of indexes, return it without changes

_column_mask(mask::Union{AbstractArray, Colon}, msa) = mask

# Annotate modifications
# ----------------------

_annotate_modification!(x, selector::Colon, fun::String, column::Bool) = nothing

function _annotate_modification!(x, selector::AbstractArray, fun::String, column::Bool)
    entity = column ? "column" : "sequence"
    verb = "has"
    if eltype(selector) <: Bool
        n_deleted = sum(.~selector)
        if n_deleted > 0
            if n_deleted > 1
                entity = "s"
                verb = "have"
            end
            str = string("$fun : $n_deleted $entity $verb been deleted.")
            annotate_modification!(x, str)
        end
    else
        cleaned_selector = unique(selector)
        n_selected = length(cleaned_selector)
        if n_selected > 0
            if n_selected > 1
                entity *= "s"
                verb = "have"
            end
            str = string("$fun : $n_selected $entity $verb been selected.")
            annotate_modification!(x, str)
            if column && n_selected > 1 && !issorted(selector) # do not store sequence order changes
                annotate_modification!(x, "$fun : column order has changed!")
            end
        end
    end
end

function _annotate_col_modification!(x, selector)
    _annotate_modification!(x, selector, "filtercolumns!", true)
end

function _annotate_seq_modification!(x, selector)
    _annotate_modification!(x, selector, "filtersequences!", false)
end

# Filter sequences
# ----------------

"It's similar to `filtersequences!` but for an `AbstractMatrix{Residue}`"
filtersequences(msa::AbstractMatrix{Residue}, mask) = msa[_sequence_mask(mask, msa), :]

"""
`filtersequences!(msa, mask[, annotate::Bool=true])`

It allows to filter `msa` sequences using a `AbstractVector{Bool}` `mask`
(It removes sequences with `false` values). `AnnotatedMultipleSequenceAlignment` annotations
are updated if `annotate` is `true` (default).
"""
function filtersequences!(msa::AnnotatedMultipleSequenceAlignment,
                          mask, annotate::Bool=true)
    boolean_vector = _sequence_mask(mask, msa)
    filtersequences!(annotations(msa), sequencenames(msa), boolean_vector)
    # Filter annotations first, names will be changed by filtersequences:
    msa.matrix = filtersequences(namedmatrix(msa), boolean_vector)
    if annotate
        _annotate_seq_modification!(msa, boolean_vector)
    end
    msa
end

function filtersequences!(msa::MultipleSequenceAlignment,
                          mask, annotate::Bool=false) # annotate is useful for calling this
                                                      # inside other functions
    msa.matrix = filtersequences(namedmatrix(msa), _sequence_mask(mask, msa))
    msa
end

function filtersequences(msa::Union{AnnotatedMultipleSequenceAlignment,
                                    MultipleSequenceAlignment}, args...)
    filtersequences!(deepcopy(msa), args...)
end

# Filter columns
# --------------

"It's similar to `filtercolumns!` but for an `AbstractMatrix{Residue}`"
function filtercolumns(msa::AbstractMatrix{Residue}, mask)
    msa[:, _column_mask(mask, msa)]
end

"""
`filtercolumns!(msa, mask[, annotate::Bool=true])`

It allows to filter MSA or aligned sequence columns/positions using a
`AbstractVector{Bool}` `mask`. Annotations are updated if `annotate` is `true` (default).
"""
function filtercolumns!(x::AnnotatedAlignedObject,
                        mask, annotate::Bool=true)
    selector = _column_mask(mask, x)
    filtercolumns!(annotations(x), selector)
    x.matrix = filtercolumns(namedmatrix(x), selector)
    if annotate
        _annotate_col_modification!(x, selector)
    end
    x
end

function filtercolumns!(x::UnannotatedAlignedObject,
                        mask, annotate::Bool=false) # annotate is useful for calling this
                                                    # inside other functions
    x.matrix = filtercolumns(namedmatrix(x), _column_mask(mask, x))
    x
end

filtercolumns(x::AbstractAlignedObject, args...) = filtercolumns!(deepcopy(x), args...)

# Util function
# -------------

"""
This function takes a vector of sequence names and a sequence id.
It returns the position of that id in the vector.
If the id isn't in the vector, It throws an error.
"""
function _get_seqid_index(names::Vector{String}, sequence_id::String)
    id_index = something(findfirst(isequal(sequence_id), names), 0)
    if id_index == 0
        throw(ErrorException("$sequence_id is not in the list of sequence names."))
    end
    id_index
end

# Reference
# ---------
"It swaps the names on the positions `i` and `j` of a `Vector{String}`"
function _swap!(names::Vector{String}, i::Int, j::Int)
    name = names[i,:]
    names[i,:] = names[j,:]
    names[j,:] = name
    names
end

"""
It swaps the sequences on the positions `i` and `j` of an MSA. Also it's possible to swap
sequences using their sequence names/identifiers when the MSA object as names.
"""
function swapsequences!(matrix::Matrix{Residue}, i::Int, j::Int)
    seq = matrix[i,:]
    matrix[i,:] = matrix[j,:]
    matrix[j,:] = seq
    return matrix
end

function swapsequences!(matrix::NamedArray, i::Int, j::Int)
    swapsequences!(getarray(matrix), i, j)
    setnames!(matrix, _swap!(sequencenames(matrix), i, j), 1)
    return matrix
end

function swapsequences!(matrix::NamedArray, i::String, j::String)
    seqnames = sequencenames(matrix)
    swapsequences!(matrix,
                   _get_seqid_index(seqnames,i),
                   _get_seqid_index(seqnames,j))
end

"""
It puts the sequence `i` (name or position) as reference (first sequence) of the MSA. This
function swaps the sequences 1 and `i`.
"""
function setreference!(msa::AnnotatedMultipleSequenceAlignment, i::Int, annotate::Bool=true)
    swapsequences!(namedmatrix(msa), 1, i)
    if annotate
        seqnames = sequencenames(msa)
        annotate_modification!(msa, string("setreference! : Using ",
                                            seqnames[1], " instead of ",
                                            seqnames[i], " as reference."))
    end
    msa
end

function setreference!(msa::MultipleSequenceAlignment, i::Int, annotate::Bool=false)
    # The annotate argument is useful for calling this inside other functions
    swapsequences!(namedmatrix(msa), 1, i)
    msa
end

function setreference!(msa::NamedResidueMatrix{T},
                       i::Int, annotate::Bool=false) where T <: AbstractArray
    swapsequences!(msa, 1, i)
end

function setreference!(msa::NamedResidueMatrix{T},
                       id::String, annotate::Bool=false) where T <: AbstractArray
    swapsequences!(msa, names(msa,1)[1], id)
end

function setreference!(msa::AbstractMultipleSequenceAlignment, id::String,
                       annotate::Bool=true)
    setreference!(msa, _get_seqid_index(sequencenames(msa), id), annotate)
end

function setreference!(msa::Matrix{Residue}, i::Int, annotate::Bool=false)
    # The annotate argument is useful for calling this inside other functions
    swapsequences!(msa, 1, i)
end

"""
Creates a new matrix of residues. This function deletes positions/columns of the MSA with
gaps in the reference (first) sequence.
"""
function adjustreference(msa::AbstractMatrix{Residue}, annotate::Bool=false)
    # The annotate argument is useful for calling this inside other functions
    msa[:, msa[1,:] .!= GAP ]
end

"""
It removes positions/columns of the MSA with gaps in the reference (first) sequence.
"""
function adjustreference!(msa::AbstractMultipleSequenceAlignment, annotate::Bool=true)
    filtercolumns!(msa, vec(getresidues(getsequence(msa,1))) .!= GAP, annotate)
end

"""
Creates a new matrix of `Residue`s (MSA) with deleted sequences and columns/positions. The
MSA is edited in the following way:

 1. Removes all the columns/position on the MSA with gaps on the reference (first) sequence
 2. Removes all the sequences with a coverage with respect to the number of
 columns/positions on the MSA **less** than a `coveragelimit`
 (default to `0.75`: sequences with 25% of gaps)
 3. Removes all the columns/position on the MSA with **more** than a `gaplimit`
 (default to `0.5`: 50% of gaps)
"""
function gapstrip(msa::AbstractMatrix{Residue}; coveragelimit::Float64=0.75,
                  gaplimit::Float64=0.5)
    msa = adjustreference(msa)
    # Remove sequences with pour coverage of the reference sequence
    if ncolumns(msa) != 0
        msa = filtersequences(msa, vec(coverage(msa) .>= coveragelimit))
    else
        throw("There are not columns in the MSA after the gap trimming")
    end
    if nsequences(msa) != 0
        msa = filtercolumns(msa, vec(columngapfraction(msa) .<= gaplimit))
    else
        throw("There are not sequences in the MSA after coverage filter")
    end
    msa
end

"""
This functions deletes/filters sequences and columns/positions on the MSA on the following
order:

1. Removes all the columns/position on the MSA with gaps on the reference (first) sequence.
2. Removes all the sequences with a coverage with respect to the number of
columns/positions on the MSA **less** than a `coveragelimit`
(default to `0.75`: sequences with 25% of gaps).
3. Removes all the columns/position on the MSA with **more** than a `gaplimit`
(default to `0.5`: 50% of gaps).
"""
function gapstrip!(msa::AbstractMultipleSequenceAlignment,
                   annotate::Bool=isa(msa,AnnotatedAlignedObject);
                   coveragelimit::Float64=0.75, gaplimit::Float64=0.5)
    if annotate
        annotate_modification!(msa,
            string("gapstrip! : Deletes columns with gaps in the first sequence."))
    end
    adjustreference!(msa, annotate)
    # Remove sequences with poor coverage of the reference sequence
    if ncolumns(msa) != 0
        if annotate
            annotate_modification!(msa,
                string("gapstrip! : Deletes sequences with a coverage less than ",
                coveragelimit))
        end
        filtersequences!(msa, vec(coverage(msa) .>= coveragelimit), annotate)
    else
        throw("There are not columns in the MSA after the gap trimming")
    end
    # Remove columns with a porcentage of gap greater than gaplimit
    if nsequences(msa) != 0
        if annotate
            annotate_modification!(msa,
                string("gapstrip! : Deletes columns with more than ", gaplimit, " gaps."))
        end
        filtercolumns!(msa, vec(columngapfraction(msa) .<= gaplimit), annotate)
    else
        throw("There are not sequences in the MSA after coverage filter")
    end
    msa
end

"""
This function is similar to `gapstrip!' but does not filter on reference

This functions deletes/filters sequences and columns/positions on the MSA on the following
order:

1. Removes all the sequences with a coverage with respect to the number of
columns/positions on the MSA **less** than a `coveragelimit`
(default to `0.75`: sequences with 25% of gaps).
2. Removes all the columns/position on the MSA with **more** than a `gaplimit`
(default to `0.5`: 50% of gaps).
"""
function gapstripnoref!(msa::AbstractMultipleSequenceAlignment,
                   annotate::Bool=isa(msa,AnnotatedAlignedObject);
                   coveragelimit::Float64=0.75, gaplimit::Float64=0.5)

    # Remove sequences with poor coverage of the reference sequence
    if ncolumns(msa) != 0
        if annotate
            annotate_modification!(msa,
                string("gapstripnoref! : Deletes sequences with a coverage less than ",
                coveragelimit))
        end
        filtersequences!(msa, vec(coverage(msa) .>= coveragelimit), annotate)
    else
        throw("There are not columns in the MSA after the gap trimming")
    end
    # Remove columns with a porcentage of gap greater than gaplimit
    if nsequences(msa) != 0
        if annotate
            annotate_modification!(msa,
                string("gapstripnoref! : Deletes columns with more than ", gaplimit, " gaps."))
        end
        filtercolumns!(msa, vec(columngapfraction(msa) .<= gaplimit), annotate)
    else
        throw("There are not sequences in the MSA after coverage filter")
    end
    msa
end


text:
# # Introduction to DataFrames.jl

# Loosely based on [this
# tutorial](https://juliaai.github.io/DataScienceTutorials.jl/data/dataframe/)
# from Data Science Tutorials.

# DataFrames.jl **cheatsheets**:
# [english](https://www.ahsmart.com/assets/pages/data-wrangling-with-data-frames-jl-cheat-sheet/DataFramesCheatSheet_v1.x_rev1.pdf),
# [chinese](https://www.ahsmart.com/assets/pages/data-wrangling-with-data-frames-jl-cheat-sheet/DataFramesCheatSheet_v1.x_rev1_zh.pdf)

# ## Setup

# Inspect Julia version:

VERSION

# The following instantiates a package environment.

# The package environment has been created using **Julia 1.6** and may not
# instantiate properly for other Julia versions.

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# ## Loading some demonstration data

# We'll be using [OpenML](https://www.openml.org/home) to grab datasets.

using OpenML
using DataFrames

table = OpenML.load(42638); # Titanic data set
typeof(table)

# This is not a `DataFrame`. However, like
# [many](https://github.com/JuliaData/Tables.jl/blob/main/INTEGRATIONS.md)
# other tabular data containers in Julia, it can be converted to one in the
# obvious way:

df = DataFrame(table);

# A `DataFrame` is essentially just a wrapper around a number of
# vectors with names, conceptualized as a table with the vectors as
# columns.

# Lets' look the first few rows (observations) of `df`:

first(df, 4)

# Here's a summary of stats for each column:

describe(df)

# To get just the column names, use `names(df)`.

# ## Selecting a single entry:

number = df[5, :age]   # or df[5, 4]

# ## Selecting a single row:

# We use the wild-card `:`

row = df[5, :]

#-

number = row.age

# ## Selecting multiple rows

# By index:

small_df = df[3:7, :]

# By applying some criterion:

df2 = filter(df) do row
    row.age > 60 && row.survived == "1"
end

# (Alternatively, one can use `subset` or Boolean indexing; see the
# [cheatsheet](https://www.ahsmart.com/assets/pages/data-wrangling-with-data-frames-jl-cheat-sheet/DataFramesCheatSheet_v1.x_rev1.pdf).)

# ## Selecting a single column

# Here's how I get the column named "Age":

df2.age

# Since this is just a vector, I can restrict the rows as usual:

df2.age[2:5]  # but same as df2[2:5, :Age]

# ## Selecting multiple columns - `select`

# The following wraps `df.age` as a single-column `DataFrame`:

select(df2, :age)

# ... which generalizes to multiple columns:

select(df2, [:age, :sex])

#-

select(df2, Not(:age))

#-

select(df2, Between(:sex, :fare))

# ## Copies or views of the data?

# So far, all the selection mechanisms discussed so far return a new
# `DataFrame` object, with these exceptions:

# * selecting a single row, as in `df[3,:]` - this returns a `DataFrameRow`.
# * selecting a single column as a vector, as in `df.age` - this returns a vector

# These exceptions are *views* of the original `DataFrame`, in the
# sense that no data is copied and mutating the view mutates the
# original `DataFrame`:

v = df.age;
v[1] = 1000
df[1, :age]

# Mutating the output of all the other selection mechanisms will not
# touch the original `DataFrame`, as data is copied in those cases.

# Use `select!` to make "in-place" column selections which mutate the
# original `DataFrame`. For example, the following permanently drops a
# column from `df`:

select!(df, Not(:cabin));
first(df, 2)

# In addition to row and column views of a `DataFrame` we can construct
# larger views, of type `SubDataFrame`, which share most of the
# behviour of an ordinary `DataFrame`, except they only *point* to
# data in the parent `DataFrame`, rather than duplicate it:

df_copy = df[1:3, [:fare, :age]]       # this is copy
df_view = @view df[1:3, [:fare, :age]] # this is a view

#-

df_view[1, :age] = 4000
df[1, :age]

#-

df_copy[1, :age]

# ## Describing the data

# To broaden the summary provided by `describe(df)` we can pass a
# number of symbols to indicate the statistics we are after:

describe(df, :min, :max, :mean, :median, :std)

# The following are all supported:
# * `mean`, `std`, `min`, `max`, `median`, `first`, `last` are all fairly self explanatory
# * `q25`, `q75` are respectively for the 25th and 75th percentile,
# * `eltype`, `nunique`, `nmissing` can also be used

# You can also pass your custom function with a pair `function =>
# :name` for instance:

using Statistics # to get functions like `mean` and `std`
foo(v) = mean(abs.(v))
d = describe(df, :mean, :median, foo => :mae)
first(d, 3)

# Note that the object returned by `describe` is itself a `DataFrame`:

select(d, [:variable, :mean])

# ## Materializing as a matrix
#
# To convert the content of the dataframe as one big matrix do this:
# use `convert`:

mat = Matrix(df)
mat[1:3, 1:3]

# ## Adding columns

# Adding a column to a dataframe is very easy:

df.weird = df.fare ./ df.age;
first(df, 3)

# That's it! Remember also that you can drop columns or make
# subselections with `select` and `select!`.

# The remainder of this tutorial is an abbreviated version of a
# [DataFrames.jl
# tutorial](http://juliadata.github.io/DataFrames.jl/latest/man/split_apply_combine/).

iris = OpenML.load(61) |> DataFrame;
first(iris, 3)

# ## `groupby`
#
# The `groupby` function allows to form "sub-dataframes" corresponding
# to groups of rows.  This can be very convenient to run specific
# analyses for specific groups without copying the data.

# The basic usage is `groupby(df, cols)` where `cols` specifies one or
# several columns to use for the grouping.

# Consider a simple example: in `iris` there is a `Species` column
# with 3 species:

unique(iris.class)

# We can form views for each of these:

gdf = groupby(iris, :class);

# The `gdf` object now corresponds to **views** of the original
# dataframe for each of the 3 species; the first species
# is `"Iris-setosa"` with:

subdf_setosa = gdf[1];
typeof(subdf_setosa)

#-

describe(subdf_setosa, :min, :mean, :max)

# Recall that this means modifying `subdf_setosa` also modifies its
# parent `iris`.

# Do `?groupby` for more information.

# ## `combine`

# The `combine` function allows to derive a new dataframe out of
# transformations of an existing one.  Here's an example taken from
# the official doc (see `?combine`):

df = DataFrame(a=1:3, b=4:6)

#-

combine(df, :a => sum, nrow)

# What happened here is that the derived DataFrame has two columns
# obtained respectively by (1) computing the sum of the first column
# and (2) applying the `nrow` function on the `df`.
#
# The transformation can produce one or several values, and `combine` will
# try to accomodate this with appropriate packing:

foo(v) = v[1:2]
combine(df, :a => maximum, :b => foo)

# Here the maximum value of `a` is copied twice so that the two
# columns have the same number of rows.

bar(v) = v[end-1:end]
combine(df, :a => foo, :b => bar)

# ## `combine` with `groupby`
#
# Combining `groupby` with `combine` is very useful.  For instance you
# might want to compute statistics across groups for different
# variables:

combine(groupby(iris, :class), :petallength => mean)

# Let's break this operatioin down:
#
# 1. The `groupby(iris, :class)` creates groups using the `:class` column (which has values `Iris-setosa`, `Iris-versicolor`, `Iris-virginica`)
# 2. The `combine` creates a derived dataframe by applying the `mean` function to the `:petallength` column
# 3. Since there are three groups, we get one column (mean of `petallength`) and three rows (one per group).
#
#
# You can do this for several columns/statistics at the time and give
# new column names to the results:

gdf = groupby(iris, :class)
combine(gdf, :petallength => mean => :MPL, :petallength => std => :SPL)

# So here we assign the names `:MPL` and `:SPL` to the derived
# columns.  If you want to apply something on all columns apart from
# the grouping one, using `names` and `Not` comes in handy:

combine(gdf, names(iris, Not(:class)) .=> std)

# where

names(iris, Not(:class))

# and note the use of `.` in `.=>` to indicate that we broadcast the
# function over each column.


text:
module Plugin

import Pkg
import UUIDs
import Pkg.TOML
using Suppressor: @suppress_err

export plugins, load

struct PluginIterator
    name::String
    ctx::Pkg.Types.Context
    deps::Vector{UUIDs.UUID}
end

struct PluginEntry
    package::Symbol
    name::String
    expr::Expr
end

function load(plugin::PluginEntry)
    importexpr = Expr(:using, Expr(:., plugin.package))
    suppressed = Expr(:macrocall, Symbol("@suppress_err"), " ", importexpr)

    if plugin.expr.head == :toplevel
        code = copy(plugin.expr)
        insert!(code.args, 1, suppressed)
    else
        code = :($suppressed; $(plugin.expr))
    end

    return eval(code)
end

Base.IteratorSize(::Type{PluginIterator}) = Base.SizeUnknown()

function plugins(name::String)
    ctx = Pkg.Types.Context()
    deps = collect(keys(ctx.env.manifest))
    PluginIterator(name, ctx, deps)
end

function plugins(entry::Pkg.Types.PackageEntry, uuid::UUIDs.UUID)
    entry.tree_hash !== nothing || return nothing
    pkg_dir = Pkg.Operations.find_installed(entry.name, uuid, entry.tree_hash)
    toml_path = joinpath(pkg_dir, "Project.toml")
    isfile(toml_path) || return nothing
    get(TOML.parsefile(toml_path), "plugins", nothing)
end

function Base.iterate(iter::PluginIterator, state = (1, []))
    (count, cases) = state
    if length(cases) > 0
        return cases[1], (count, cases[2:end])
    end

    (entry, target_plugins) = (nothing, nothing)

    while target_plugins === nothing || length(target_plugins) == 0
        if length(iter.deps) < count
            return nothing
        end

        uuid = iter.deps[count]
        entry = iter.ctx.env.manifest[uuid]
        results = plugins(entry, uuid)

        if results !== nothing
            target_plugins = get(results, iter.name, nothing)
        end
        count = count + 1
    end

    cases = [PluginEntry(Symbol(entry.name), k, Meta.parse(v)) for (k,v) in target_plugins]

    return cases[1], (count, cases[2:end])
end

end # module


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    sdfilewriter

function atomblock(io::IO, mol::GraphMol)
    for (i, atom) in enumerate(nodeattrs(mol))
        # TODO: SMILES coords
        (x, y, z) = atom.coords
        xyzsym = @sprintf "%10.4f%10.4f%10.4f %-3s" x y z string(atom.symbol)
        println(io, "$(xyzsym) 0  0  0  0  0  0  0  0  0  0  0  0")
    end
    return
end

function bondblock(io::IO, mol::GraphMol)
    for (i, (u, v)) in enumerate(edgesiter(mol))
        bond = edgeattr(mol, i)
        # TODO: SmilesBond
        uv = @sprintf "%3d%3d%3d%3d  0  0  0" u v bond.order bond.notation
        println(io, uv)
    end
    return
end

function propertyblock(io::IO, mol::GraphMol)
    charges = Tuple{Int,Int}[]
    radicals = Tuple{Int,Int}[]
    masses = Tuple{Int,Float64}[]
    for (i, atom) in enumerate(nodeattrs(mol))
        atom.charge == 0 || push!(charges, (i, atom.charge))
        atom.multiplicity == 1 || push!(radicals, (i, atom.multiplicity))
        atom.mass === nothing || push!(masses, (i, atom.mass))
    end
    if !isempty(charges)
        head = @sprintf "M  CHG%3d" length(charges)
        print(io, head)
        for (i, chg) in charges
            rcd = @sprintf " %3d %3d" i chg
            print(io, rcd)
        end
        println(io)
    end
    if !isempty(radicals)
        head = @sprintf "M  RAD%3d" length(radicals)
        print(io, head)
        for (i, rad) in radicals
            rcd = @sprintf " %3d %3d" i rad
            print(io, rcd)
        end
        println(io)
    end
    if !isempty(masses)
        head = @sprintf "M  ISO%3d" length(masses)
        print(io, head)
        for (i, iso) in masses
            rcd = @sprintf " %3d %3d" i iso
            print(io, rcd)
        end
        println(io)
    end
    return
end

function datablock(io::IO, mol::GraphMol)
    for (key, val) in mol.attributes
        println(io, "> <$(string(key))>")
        println(io, string(val))
        println(io, "")
    end
    println(io, raw"$$$$")
    return
end

function molblock(io::IO, mol::GraphMol)
    # TODO: 2D coords for SmilesMol
    println(io)
    println(io, "MolecularGraph.jl version $(version())")
    println(io)
    ncnt = nodecount(mol)
    ecnt = edgecount(mol)
    header = @sprintf "%3d%3d  0  0  0  0  0  0  0  0999 V2000" ncnt ecnt
    println(io, header)
    atomblock(io, mol)
    bondblock(io, mol)
    propertyblock(io, mol)
    println(io, "M  END")
    datablock(io, mol)
    return
end

"""
    sdfilewriter(io::IO, mols)
    sdfilewriter(filename::AbstractString, mols)

Write molecule data to the output stream as a SDFile format file.
"""
sdfilewriter(io::IO, mols) = molblock.((io,), mols)
sdfilewriter(filename::AbstractString, mols
    ) = sdfilewriter(open(path, "w"), mols)


text:
module ConvHMM

using LinearAlgebra
using Distributions
using Parameters
using Random: shuffle
using Requires

import Distributions: logpdf, loglikelihood

include("utils.jl")
include("wavelets.jl")
include("convolvedhmm.jl")
include("likelihoods.jl")
include("approximator.jl")
include("forward.jl")

# optionally load GeoStats.jl API
function __init__()
  @require GeoStatsBase="323cb8eb-fbf6-51c0-afd0-f8fba70507b2" include("geostats.jl")
end

export
  ConvolvedHMM,
  logpdf,
  loglikelihood,
  sample,
  mcmc,

  # wavelets
  kernelmatrix,
  Ricker,
  Ormsby,

  # forward-backward
  forward,
  viterbi,

  # likelihood approximations
  LikelihoodApproximation,
  Truncation,
  Projection
end


text:
function _clear_cache(network::AbstractGibbsNetwork, sol::Solution)
    i, j = node_from_index(network, length(first(sol.states))+1)
    if j != network.ncols return end
    delete!(memoize_cache(mps), (network, i))
    delete!(memoize_cache(dressed_mps), (network, i))
    delete!(memoize_cache(mpo), (network, i-1))
    delete!(memoize_cache(_mpo), (network, i-1))
    lec = memoize_cache(left_env)
    delete!.(Ref(lec), filter(k->k[2]==i, keys(lec)))
    rec = memoize_cache(right_env)
    delete!.(Ref(rec), filter(k->k[2]==i, keys(rec)))
end

export
    PEPSNetwork,
    generate_boundary,
    node_from_index,
    conditional_probability

struct PEPSNetwork <: AbstractGibbsNetwork{NTuple{2, Int}, NTuple{2, Int}}
    factor_graph::LabelledGraph{T, NTuple{2, Int}} where T
    network_graph::LabelledGraph{S, NTuple{2, Int}} where S
    vertex_map::Function
    m::Int
    n::Int
    nrows::Int
    ncols::Int
    β::Real
    bond_dim::Int
    var_tol::Real
    sweeps::Int
    gauges
    tensor_spiecies
    columns_MPO::NTuple{N, Union{Rational{Int}, Int}} where N
    layers_MPS::NTuple{M, Union{Rational{Int}, Int}} where M
    layers_left_env::NTuple{K, Union{Rational{Int}, Int}} where K
    layers_right_env::NTuple{L, Union{Rational{Int}, Int}} where L

    function PEPSNetwork(
        m::Int,
        n::Int,
        factor_graph::LabelledGraph,
        transformation::LatticeTransformation;
        β::Real,
        bond_dim::Int=typemax(Int),
        var_tol::Real=1E-8,
        sweeps::Int=4,
        columns_MPO = (-1//2, 0),  # from left to right
        # with gauges
        # layers_MPS=(1//6, 0, -3//6, -4//6),  # from bottom to top
        # layers_left_env=(1//6,),
        # layers_right_env=(0, -3//6)
        # with gauges but differently
        layers_MPS=(4//6, 3//6, 0, -1//6),  # from bottom to top
        layers_left_env=(4//6, 3//6),
        layers_right_env=(0, -3//6)
    )
        vmap = vertex_map(transformation, m, n)
        ng = peps_lattice(m, n)
        nrows, ncols = transformation.flips_dimensions ? (n, m) : (m, n)

        if !is_compatible(factor_graph, ng)
            throw(ArgumentError("Factor graph not compatible with given network."))
        end

        _types = (:site, :central_h, :central_v, :gauge_h)
        _gauges = Dict{Tuple{Rational, Int}, Vector{Float64}}()
        _tensor_spiecies = Dict()

        network = new(factor_graph, ng, vmap, m, n, nrows, ncols, β, bond_dim,
                      var_tol, sweeps, _gauges, _tensor_spiecies,
                      columns_MPO, layers_MPS, layers_left_env, layers_right_env
                )
        for type ∈ _types
            push!(network.tensor_spiecies, tensor_assignment(network, type)...)
        end

        update_gauges!(network, :id)
        network
    end
end

function projectors(network::PEPSNetwork, vertex::NTuple{2, Int})
    i, j = vertex
    neighbours = ((i, j-1), (i-1, j), (i, j+1), (i+1, j))
    projector.(Ref(network), Ref(vertex), neighbours)
end

node_index(peps::AbstractGibbsNetwork, node::NTuple{2, Int}) = peps.ncols * (node[1] - 1) + node[2]

node_from_index(peps::AbstractGibbsNetwork, index::Int) =
    ((index-1) ÷ peps.ncols + 1, mod1(index, peps.ncols))

function boundary(peps::PEPSNetwork, node::NTuple{2, Int})
    i, j = node
    x = (-4, -2)
    vcat(
        [
            [(x, x), ((i, k), (i+1, k))] for k ∈ 1:j-1
        ]...,
            (x, x), ((i, j-1), (i, j)),
        [
            [(x, x), ((i-1, k), (i, k))] for k ∈ j:peps.ncols
        ]...
    )
end

function conditional_probability(peps::PEPSNetwork, w::Vector{Int})
    i, j = node_from_index(peps, length(w)+1)
    ∂v = boundary_state(peps, w, (i, j))

    L = left_env(peps, i, ∂v[1:2*j-1])
    R = right_env(peps, i, ∂v[2*j+3 : 2*peps.ncols+2])
    A = reduced_site_tensor(peps, (i, j), ∂v[2*j], ∂v[2*j+2])

    ψ = dressed_mps(peps, i)
    M = ψ[2 * j]

    @tensor prob[σ] := L[x] * M[x, d, y] * A[r, d, σ] *
                       R[y, r] order = (x, d, r, y)

    normalize_probability(prob)
end

function bond_energy(
    network::AbstractGibbsNetwork,
    u::NTuple{2, Int},
    v::NTuple{2, Int},
    σ::Int
)
    fg_u, fg_v = network.vertex_map(u), network.vertex_map(v)
    if has_edge(network.factor_graph, fg_u, fg_v)
        pu, en, pv = get_prop.(Ref(network.factor_graph), Ref(fg_u), Ref(fg_v), (:pl, :en, :pr))
        energies = (pu * (en * pv[:, σ:σ]))'
    elseif has_edge(network.factor_graph, fg_v, fg_u)
        pv, en, pu = get_prop.(Ref(network.factor_graph), Ref(fg_v), Ref(fg_u), (:pl, :en, :pr))
        energies = (pv[σ:σ, :] * en) * pu
    else
        energies = zeros(length(local_energy(network, u)))
    end
    vec(energies)
end

function update_energy(network::PEPSNetwork, σ::Vector{Int})
    i, j = node_from_index(network, length(σ)+1)
    bond_energy(network, (i, j), (i, j-1), local_state_for_node(network, σ, (i, j-1))) +
    bond_energy(network, (i, j), (i-1, j), local_state_for_node(network, σ, (i-1, j))) +
    local_energy(network, (i, j))
end

export FusedNetwork, boundary,
    conditional_probability,
    update_energy,
    projectors

struct FusedNetwork <: AbstractGibbsNetwork{NTuple{2, Int}, NTuple{2, Int}}
    factor_graph::LabelledGraph{T, NTuple{2, Int}} where T
    network_graph::LabelledGraph{S, NTuple{2, Int}} where S
    vertex_map::Function
    m::Int
    n::Int
    nrows::Int
    ncols::Int
    β::Real
    bond_dim::Int
    var_tol::Real
    sweeps::Int
    gauges
    tensor_spiecies
    columns_MPO::NTuple{N, Union{Rational{Int}, Int}} where N
    layers_MPS::NTuple{M, Union{Rational{Int}, Int}} where M
    layers_left_env::NTuple{K, Union{Rational{Int}, Int}} where K
    layers_right_env::NTuple{L, Union{Rational{Int}, Int}} where L

    function FusedNetwork(
        m::Int,
        n::Int,
        factor_graph::LabelledGraph,
        transformation::LatticeTransformation;
        β::Real,
        bond_dim::Int=typemax(Int),
        var_tol::Real=1E-8,
        sweeps::Int=4,
        columns_MPO = (-1//2, 0),  # from left to right
        # with gauges
        # layers_MPS=(1//6, 0, -3//6, -4//6),  # from bottom to top
        # layers_left_env=(1//6,),
        # layers_right_env=(0, -3//6)
        # with gauges
         layers_MPS=(4//6, 3//6, 0, -1//6),  # from bottom to top
         layers_left_env=(4//6, 3//6),
         layers_right_env=(0, -3//6)
        # without gauges
        # layers_MPS=(0, -3//6),  # from bottom to top
        # layers_left_env=(),
        # layers_right_env=(0, -3//6)
        # no gauges
        # layers_MPS=(3//6, 0),  # from bottom to top
        # layers_left_env=(3//6,),
        # layers_right_env=(0, -3//6)
    )
    vmap = vertex_map(transformation, m, n)
    ng = cross_lattice(m, n)
    nrows, ncols = transformation.flips_dimensions ? (n, m) : (m, n)

    if !is_compatible(factor_graph, ng)
        throw(ArgumentError("Factor graph not compatible with given network."))
    end

    _types = (:site, :central_h, :central_v, :virtual, :central_d, :gauge_h)
    _gauges = Dict()

    _tensor_spiecies = Dict()

    network = new(factor_graph, ng, vmap, m, n, nrows, ncols, β, bond_dim,
                  var_tol, sweeps, _gauges, _tensor_spiecies,
                  columns_MPO, layers_MPS, layers_left_env, layers_right_env
            )

    for type ∈ _types
        push!(network.tensor_spiecies, tensor_assignment(network, type)...)
    end

    update_gauges!(network, :id)
    network
    end
end

function projectors(network::FusedNetwork, vertex::NTuple{2, Int})
    i, j = vertex
    neighbours = (
                    (
                        (i+1, j-1), (i, j-1), (i-1, j-1)
                    ),
                    (i-1, j),
                    (
                        (i+1, j+1), (i, j+1), (i-1, j+1)
                    ),
                    (i+1, j)
                )
    projector.(Ref(network), Ref(vertex), neighbours)
end

function boundary(peps::FusedNetwork, node::NTuple{2, Int})
    i, j = node
    vcat(
        [
            [((i, k-1), (i+1, k), (i, k), (i+1, k-1)), ((i, k), (i+1, k))] for k ∈ 1:j-1
        ]...,
        (
            (i, j-1), (i+1, j)
        ),
        (
            (i, j-1), (i, j)
        ),
        (
            (i-1, j-1), (i, j)
        ),
        (
            (i-1, j), (i, j)
        ),
        [
            [((i-1, k-1), (i, k), (i-1, k), (i, k-1)), ((i-1, k), (i, k))] for k ∈ j+1:peps.ncols
        ]...
    )
end

function conditional_probability(peps::FusedNetwork, w::Vector{Int})
    i, j = node_from_index(peps, length(w)+1)
    ∂v = boundary_state(peps, w, (i, j))

    L = left_env(peps, i, ∂v[1:2*j-2])
    R = right_env(peps, i, ∂v[2*j+3 : 2*peps.ncols+2])
    A = reduced_site_tensor(peps, (i, j), ∂v[2*j-1], ∂v[2*j], ∂v[2*j+1], ∂v[2*j+2])

    ψ = dressed_mps(peps, i)
    MX, M = ψ[2*j-1], ψ[2*j]

    @tensor prob[σ] := L[x] * MX[x, m, y] * M[y, l, z] * R[z, k] *
                        A[k, l, m, σ] order = (x, y, z, k, l, m)

    normalize_probability(prob)
end

function update_energy(network::FusedNetwork, σ::Vector{Int})
    i, j = node_from_index(network, length(σ)+1)
    bond_energy(network, (i, j), (i, j-1), local_state_for_node(network, σ, (i, j-1))) +
    bond_energy(network, (i, j), (i-1, j), local_state_for_node(network, σ, (i-1, j))) +
    bond_energy(network, (i, j), (i-1, j-1), local_state_for_node(network, σ, (i-1, j-1))) +
    bond_energy(network, (i, j), (i-1, j+1), local_state_for_node(network, σ, (i-1, j+1))) +
    local_energy(network, (i, j))
end


text:
#=
Cannot reproduce file dialog examples #189 - Closed
rowlesmr opened this issue on Oct 1, 2015 · 6 comments
 Closed
Cannot reproduce file dialog examples
#189
rowlesmr opened this issue on Oct 1, 2015 · 6 comments
Comments
@rowlesmr rowlesmr commented on Oct 1, 2015
I have tried these in Julia v0.3.11 & Gtk 0.8.6 and Julia v0.4.0-rc3 & Gtk 0.9.2

I have copied the code from the readme.md file under the section "File dialogs" right at the end of the file.

I am running it in REPL.I have entered using Gtk and using Gtk.ShortNames.
=#

open_dialog("Pick some files", multiple=true) #gives a file dialog, but also gives a warning: GLib-GObject-WARNING **: g_object_set_property: object class 'GtkFileChooseDialog' has no property 'multiple'. I am not able to choose multiple files.

open_dialog("Pick a file", filters=("*.jl",)) #gives an error: ERROR: 'g_type' has no method matching g_type(::Type{(ASCIIString,)}). No dialog is opened. Removing the final , removes this error and opens a dialog box, but a warning is given: GLib-GObject-WARNING **: g_object_set_property: object class 'GtkFileChooseDialog' has no property 'filters'. The list of files is not filtered.

open_dialog("Pick a file", filters=(@FileFilter(mimetype="text/csv"),)) #gives an error: ERROR: @FileFilter not defined.

#=
What am I doing wrong?

Contributor
@vtjnash vtjnash commented on Oct 1, 2015
what OS? what version of Gtk/GLib are you using?

Author
@rowlesmr rowlesmr commented on Oct 1, 2015
Windows 7, service pack 1, 64 bit
Julia v0.3.11 & Gtk 0.8.6 and Julia v0.4.0-rc3 & Gtk 0.9.2

Gtk was installed by Pkg.add("Gtk")

Contributor
@lobingera lobingera commented on Oct 2, 2015
I can reproduce all mentioned errors/warnings with
=#
versioninfo()
#Julia Version 0.4.0-rc3
#=Commit 483d548* (2015-09-27 20:34 UTC)
Platform Info:
  System: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i3-2120 CPU @ 3.30GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

lobi@orange4:~/juliarepo$ ../julia04/julia
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.4.0-rc3 (2015-09-27 20:34 UTC)
 _/ |\__'_|_|_|\__'_|  |
|__/                   |  x86_64-linux-gnu
=#
using Gtk

using Gtk.ShortNames # (non-)Official use of ShortNames

open_dialog("Pick some files", multiple=true)
#=WARNING: Base.Uint8 is deprecated, use UInt8 instead.

(julia:9963): GLib-GObject-WARNING **: g_object_set_property: object class 'GtkFileChooserDialog' has no property named 'multiple'
"""=#

julia> open_dialog("Pick a file", filters=("*.jl",))

#=ERROR: MethodError: `g_type` has no method matching g_type(::Type{Tuple{ASCIIString}})
 in setindex! at /home/lobi/.julia/v0.4/Gtk/src/GLib/gvalues.jl:64
 in gvalue at /home/lobi/.julia/v0.4/Gtk/src/GLib/gvalues.jl:18
 in setproperty! at /home/lobi/.julia/v0.4/Gtk/src/GLib/gvalues.jl:163
 in GtkFileChooserDialogLeaf at /home/lobi/.julia/v0.4/Gtk/src/GLib/gtype.jl:209
 in GtkFileChooserDialogLeaf at /home/lobi/.julia/v0.4/Gtk/src/selectors.jl:28
 in open_dialog at /home/lobi/.julia/v0.4/Gtk/src/selectors.jl:73 (repeats 2 times)
=#

open_dialog("Pick a file", filters=("*.jl"))

#=
(julia:9963): GLib-GObject-WARNING **: g_object_set_property: object class 'GtkFileChooserDialog' has no property named 'filters'
""
=#
open_dialog("Pick a file", filters=(@FileFilter(mimetype="text/csv"),))
#=
ERROR: UndefVarError: @FileFilter not defined
=#
@FileFilter(; name = nothing, pattern = "", mimetype = "")
#=
ERROR: UndefVarError: @FileFilter not defined
@silvernode silvernode commented on Oct 2, 2015
On Arch Linux, Gtk loads (using Gtk) but I have all sorts of errors. @window is undefined, as well as @frame. getting rid of the '@' symbol fixes window but not Frame. Frame is still undefined. I haven't been able to get the examples to work on any operating system. I don't know whether I am missing something or Gtk just does not work with Julia.

Contributor
@lobingera lobingera commented on Oct 2, 2015
@silvernode or Gtk just does not work with Julia, i'm using Gtk.jl for quite some time right now and apart from some dependency hicups and general julia restructuring, i'm fine. I think part of the issues above link to julia restructuring on tuples and e.g. String->AbstractString.

@vtjnash vtjnash closed this in JuliaGraphics/Gtk.jl@3f31e5e on Oct 2, 2015
Contributor
@vtjnash vtjnash commented on Oct 2, 2015
@silvernode you need using Gtk.ShortNames to use the abbreviations, otherwise everything uses the full name from the Gtk docs.
=#

#Solution:
#@rowlesmr I've fixed the examples. I have reserved keyword arguments in Gtk.jl to be used for the Gtk properties (e.g. https://developer.gnome.org/gtk3/stable/GtkFileChooser.html#GtkFileChooser--select-multiple), so the examples needed to be corrected to reflect this.
#me:links leads to no where

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
"""
Compute the dark matter density profile at the point (x, vr, vt).  Where x is the radius divided by the scale radius (r / rs), vr is the radial velocity, and vt is the tangential velocity.
"""
function profile_density(x::G, vr::G, vt::G, p::NFWParameters{G}; N::G = 1000.0, αt::G = 1.0) where G <: AbstractFloat
    if x < 0.0
        return 0.0
    end

    #Compute Energy
    E = 0.5 * (vt * vt + vr * vr) + gravitational_potential(x, p)
    if E < p.Φlim && E >= 0.0
        hE = N * E^p.a * (E^p.q + p.adjEc^p.q)^(p.d/p.q) * (p.Φlim - E)^p.e
    else
        return 0.0
    end

    #Compute angular momentum
    J = abs(x * p.rs * vt) #vt = v * sin(theta)
    if p.b <= 0.0
        #gJ = 1.0 / (1.0 + (J / p.adjJb)^(-p.b))
        gJ = (1.0 + (J / p.adjJb)^(-p.b / αt))^-αt
    else
        #gJ = 1.0 + (J / p.adjJb)^p.b
        gJ = (1.0 + (J / p.adjJb)^(p.b / αt))^αt
    end

    return hE * gJ
end

function log_profile_density(x::G, vr::G, vt::G, p::NFWParameters{G}; N::G = 1000.0, αt::G = 1.0) where G <: AbstractFloat
    if x < 0.0
        return -Inf
    end

    #Compute Energy
    E = 0.5 * (vt * vt + vr * vr) + gravitational_potential(x, p)
    lE = log(E)
    if E < p.Φlim && E >= 0.0
        lhE = log(N) + p.a * lE + (p.d/p.q) * log(E^p.q + p.adjEc^p.q) + p.e * log(p.Φlim - E)
    else
        return -Inf
    end

    #Compute angular momentum
    J = abs(x * p.rs * vt) #vt = v * sin(theta)
    if p.b <= 0.0
        #gJ = 1.0 / (1.0 + (J / p.adjJb)^(-p.b))
        lgJ = -αt * log(1.0 + (J / p.adjJb)^(-p.b / αt))
    else
        #gJ = 1.0 + (J / p.adjJb)^p.b
        lgJ = αt * log(1.0 + (J / p.adjJb)^(p.b / αt))
    end

    return lhE + lgJ
end

function profile_density(x::G, vr::G, vt::G, p::SFWParameters{G}; N::G = 1000.0, αt::G = 1.0) where G <: AbstractFloat
    if x < 0.0
        return 0.0
    end

    #Compute Energy
    E = 0.5 * (vt * vt + vr * vr) + gravitational_potential(x, p)
    if E < p.Φlim && E >= 0.0
        hE = N * E^p.a * (E^p.q + p.adjEc^p.q)^(p.d/p.q) * (p.Φlim - E)^p.e
    else
        return 0.0
    end

    #Compute angular momentum
    J = abs(x * p.rs * vt) #vt = v * sin(theta)
    if p.b <= 0.0
        #gJ = 1.0 / (1.0 + (J / p.adjJb)^(-p.b))
        gJ = (1.0 + (J / p.adjJb)^(-p.b / αt))^-αt
    else
        #gJ = 1.0 + (J / p.adjJb)^p.b
        gJ = (1.0 + (J / p.adjJb)^(p.b / αt))^αt
    end

    return hE * gJ
end

function log_profile_density(x::G, vr::G, vt::G, p::SFWParameters{G}; N::G = 1000.0, αt::G = 1.0) where G <: AbstractFloat
    if x < 0.0
        return -Inf
    end

    #Compute Energy
    E = 0.5 * (vt * vt + vr * vr) + gravitational_potential(x, p)
    lE = log(E)
    if E < p.Φlim && E >= 0.0
        lhE = log(N) + p.a * lE + (p.d/p.q) * log(E^p.q + p.adjEc^p.q) + p.e * log(p.Φlim - E)
    else
        return -Inf
    end

    #Compute angular momentum
    J = abs(x * p.rs * vt) #vt = v * sin(theta)
    if p.b <= 0.0
        lgJ = -αt * log(1.0 + (J / p.adjJb)^(-p.b / αt))
    else
        lgJ = αt * log(1.0 + (J / p.adjJb)^(-p.b / αt))
    end

    return lhE + lgJ
end


text:
/**
 * This class models a cloud storage API.
 */
class Cloud {
    int{L} cloud;

    /**
     * Put a value into the cloud.
     */
    void put(int x) {
	this.cloud = x;
    }
    /**
     * Put the only value the cloud stores.
     */
    int get() {
	return this.cloud;
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
# diffusion variance

abstract type AbstractDiffusion end

mutable struct UnitDiffusion <: AbstractDiffusion end


text:
"""
# Linear programming and reductions

- An introduction to linear programming
    - Example: profit maximization
    - Solving linear programs
    - More products
    - A magic trick called duality
    - Example: production planning
    - Example: optimum bandwidth allocation
    - Reductions
    - Variants of linear programming
    - Matrix-vector notation
- Flows in networks
    - Shipping oil
    - Maximizing flow
    - A closer look at the algorithm
    - A certificate of optimality
    - Max-flow min-cut theorem
    - Efficiency
- Bipartite matching
- Duality
    - Duality theorem
    - Visualizing duality
- Zero-sum games
- The simplex algorithm
    - Vertices and neighbors in n-dimensional space
    - The algorithm
    - Loose ends
    - The starting vertex
    - Degeneracy
    - Unboundedness
    - The running time of simplex
    - Gaussian elimination
    - Linear programming in polynomial time
- Postscript: circuit evaluation
"""
module Chapter7

end # module


text:
#include all of the test files.

@time begin include("test_deviations.jl") end

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
module ConvexTests

export do_tests

using TableTestSets
using Convex
using Convex.ProblemDepot: foreach_problem, PROBLEMS
using Test, Pkg, InteractiveUtils, Dates
using TimerOutputs

const DOCS_SRC = Ref(joinpath(@__DIR__, "..", "docs", "src"))

# To compile some things
function dummy_problem(opt; T = Float64)
    x = Variable(3)
    p = minimize(2*x[2] + x[1] - x[3], [x >= 1]; numeric_type = T)
    solve!(p, opt)
end

# Modified from Convex.jl's ProblemDepot code to add TimerOutputs logging
function _run_tests( handle_problem!::Function,
                    problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing;
                    exclude::Vector{Regex} = Regex[], T=Float64, atol=1e-3, rtol=0.0, to)
    exclude = copy(exclude) # rebind the local variable
    push!(exclude, r"benchmark")
    for class in keys(PROBLEMS)
        any(occursin.(exclude, Ref(class))) && continue
        @timeit to "$class" begin
            @testset "$class" begin
                foreach_problem(class, problems; exclude=exclude) do name, problem_func
                    @timeit to "$name" begin
                        @testset "$name" begin
                            problem_func(handle_problem!, Val(true), atol, rtol, T)
                        end
                    end
                end
            end
        end
    end
end

function formatted_seconds(t)
    if t < 1
        input = Dates.Millisecond(ceil(Int, t*1000))
    else
        input = Dates.Second(ceil(Int, t))
    end
    Dates.CompoundPeriod(input) |> canonicalize
end

function do_tests(name, opt; variant="", first = true, last = true, description = "",T=Float64, exclude = Regex[], kwargs...)
    if first
        t1 = @elapsed dummy_problem(opt; T=T)
        t2 = @elapsed dummy_problem(opt; T=T)
        compilation_time = formatted_seconds(t1-t2)
    end

    to = TimerOutput()
    results, t = @timed begin
        @testset TableTestSet "$name tests" begin
            _run_tests(; to=to, T=T, exclude=exclude, kwargs...) do p
                solve!(p, opt)
            end
        end
    end
    duration = formatted_seconds(t)
    filename = joinpath(DOCS_SRC[], "$(name).md")

    open(filename, write=true, append=!first) do io
        if first
            println(io, """
            Table of contents:

            ```@contents
            Pages = ["$(name).md"]
            Depth = 4
            ```
            """)
            println(io)
            println(io, "Compilation warmup gives an estimate of $(compilation_time) of compilation time.")
            println(io)
        end

        println(io, "## $name $variant")

        datestr = Dates.format(Dates.now(Dates.UTC), "U d, Y at HH:MM")

        println(io, "These tests were run on $(datestr) (UTC).")
        println(io)
        if !isempty(description)
            println(io, description)
        end
        println(io)
        if !isempty(exclude)
            println(io, "Excluded problems and classes of problems:")
            println(io, "```julia")
            println(io, exclude)
            println(io, "```")
        else
            println(io, "No problems excluded.")
        end
        println(io)
        println(io, "### Tests")
        println(io)
        println(io, "Tests took $(duration) to run (after warmup).")
        println(io)
        println(io, "```@raw html")
        html_table(io, results; standalone = false)
        println(io, "```")
        println(io)
        println(io, "### Errors")
        println(io)
        println(io, "```julia")
        TableTestSets.print_test_errors(io, results)
        println(io, "```")
        println(io)
        println(io)
        println(io, "### Timing information")
        println(io, "```julia")
        print_timer(io, to)
        println(io)
        println(io, "```")
        println(io)
        if last
            println(io, "## Version information")
            println(io, "`versioninfo()`:")
            println(io, "```julia")
            versioninfo(io)
            println(io, "```")
            println(io)
            println(io, "Manifest:")
            println(io, "```julia")
            redirect_stdout(io) do
                Pkg.status(;mode = Pkg.PKGMODE_MANIFEST)
            end
            println(io, "```")
        end

    end
    nothing
end

end # module