text:
# Shows are the same as podcasts on Spotify

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-a-show

"""
    show_get(show_id; market="US")

**Summary**: Get a Spotify catalog information for a single show identified by it's unique Spotify ID.

# Arguments
- `show_id` : The Spotify ID for the show

# Optional keywords
- `market` : An ISO 3166-1 alpha-2 country code. If a country code is specified, only shows
             and episodes that are available in that market will be returned.

# Example
```julia-repl
julia> Spotify.show_get_single("2MAi0BvDc6GTFvKFPXnkCL")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 18 entries:
  :available_markets    => ["AD", "AE", "AG", "AL", "AM", "AR", "AT", "AU", "BA", "BB"  …  "TV", "TW", "US", "UY", "VC", …
  :copyrights           => Union{}[]
  :description          => "Conversations about science, technology, history, philosophy and the nature of intelligence, …
  :episodes             => {…
```
"""
function show_get_single(show_id; market::String="US")
    return spotify_request("shows/$show_id?market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-multiple-shows

"""
    show_get_multiple(ids; market::String="US")

**Summary**: Get Spotify catalog information for several shows based on their Spotify IDs.

# Arguments
- `ids` : A comma-separated list of the Spotify IDs for the shows. Maximum: 50 IDs.

# Optional keywords
- `market` : An ISO 3166-1 alpha-2 country code. If a country code is specified, only shows
             and episodes that are available in that market will be returned.

# Example
```julia-repl
julia> Spotify.show_get_multiple("2MAi0BvDc6GTFvKFPXnkCL, 4rOoJ6Egrf8K2IrywzwOMk")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 1 entry:
  :shows => JSON3.Object[{…
```
"""
function show_get_multiple(ids; market::String="US")
    return spotify_request("shows/?ids=$ids&market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-a-shows-episodes
"""
    show_get_episodes(show_id; market="US", limit=20, offset=0)

**Summary**: Get Spotify catalog information about a show’s episodes. Optional parameters
             can be used to limit the number of episodes returned.

# Arguments
- `show_id` : The Spotify ID for the show

# Optional keywords
- `market::String` : An ISO 3166-1 alpha-2 country code. If a country code is specified,
                     only episodes that are available in that market will be returned.
                     Default is set to "US".
- `limit::Int64` : Maximum number of items to return, default is set to 20. (0 < limit <= 50)
- `offset::Int64` : Index of the first item to return, default is set to 0

# Example
```julia-repl
julia> Spotify.show_get_episodes("2MAi0BvDc6GTFvKFPXnkCL")[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 7 entries:
  :href     => "https://api.spotify.com/v1/shows/2MAi0BvDc6GTFvKFPXnkCL/episodes?offset=0&limit=20&market=US"
  :items    => JSON3.Object[{…
  :limit    => 20
```
"""
function show_get_episodes(show_id; market::String="US", limit::Int64=20, offset::Int64=0)
    return spotify_request("shows/$show_id/episodes?market=$market&limit=$limit&offset=$offset")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-users-saved-shows

"""
    show_get_saved(;limit::Int64=20, offset::Int64=0)

**Summary**: Get a list of shows saved in the current Spotify user's library. Optional parameters can
             be used to limit the number of shows returned.

# Optional keywords
- `limit::Int64` : Maximum number of items to return, default is set to 20. (0 < limit <= 50)
- `offset::Int64` : Index of the first item to return, default is set to 0

# Example
```julia-repl
julia> Spotify.show_get_saved()[1]
[ Info: We try the request without checking if current grant includes scope user-library-read.
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 7 entries:
  :href     => "https://api.spotify.com/v1/me/shows?offset=0&limit=20"
  :items    => Union{}[]
  :limit    => 20
```
"""
function show_get_saved(;limit::Int64=20, offset::Int64=0)
    return spotify_request("me/shows?limit=$limit&offset=$offset"; scope = "user-library-read")
end

## https://developer.spotify.com/documentation/web-api/reference/#/operations/check-users-saved-shows

"""
    show_get_contains(ids)

**Summary**: Check if one or more shows is already saved in the current Spotify user's library.

# Arguments
- `ids` : A comma-separated list of the Spotify IDs for the shows. Maximum: 50 IDs.
"""
function show_get_contains(ids)
    return spotify_request("me/shows/contains?ids=$ids"; scope = "user-library-read")
end


text:
function findSurface(alpha,beta,M,a,b,k,l;camera::NTuple{2,Float64}=(-30, 30), plotSlice::Bool=false, saveFig::Bool=false)

    const a_prime = a*cosd(alpha)
    const k_prime = k*cosd(beta)

end

function

ext: jl


text:
using ModelObjectsLH, ModelParams, Test

function make_bounded_test_vector(vLength :: Integer, increasing :: Symbol)
    ownId = ObjectId(:test1);
    lb = 3.0;
    ub = 7.0;
    if increasing == :nonmonotone
        if vLength > 1
            dxV = collect(LinRange(lb + 0.1, ub - 0.1, vLength));
        else
            dxV = [lb + 1.5];
        end
    else
        if vLength > 1
            dxV = collect(range(0.4, 0.2, length = vLength));
        else
            dxV = [0.4];
        end
    end
    iv = BoundedVector(ownId, ParamVector(ownId), increasing, lb, ub, dxV);
    set_pvector!(iv);
    return iv
end

@testset "BoundedVector" begin
    for vLength in [1, 3]
        for increasing in [:increasing, :decreasing, :nonmonotone]
            iv = make_bounded_test_vector(vLength, increasing);
            n = length(iv);
            xV = values(iv);
            @test ModelParams.check_values(iv, xV)

            for j = 1 : n
                @test values(iv, j) == xV[j]
            end

            xNewV = xV .+ 0.05;
    		fix_values!(iv, xNewV);
            x2V = values(iv);
            @test isapprox(xNewV, x2V)
        end
    end
end

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


text:
"""
$(TYPEDSIGNATURES)

Wrap a [`AnyElement`](@ref). Use until specialized/cooked.

From the 'primer': built-in sorts of Symmetric Nets are the following:
  Booleans, range of integers, finite enumerations, cyclic enumerations and dots
"""
struct BuiltInSort <: AbstractSort
    dict::AnyElement
end

"""
$(TYPEDSIGNATURES)

Wrap a [`AnyElement`](@ref). Use until specialized/cooked.
"""
struct MultisetSort <: AbstractSort
    dict::AnyElement
end

"""
$(TYPEDSIGNATURES)

Wrap a [`AnyElement`](@ref). Use until specialized/cooked.
Should contain an ordered collection of sorts.
"""
struct ProductSort <: AbstractSort
    dict::AnyElement
end

"""
$(TYPEDSIGNATURES)

Wrap a [`AnyElement`](@ref). Use until specialized/cooked.
"""
struct UserSort <: AbstractSort
    dict::AnyElement
end

"""
Return default  based on `PNTD`. Has meaning of empty, as in `zero`.
"""
function default_sort end
default_sort(::PNTD) where {PNTD <: PnmlType} = zero(Integer) #!
default_sort(::PNTD) where {PNTD <: AbstractContinuousCore} = zero(Float64) #!
default_sort(::PNTD) where {PNTD <: AbstractHLCore} = Sort() #!
# For a HL Net the sort is Dot.

"""
$(TYPEDEF)
$(TYPEDFIELDS)

Part of the many-sorted algebra attached to nodes on a Petri Net Graph.
Is content of a <structure> element of a High-Level label.
"""
struct Sort #TODO
    dict::PnmlDict #TODO AnyElement for bring-up? What should be here?
    #TODO xml
end

Sort() = Sort(PnmlDict())
convert(::Type{Maybe{Sort}}, pdict::PnmlDict) = Sort(pdict)


text:
@testset "@loadOrGenerate" begin
    using YAXArrayBase
    using Zarr

    data_a = collect(reshape(rand(20), 4, 5))
    axlist = [RangeAxis("XVals", 1.0:4.0), CategoricalAxis("YVals", [1, 2, 3, 4, 5])]
    props = Dict("att1" => 5, "att2" => "Hallo")

    @loadOrGenerate a => "a.zarr" begin
        a = YAXArray(axlist, data_a, props)
    end

    @test isdir("a.zarr")

    data_a2 = collect(reshape(rand(20), 4, 5))
    @loadOrGenerate a2 => "a.zarr" begin
        a2 = YAXArray(axlist, data_a2, props)
    end

    @test all(a.data .=== a2[:, :])

    data_b = collect(reshape(rand(20), 4, 5))
    data_c = collect(reshape(rand(20), 4, 5))
    @loadOrGenerate b => "b.zarr" c => "c.zarr" begin
        b = YAXArray(axlist, data_b, props)
        c = YAXArray(axlist, data_c, props)
    end

    data_b2 = collect(reshape(rand(20), 4, 5))
    data_c2 = collect(reshape(rand(20), 4, 5))
    @loadOrGenerate b2 => "b.zarr" c2 => "c.zarr" begin
        b2 = YAXArray(axlist, data_b2, props)
        c2 = YAXArray(axlist, data_c2, props)
    end

    @test all(b.data .=== b2[:, :])
    @test all(c.data .=== c2[:, :])

    d_dir = "d.zarr"
    data_d = collect(reshape(rand(20), 4, 5))
    @loadOrGenerate d => d_dir begin
        d = YAXArray(axlist, data_d, props)
    end

    @test isdir(d_dir)

    data_d2 = collect(reshape(rand(20), 4, 5))
    @loadOrGenerate d2 => d_dir begin
        d2 = YAXArray(axlist, data_d2, props)
    end

    @test all(d.data .=== d2[:, :])

    rm("a.zarr", recursive = true)
    rm("b.zarr", recursive = true)
    rm("c.zarr", recursive = true)
    rm("d.zarr", recursive = true)
end

max_stars_repo_name: meggart/YAXArrays

max_forks_repo_name: meggart/YAXArrays


text:
function laszip_error(laszip_obj::Ptr{Cvoid})
    errstr = Ref(Cstring(C_NULL))
    laszip_get_error(laszip_obj, errstr)
    if errstr[] != C_NULL
        error(unsafe_string(errstr[]))
    end
    nothing
end

macro check(obj, ex)
    return :( $(esc(ex)) == 0 ? nothing : laszip_error($(esc(obj))) )
end

function readstring(s::NTuple{T,UInt8}) where T
    lastchar = findlast(s .!= 0)
    if isnothing(lastchar)
        return ""
    else
        return join(Char.(s[1:lastchar]))
    end
end

function writestring(str::AbstractString, nb::Integer)
    s = ascii(str)
    n = length(s)
    npad = nb - n
    if npad < 0
        error("string too long")
    elseif npad == 0
        NTuple{nb,UInt8}(codeunits(ascii(s)))
    else
        writestr = string(s * "\0"^npad)
        NTuple{nb,UInt8}(codeunits(ascii(writestr)))
    end
end

max_stars_repo_name: evetion/LazIO

max_issues_repo_name: evetion/LazIO

max_forks_repo_name: evetion/LazIO


text:
module Tmp
export say_hello

say_hello()=println("Hi")

# end of module
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
export SnfMcmcInsertDelete, SnfMcmcFlip, SnfMcmcSysGibbs, SnfMcmcRandGibbs
export SnfMcmcSampler

abstract type SnfMcmcSampler end

struct SnfMcmcInsertDelete <: SnfMcmcSampler
    ν::Int
    nz_edges::Vector{Int}
    desired_samples::Int
    burn_in::Int
    lag::Int
end

function SnfMcmcInsertDelete(
    ;
    ν::Int=3,
    desired_samples::Int=1000,
    burn_in::Int=0,
    lag::Int=1
    )
    return SnfMcmcInsertDelete(ν, Int[], desired_samples, burn_in, lag)
end

struct SnfMcmcFlip <: SnfMcmcSampler
    ν::Int
    desired_samples::Int
    burn_in::Int
    lag::Int
end

struct SnfMcmcSysGibbs <: SnfMcmcSampler
    ν::Int
    desired_samples::Int
    burn_in::Int
    lag::Int
end

function SnfMcmcSysGibbs(
    ;
    ν::Int=3,
    desired_samples::Int=1000,
    burn_in::Int=0,
    lag::Int=1
    )
    return SnfMcmcSysGibbs(ν, desired_samples, burn_in, lag)
end

struct SnfMcmcRandGibbs <: SnfMcmcSampler
    ν::Int  # Neighbourhood size parameter for component updates
    desired_samples::Int
    burn_in::Int
    lag::Int
end

function SnfMcmcRandGibbs(
    ;
    ν::Int=3,
    desired_samples::Int=1000,
    burn_in::Int=0,
    lag::Int=1
    )
    return SnfMcmcRandGibbs(ν, desired_samples, burn_in, lag)
end

max_stars_repo_path: src/Graph Models/SNF/SNF_samplers.jl

max_issues_repo_path: src/Graph Models/SNF/SNF_samplers.jl

max_forks_repo_path: src/Graph Models/SNF/SNF_samplers.jl


text:
function completion(string)
    REPLCompletions.completions(string, length(string))[1]
end

function incomplete(string)
    e = Base.parse_input_line(string)
    isa(e,Expr) && e.head === :incomplete
end

max_stars_repo_name: KevCaz/JuliaCall

max_issues_repo_name: tinagui/JuliaCall

max_forks_repo_name: tinagui/JuliaCall


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
import JuliaAcademyData; JuliaAcademyData.activate("Foundations of machine learning")
using JuliaAcademyData;
# ## Functions
#
#
# In the last notebook, we talked about modeling data with functions. A **function** is one of the most fundamental concepts in computing (and also in mathematics).
#
# A function is a piece of a program that receives **input arguments**, processes them by doing certain calculations on them, and returns **outputs**.
#
# For example, we might have a function `g` that takes a number as an input and returns the square of that number as an output. How can we define this function `g` on a computer? Julia gives us a few different ways to do this.

#-

# ### Defining functions

#-

# Firstly, we could write `g` as follows:

g(x) = x^2

#-

g(2), g(3.5)

#-

g(3.0)

# Alternatively, we could declare this function using the `function` and `end` keywords:

function g1(x)
    return x^2
end

#-

g1(2), g1(3.5)

# The third way we could have declared this function is as an "anonymous" or "lambda" function. "Anonymous" functions are functions that truly don't need names! For example, we could have declared a function that squares its input as

(x->x^2)(2)

# Now that we've done that, we can't access the function `x -> x^2` again because we have no name to call! That seems a little silly, doesn't it?
#
# Actually, there are times where functions without names are especially handy.  Most commonly they're used arguments to "higher-order" functions: these are functions which take _other functions_ as arguments. For example:

map(sqrt, [1, 2, 3]) # equivalent to [sqrt(1), sqrt(2), sqrt(3)]

#-

map(x->x^2, [1, 2, 3])

# Of course, you can also assign anonymous function a name, but in general it's better style to use the non-anonymous equivalents above if you actually want to use a name.

g2 = x -> x^2
g2(3.5), g2("I ♡ Julia ") # Use \heartsuit + TAB to get the ♡ character

# ## An important sigmoidal function

#-

# A particular function that is used a lot in machine learning is a so-called "sigmoidal" function (meaning a function that is S-shaped, i.e. the graph of the function looks like an `S`).
#
# The sigmoid function that we will use is given the name $\sigma$, and is defined by the following mathematical expression:
#
# $$\sigma(x) := \frac{1}{1 + \exp(-x)}.$$

#-

# #### Exercise 1
#
# Use the first syntax given above to define the function `σ` in Julia. Note that Julia actually allows us to use the symbol σ as a variable name! To do so, type `\sigma<TAB>` in the code cell.

#-

# #### Solution

σ(x) = 1 / (1 + exp(-x)) # \sigma + <tab>

#-

σ(1)

#-

function σ1(x)
    y = 1/(1+ℯ^(-x))
    return y
end

#-

σ1(1)

#-

(x->1/(1+ℯ^(-x)))(1)

# ## Plotting functions

#-

# Let's draw the function σ to see what it looks like. Throughout this course, we'll use the Julia package `Plots.jl` for all of the graphics. This package provides a flexible syntax for plotting, in which options to change attributes like the width of the lines used in the figure are given as named keyword arguments.
#
# In addition, it allows us to use different "backends", which are the other libraries that actually carry out the plotting following the instructions from `Plots.jl`.

using Plots
gr()   # use the GR "backend" (plotting library)

#-

plot(σ, -5, 5)

hline!([0, 1], ls=:dash, lw=3)  # add horizontal lines at 0 and 1, with dashed style and linewidth 3
vline!([0], ls=:dash, lw=3)     # add a vertical line at 0

# We can think of $\sigma$ as a smooth version of a step or threshold function (often called a "Heaviside" function). To see this, let's modify the steepness of the jump in $\sigma$ and compare it to the Heaviside function; we'll see how all this works in more detail later:

#-

heaviside(x) = x < 0 ? 0.0 : 1.0

#-

## try manipulating the value of `w` between 0 to :
w = 10.0
plot(x -> σ(w*x), -5, 5, label="sigma", lw=2)
plot!(heaviside, ls=:dash, label="step")

# This particular function takes any real number as input, and gives an output between $0$ and $1$. It is continuous and smooth.

#-

# #### Exercise 2
#
# Declare the sigmoid function above as an anonymous function with a different name.

#-

# #### Solution

#-

# Possible answers include
#
# ```julia
# x -> 1 / (1 + exp(-x))
# ```
#
# and
#
# ```julia
# σ = x -> 1 / (1 + exp(-x))
# ```

#-

# ### Mutating functions: `...!`
#
# To generate our plot of σ above, we used some functions that end with `!`. What does a `!` at the end of a function name mean in Julia?

#-

# Functions that change or modify their inputs are called **mutating functions**. But wait, don't all functions do that?
#
# No, actually. Functions typically take *inputs* and use those *inputs* to generate *outputs*, but the inputs themselves usually don't actually get changed by a function. For example, copy and execute the following code:
#
# ```julia
# v1 = [9, 4, 7, 11]
# v2 = sort(v1)
# ```

v1 = [9, 4, 7, 11]
v2 = sort(v1)

# `v2` is a sorted version of `v1`, but after calling `sort`, `v1` is still unsorted.
#
# However, now trying adding an exclamation point after `sort` and executing the following code:
#
# ```julia
# sort!(v1)
# ```

sort!(v1)

# Look at the values in `v1` now!

v1

# This time, the original vector itself was changed (mutated), and is now sorted. Unlike `sort`, `sort!` is a mutating function. Did the `!` make `sort!` mutating? Well, no, not really. In Julia, `!` indicates mutating functions by convention. When the author of `sort!` wrote `sort!`, they added a `!` to let you to know that `sort!` is mutating, but the `!` isn't what makes a function mutating or non-mutating in the first place.
#
# #### Exercise
#
# Some of our plotting commands end with `!`. Copy and execute the following code:
#
# ```julia
# r = -5:0.1:5
# g(x) = x^2
# h(x) = x^3
# plot(r, g, label="g")
# plot!(r, h, label="h")
# ```
#
# Then change the code slightly to remove the `!` after `plot!(r, h)`. How does this change your output? What do you think it means to add `!` after plotting commands?

#-

# #### Solution

r = -5:0.1:5
g(x) = x^2
h(x) = x^3
p1 = plot(r, g, label="g")
p2 = plot(r, h, label="h")

#-

p1

# ```julia
# plot(r, g)
# plot!(r, h)
# ```
# creates an overlay of `g` and `h`, whereas
#
# ```julia
# plot(r, g)
# plot(r, h)
# ```
#
# creates one plot for `g` and one for `h`.
#
# When we add a `!` after a plotting command, we are mutating or updating an *existing* plot.

#-

# ## Pointwise application of functions, `f.(x, y)` and `x .+ y` (known as "broadcasting")
#
# We saw in a previous notebook that we needed to add `.` after the names of some functions, as in
#
# ```julia
# green_amount = mean(Float64.(green.(apple)))
# ```
#
# What are those extra `.`s really doing?
#
# When we add a `.` after a function's name, we are telling Julia that we want to "**broadcast**" that function over the inputs passed to the function. This means that we want to apply that function *element-wise* over the inputs; in other words, it will apply the function to each element of the input, and return an array with the newly-calculated values.
#
# For example, copy and execute the following code:
# ```julia
# g.(r)
# ```
# Since the function `g` squares it's input, this squares all the elements of the range `r`.
#
# What happens if instead we just call `g` on `r` via
#
# ```julia
# g(r)
# ```
# ? Try this out and see what happens.

x = [1 2 3;4 5 6]
f(t) = √t + 5
f.(x)

# You should see an error message after calling `g(r)`, which says that Julia cannot multiply two vectors. When we call `g(r)`, we ask Julia to multiply `r` by `r`. When we call `g.(r)`, we ask Julia to multiply *each element* in `r` by itself.

#-

# #### Exercise 3
#
# Copy and execute the following code to get the type of the object `numbers = [1, 2, "three", 4.0]`:
#
# ```julia
# numbers = [1, 2, "three", 4.0]
# typeof(numbers)
# ```
#
# What is the type of `numbers`?

#-

# #### Solution

#-

# Its type is `Array{Any,1}`. This means that `numbers` is a
# 1-dimensional array that has elements of abstract type `Any`. We can think of the `Any` type as a classification that includes *all* concrete types, for example, `Float64`, `Int32`, `Bool`, `String`, `Char`, etc.

#-

# #### Exercise 4
#
# Broadcast `typeof` over `numbers` to see what the types of the elements stored inside `numbers` are.

#-

# #### Solution

#-

# We want to run
#
# ```julia
# typeof.(numbers)
# ```
#
# which will return
#
# ```julia
# 4-element Array{DataType,1}:
#  Int64
#  Int64
#  String
#  Float64
# ```
#
# The output is an array that contains one entry
# for each entry in `numbers`. For example, the output array tells us that `numbers[1]` is an `Int64` and `numbers[3]` is a `String`. Note that while the `numbers` array takes abstract type `Any`, the values stored in `numbers` have concrete (specific) types.

#-

# Note: Alternatively, we could have looked at the `typeof` each element in the array `numbers` using a `for` loop:
#
# ```julia
# for n in numbers
#     println(typeof(n))
# end
# ```

#-

# #### Exercise 5
#
# Write a `for` loop that applies `g` to each of the elements of `r` and prints the results. Verify that the numbers printed by this `for` loop are equal to the entries of the output of `g.(r)`.

#-

# #### Solution
#
# The requested `for` loop is
#
# ```julia
# r = -10:10
# for entry in r
#     println(g(entry))
# end
# ```
#
# which prints values that are equal to the entries of
#
# ```julia
# g.(r)
# ```
#
# We can check this via
#
# ```julia
# broadcasted_g = g.(r)
# i = 1
# for entry in r
#     println((g(entry) == broadcasted_g[i])
#     i += 1
# end
# ```
#
# which should print `true` 21 times (once for each entry in `broadcasted_g`).
#
# If you're familiar with conditional evaluation, we could instead do something like
#
# ```julia
# broadcasted_g = g.(r)
#
# # Create a boolean flag that gets updated to `false` if we detect a wrong answer:
# answers_match = true
#
# i = 1
# for entry in r
#     # Update our boolean flag if a wrong answer is detected
#     if (g(entry) != broadcasted_g[i])
#         answers_match = false
#     end
#     i += 1
# end
#
# # Use a ternary operator to report if any wrong answers were detected.
# answers_match ? println("Our `for` loop worked!") : println("Our `for` loop didn't give the same answers as broadcasting `g`.")
# ```

#-

# #### Exercise 6
#
# Define a range `xs` between -5 and 5 with steps of 0.5.
# Apply the $\sigma$ function pointwise to this range and define `ys` as the result.
# What does the result look like? Plot these as points and join them with lines.
#
# Make the plot interactive where you can vary the step size. Fix the range of the plot in the `x` and `y` directions using the functions `xlims!` and `ylims!`.

#-

# #### Solution

xs = -5:0.5:5
ys = σ.(xs)

#-

scatter(xs, ys)
plot!(xs, ys)

#-

## Try manipulating `stepsize` between 0 and 1:
stepsize = 0.5
xs = -5.0:stepsize:5.0
ys = σ.(xs)

scatter(xs, ys)
plot!(xs, ys)

xlims!(-5, 5)
ylims!(0, 1)

max_stars_repo_path: Courses/Foundations of machine learning/0300.Representing-data-with-models-quiz.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
using QuantumOptics
using Test, IonSim
using IonSim.analytical
using Suppressor

@suppress_err begin

@testset "Dynamics -- Rabi" begin
    # carrier
    C = Ca40(["S-1/2", "D-1/2"])
    L = Laser()
    chain = LinearChain(
            ions=[C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1])
            )
    T = Trap(configuration=chain, B=4e-4, Bhat=ẑ, δB=0, lasers=[L])
    mode = T.configuration.vibrational_modes.z[1]
    L.k = (x̂ + ẑ)/√2
    L.ϵ = ŷ
    Δf = transition_frequency(T, 1, ("S-1/2", "D-1/2"))
    L.Δ = Δf
    Ω = (rand() + 0.1) * 1e6
    Efield_from_rabi_frequency!(Ω, T, 1, 1, ("S-1/2", "D-1/2"))

    h = hamiltonian(T)
    tspan = 0:1e-3:4
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[0], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)

    @test isapprox(real.(ex), @.(sin(2π * Ω/2 * tout)^2), rtol=1e-6)

    # add detuning
    Δ = (rand() + 0.5) * 10e5
    L.Δ = Δf + Δ
    h = hamiltonian(T)
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[0], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)
    @test isapprox(real.(ex), @.((Ω^2/(Ω^2 + Δ^2)) * sin(2π * √(Ω^2 + Δ^2)/2 * tout)^2), rtol=1e-6)

    # add detuning using ion's stark_shift field
    L.Δ = Δf
    C.stark_shift["S-1/2"] = -Δ/2
    C.stark_shift["D-1/2"] = Δ/2
    h = hamiltonian(T)
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[0], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)
    @test isapprox(real.(ex), @.((Ω^2/(Ω^2 + Δ^2)) * sin(2π * √(Ω^2 + Δ^2)/2 * tout)^2), rtol=1e-6)

    # hot carrier
    zero_stark_shift(C)
    mode.N = 100
    ψi_ion = C["S-1/2"] ⊗ C["S-1/2"]'
    n̄ = rand(1:20)
    ψi_mode = thermalstate(mode, n̄)
    ψi = ψi_ion ⊗ ψi_mode
    tspan = 0:1e-3:4
    h = hamiltonian(T, lamb_dicke_order=0)
    @time tout, sol = timeevolution.schroedinger_dynamic(tspan, ψi, h)
    ex = expect(ionprojector(T, "D-1/2"), sol)
    η = get_η(mode, L, C)
    @test isapprox(real.(ex), rabi_flop(tout, Ω/2, η, n̄), rtol=1e-6)

    # sideband transitions

    ## under an RWA, RSB transitions, when starting in motional ground state,
    ## should be suppressed
    L.Δ = Δf - mode.ν
    h = hamiltonian(T, rwa_cutoff=1e5)
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[0], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)
    @test sum(ex) == 0

    ## a BSB should have a coupling strength ηΩ√(n+1)
    ## a RSB should have a coupling strength ηΩ√n
    L.Δ = Δf - mode.ν
    h = hamiltonian(T, rwa_cutoff=1e5)
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[1], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)
    @test isapprox(real.(ex), @.(sin(2π * η * Ω/2 * tout)^2), rtol=1e-6)

    L.Δ = Δf + mode.ν
    h = hamiltonian(T, rwa_cutoff=1e5)
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[1], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)
    @test isapprox(real.(ex), @.(sin(2π * sqrt(2) * η * Ω/2 * tout)^2), rtol=1e-6)
end

@testset "Dynamics -- time-dependent-parameters" begin
    # could use some better tests

    # use ϕ(t) to set laser frequency to resonance
    C = Ca40(["S-1/2", "D-1/2"])
    L = Laser()
    chain = LinearChain(
            ions=[C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1])
            )
    T = Trap(configuration=chain, B=4e-4, Bhat=ẑ, δB=0, lasers=[L])
    mode = T.configuration.vibrational_modes.z[1]
    L.k = (x̂ + ẑ)/√2
    L.ϵ = ŷ
    Δf = transition_frequency(T, 1, ("S-1/2", "D-1/2"))
    L.ϕ = t -> Δf * t
    Ω = (rand() + 0.1) * 1e6
    Efield_from_rabi_frequency!(Ω, T, 1, 1, ("S-1/2", "D-1/2"))

    h = hamiltonian(T)
    tspan = 0:1e-3:4
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[0], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)

    @test isapprox(real.(ex), @.(sin(2π * Ω/2 * tout)^2), rtol=1e-6)

    # set Ω(t) to a step function
    L.Δ = Δf
    L.ϕ = 0
    E = Efield_from_pi_time(Ω, T, 1, 1, ("S-1/2", "D-1/2"))
    L.E = t -> t < 1 ? 0 : E
    h = hamiltonian(T)
    tspan = 0:1e-3:3
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[0], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)

    delayed_sin(t) = t < 1 ? 0 : sin(2π * Ω/2 * (t-1))
    @test isapprox(real.(ex), @.(delayed_sin(tout)^2), rtol=1e-6)

    # check that δν(t) shifts sideband frequencies appropriately
    tspan = 0:1e-3:30
    mode.δν = t -> 20e3 * 1e-6
    L.Δ = Δf + mode.ν + 20e3
    h = hamiltonian(T, rwa_cutoff=1e5)
    tout, sol = timeevolution.schroedinger_dynamic(
            tspan, ionstate(T, "S-1/2") ⊗ mode[1], h
        )
    ex = expect(ionprojector(T, "D-1/2"), sol)
    η = get_η(mode, L, C)
                                                    # η -> η(t) also
    @test isapprox(real.(ex), @.(sin(2π * sqrt(2) * (η / sqrt(1.02)) * Ω/2 * tout)^2), rtol=1e-6)

    # δB(t)
    ## add test
end

@testset "molmer-sorensen" begin
    # slow gate
    C = Ca40(["S-1/2", "D-1/2"])
    L1 = Laser()
    L2 = Laser()
    chain = LinearChain(ions=[C, C], com_frequencies=(x=3e6,y=3e6,z=1e6), vibrational_modes=(;z=[1]))
    T = Trap(configuration=chain, B=4e-4, Bhat=(x̂ + ẑ)/√2, lasers=[L1, L2])
    mode = T.configuration.vibrational_modes.z[1]
    # Set the laser parameters
    ϵ = 40e3
    d = 80  # corrects for AC stark shift from single-photon coupling to sidebands
    Δf = transition_frequency(T, 1, ("S-1/2", "D-1/2"))
    L1.Δ = Δf + mode.ν + ϵ - d
    L1.k = ẑ
    L1.ϵ = x̂

    L2.Δ = Δf - mode.ν - ϵ + d
    L2.k = ẑ
    L2.ϵ = x̂

    mode.N = 15
    η = abs(get_η(mode, L1, C))
    Ω = √(1e3 * ϵ) / η  # This will give a 1kHz MS strength, since coupling goes like (ηΩ)^2/ϵ

    Efield_from_rabi_frequency!(Ω, T, 1, 1, ("S-1/2", "D-1/2"))
    Efield_from_rabi_frequency!(Ω, T, 2, 1, ("S-1/2", "D-1/2"))
    ψi = ionstate(T, "S-1/2", "S-1/2") ⊗ mode[0];  # initial state
    h = hamiltonian(T, rwa_cutoff=5e5)
    tspan = 0:0.25:1000
    @time tout, sol = timeevolution.schroedinger_dynamic(tspan, ψi, h);
    SS = expect(ionprojector(T, "S-1/2", "S-1/2"), sol)
    DD = expect(ionprojector(T, "D-1/2", "D-1/2"), sol)
    ex = analytical.two_ion_ms(tspan, 1e-6Ω, 1e-6mode.ν, 1e-6mode.ν + 1e-6ϵ, η, 0)
    @test isapprox(ex[1], SS, rtol=1e-2)
    @test isapprox(ex[2], DD, rtol=1e-2)

    # add checks for fast and hot gates
end
end  # end suppress


text:
include("PrimitiveVars.jl")
import Base.copy

mutable struct HydroSnapshot
    grid::Array{Float64, 1}
    cells::Array{PrimitiveVars, 1}
end

function midGrid(grid)::Array{Float64, 1}
    return [0.5*(grid[n]+grid[n+1])
            for n in 1:length(grid)-1]
end

function initHydroSnapshot(grid, d_func, p_func, v_func)
    cells = [PrimitiveVars(d_func(x),
                           p_func(x),
                           v_func(x))
             for x in midGrid(grid)]
    return HydroSnapshot(grid, cells)
end

function copy(source::HydroSnapshot)::HydroSnapshot
    return HydroSnapshot(source.grid, source.cells)
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
@testset "utils.jl" begin

import HTTP.Parsers
import HTTP.URIs

@test HTTP.Strings.escapehtml("&\"'<>") == "&"'<>"

@test URIs.isurlchar('\u81')
@test !URIs.isurlchar('\0')

for c = '\0':'\x7f'
    if c in ('.', '-', '_', '~')
        @test Parsers.ishostchar(c)
        @test Parsers.ismark(c)
        @test Parsers.isuserinfochar(c)
    elseif c in ('-', '_', '.', '!', '~', '*', '\'', '(', ')')
        @test Parsers.ismark(c)
        @test Parsers.isuserinfochar(c)
    else
        @test !Parsers.ismark(c)
    end
end

@test Parsers.isalphanum('a')
@test Parsers.isalphanum('1')
@test !Parsers.isalphanum(']')

@test Parsers.ishex('a')
@test Parsers.ishex('1')
@test !Parsers.ishex(']')

@test HTTP.Strings.tocameldash!("accept") == "Accept"
@test HTTP.Strings.tocameldash!("Accept") == "Accept"
@test HTTP.Strings.tocameldash!("eXcept-this") == "Except-This"
@test HTTP.Strings.tocameldash!("exCept-This") == "Except-This"
@test HTTP.Strings.tocameldash!("not-valid") == "Not-Valid"
@test HTTP.Strings.tocameldash!("♇") == "♇"
@test HTTP.Strings.tocameldash!("bλ-a") == "Bλ-A"
@test HTTP.Strings.tocameldash!("not fixable") == "Not fixable"
@test HTTP.Strings.tocameldash!("aaaaaaaaaaaaa") == "Aaaaaaaaaaaaa"
@test HTTP.Strings.tocameldash!("conTENT-Length") == "Content-Length"
@test HTTP.Strings.tocameldash!("Sec-WebSocket-Key2") == "Sec-Websocket-Key2"
@test HTTP.Strings.tocameldash!("User-agent") == "User-Agent"
@test HTTP.Strings.tocameldash!("Proxy-authorization") == "Proxy-Authorization"
@test HTTP.Strings.tocameldash!("HOST") == "Host"
@test HTTP.Strings.tocameldash!("ST") == "St"
@test HTTP.Strings.tocameldash!("X-\$PrototypeBI-Version") == "X-\$prototypebi-Version"
@test HTTP.Strings.tocameldash!("DCLK_imp") == "Dclk_imp"

for (bytes, utf8) in (
        (UInt8[], ""),
        (UInt8[0x00], "\0"),
        (UInt8[0x61], "a"),
        (UInt8[0x63, 0x61, 0x66, 0xe9, 0x20, 0x63, 0x72, 0xe8, 0x6d, 0x65], "café crème"),
        (UInt8[0x6e, 0x6f, 0xeb, 0x6c], "noël"),
        # (UInt8[0x6e, 0x6f, 0xeb, 0x6c, 0x20, 0xa4], "noël €"),
        (UInt8[0xc4, 0xc6, 0xe4], "ÄÆä"),
    )
    @test HTTP.Strings.iso8859_1_to_utf8(bytes) == utf8
end

# using StringEncodings
# println(encode("ÄÆä", "ISO-8859-15"))

end


text:
prob, opts = Problems.DoubleIntegrator()

solver = ALTROSolver(prob)
solver = ALTROSolver(prob, opts)
@test solver.opts === solver.solver_al.opts === solver.solver_al.solver_uncon.opts

# Pass in option arguments
solver = ALTROSolver(prob, opts, verbose=2, cost_tolerance=1)
@test solver.opts.verbose == 2
@test solver.opts.cost_tolerance == 1
@test solver.stats.parent == Altro.solvername(solver)

# Test other solvers
stats = solver.stats
al = Altro.AugmentedLagrangianSolver(prob, opts)
@test al.opts === opts
@test al.stats.parent == Altro.solvername(al)
al = Altro.AugmentedLagrangianSolver(prob, opts, stats)
@test al.stats === solver.stats
@test al.stats.parent == Altro.solvername(solver)
@test al.solver_uncon.stats.parent == Altro.solvername(solver)

# Try passing in a bad option
ilqr = Altro.iLQRSolver(prob, opts, something_wrong=false)
@test ilqr.opts === solver.opts

# Solve an unconstrained problem
ilqr = Altro.iLQRSolver(Problems.Cartpole(constrained=false)..., verbose=2)
b0 = benchmark_solve!(ilqr)
solver = ALTROSolver(Problems.Cartpole(constrained=false)..., verbose=2, projected_newton=false)
b1 = benchmark_solve!(solver)
@test iterations(solver) == iterations(ilqr)
@test cost(solver) ≈ cost(ilqr)
@test max_violation(solver) == 0
@test b0.allocs == 0
@test b1.allocs == 0
t0 = minimum(b0).time
t1 = minimum(b1).time
@test abs(t0-t1)/t1 < 0.1


text:
struct MyCurve end

@testset "Primitives" begin
    # Implement interface
    GeoInterfaceRFC.geomtype(::MyCurve) = GeoInterfaceRFC.LineString()
    GeoInterfaceRFC.ncoord(::GeoInterfaceRFC.LineString, geom::MyCurve) = 2
    GeoInterfaceRFC.ngeom(::GeoInterfaceRFC.LineString, geom::MyCurve) = 2
    GeoInterfaceRFC.getgeom(::GeoInterfaceRFC.LineString, geom::MyCurve, i) = [[1,2],[2,3]][i]

    # Test validity
    geom = MyCurve()
    @test test_interface_for_geom(geom)

    # Check functions
    @test GeoInterfaceRFC.npoint(geom) == 2  # defaults to ngeom
    @test_throws MethodError GeoInterfaceRFC.area(geom)

end


text:
"""
This module implements two related pieces of functionality:

* Interpolating ACSR conductor type and looking up appropriate data using current limit and base voltage information, and
* Implementing equations in IEEE Standard 738-2012 to compute transmission line thermal parameters.
"""
module LineThermalModel
export
    # ieee 738 functions
    eq2c_N_Re, eq3a_q_c1, eq3b_q_c2, eq4a_K_angle,
    eq6_T_film, eq8_q_s, eq9_theta, eq13a_mu_f, eq14a_p_f, eq15a_k_f,
    eq_omega, eq16a_H_c, eq16b_delta, eq17a_Z_c, eq17b_chi, eqtable2_C,
    eq18_Q_s, eq19_Q_se, eq20_K_solar, eq_eta_c, eq_eta_r,

    # acsr lookup/interpolation functions
    eq_mCp, acsr_table_english, acsr_table_si, acsr_interpolation, return_current_limits,
    estimate_length, ACSRSpecsEnglish, ACSRSpecsMetric

include("acsr.jl")
include("ieee738.jl")
include("powermodels.jl")
include("length.jl")

end # module


text:
module Placements

using ..Ahorn, Maple

displayName = "Placements"
group = "Placements"

drawingLayers = Ahorn.Layer[]

toolsLayer = nothing
targetLayer = nothing

material = nothing
materialName = nothing

scaleX = 1
scaleY = 1

selectionRect = nothing
previewGhost = nothing
clonedEntity = nothing

blacklistedCloneAttrs = ["id", "x", "y"]

placementLayers = String["entities", "triggers", "fgDecals", "bgDecals"]

function drawPlacements(layer::Ahorn.Layer, room::Ahorn.DrawableRoom, camera::Ahorn.Camera)
    ctx = Ahorn.getSurfaceContext(toolsLayer.surface)

    if selectionRect !== nothing && selectionRect.w > 0 && selectionRect.h > 0
        if isa(material, Ahorn.EntityPlacement) && material.placement == "rectangle"
            Ahorn.drawRectangle(ctx, selectionRect, Ahorn.colors.selection_selection_fc, Ahorn.colors.selection_selection_bc)
        end
    end

    if previewGhost !== nothing
        if Ahorn.layerName(targetLayer) == "entities"
            Ahorn.renderEntity(ctx, toolsLayer, previewGhost, room.room, alpha=Ahorn.colors.ghost_preplacement_alpha)
            Ahorn.renderEntitySelection(ctx, toolsLayer, previewGhost, room.room, alpha=Ahorn.colors.ghost_preplacement_alpha)

        elseif Ahorn.layerName(targetLayer) == "triggers"
            Ahorn.renderTrigger(ctx, toolsLayer, previewGhost, room.room, alpha=Ahorn.colors.ghost_preplacement_alpha)
            Ahorn.renderTriggerSelection(ctx, toolsLayer, previewGhost, room.room, alpha=Ahorn.colors.ghost_preplacement_alpha)

        elseif Ahorn.layerName(targetLayer) == "fgDecals" || Ahorn.layerName(targetLayer) == "bgDecals"
            Ahorn.drawDecal(ctx, previewGhost, alpha=Ahorn.colors.ghost_preplacement_alpha)
        end
    end
end

function generatePreview!(layer::Ahorn.Layer, material::Any, x, y; sx=1, sy=1, nx=x + 8, ny=y)
    if material !== nothing
        if layer.name == "entities" || layer.name == "triggers"
            # Use cache if possible, otherwise create a new entity/trigger

            placementsCache = layer.name == "entities" ? Ahorn.entityPlacementsCache : Ahorn.triggerPlacementsCache
            placements = layer.name == "entities" ? Ahorn.entityPlacements : Ahorn.triggerPlacements

            if materialName !== nothing
                return Ahorn.updateCachedEntityPosition!(placementsCache, placements, Ahorn.loadedState.map, Ahorn.loadedState.room, materialName, x, y, nx, ny)

            else
                if clonedEntity === nothing
                    global clonedEntity = Ahorn.generateEntity(Ahorn.loadedState.map, Ahorn.loadedState.room, material, x, y, nx, ny)
                end

                return Ahorn.updateEntityPosition!(clonedEntity, material, Ahorn.loadedState.map, Ahorn.loadedState.room, x, y, nx, ny)
            end

        elseif layer.name == "fgDecals" || layer.name == "bgDecals"
            texture = splitext(material)[1] * ".png"

            return Maple.Decal(texture, x, y, sx, sy)
        end
    end
end

function pushPreview!(layer::Ahorn.Layer, room::Maple.Room, preview::Any)
    if !get(Ahorn.config, "allow_out_of_bounds_placement", false)
        width, height = room.size
        x, y = Ahorn.position(preview)

        if x < 0 || x > width || y < 0 || y > height
            return false
        end
    end

    name = Ahorn.layerName(layer)

    # Update the id when we place the entity/trigger
    if name == "entities" || name == "triggers"
        preview.id = Ahorn.EntityIds.nextId()
    end

    # Make sure we don't have the same referance,
    # as the preview can be placed multiple times
    preview = deepcopy(preview)
    target = Ahorn.selectionTargets[name](room)

    Ahorn.History.addSnapshot!(Ahorn.History.RoomSnapshot("Placement ($name)", room))
    push!(target, preview)

    return true
end

function updateMaterialsEntities!()
    selectables = collect(keys(Ahorn.entityPlacements))
    sort!(selectables)

    if !isa(clonedEntity, Maple.Entity)
        wantedEntity = get(Ahorn.persistence, "placements_placements_entity", selectables[1])
        Ahorn.setMaterialList!(selectables, row -> row[1] == wantedEntity)
    end
end

function updateMaterialsTriggers!()
    selectables = collect(keys(Ahorn.triggerPlacements))
    sort!(selectables)

    if !isa(clonedEntity, Maple.Trigger)
        wantedTrigger = get(Ahorn.persistence, "placements_placements_trigger", selectables[1])
        Ahorn.setMaterialList!(selectables, row -> row[1] == wantedTrigger)
    end
end

function updateMaterialsDecals!()
    textures = Ahorn.decalTextures()

    wantedDecal = get(Ahorn.persistence, "placements_placements_decal", textures[1])
    Ahorn.setMaterialList!(textures, row -> row[1] == wantedDecal)
end

function cleanup()
    global previewGhost = nothing
    global selectionRect = nothing

    Ahorn.redrawLayer!(toolsLayer)
end

function toolSelected(subTools::Ahorn.ListContainer, layers::Ahorn.ListContainer, materials::Ahorn.ListContainer)
    wantedLayer = get(Ahorn.persistence, "placements_layer", "entities")
    Ahorn.updateLayerList!(placementLayers, row -> row[1] == wantedLayer)
    Ahorn.updateTreeView!(subTools, [])

    Ahorn.redrawingFuncs["tools"] = drawPlacements
    Ahorn.redrawLayer!(toolsLayer)
end

function layerSelected(list::Ahorn.ListContainer, materials::Ahorn.ListContainer, selected::String)
    prevLayer = Ahorn.layerName(targetLayer)
    global targetLayer = Ahorn.getLayerByName(drawingLayers, selected)

    if selected == "entities"
        updateMaterialsEntities!()

    elseif selected == "triggers"
        updateMaterialsTriggers!()

    elseif (selected == "fgDecals" || selected == "bgDecals") && prevLayer != "fgDecals" && prevLayer != "bgDecals" || isempty(materials.data)
        updateMaterialsDecals!()
    end

    global previewGhost = nothing

    Ahorn.persistence["placements_layer"] = selected
end

function subToolSelected(list::Ahorn.ListContainer, selected::String)

end

function materialSelected(list::Ahorn.ListContainer, selected::String)
    global materialName = selected
    layerName = Ahorn.layerName(targetLayer)

    if layerName == "entities"
        if haskey(Ahorn.entityPlacements, selected)
            global clonedEntity = nothing
            global material = Ahorn.entityPlacements[selected]
            Ahorn.persistence["placements_placements_entity"] = selected
        end

    elseif layerName == "triggers"
        if haskey(Ahorn.triggerPlacements, selected)
            global clonedEntity = nothing
            global material = Ahorn.triggerPlacements[selected]
            Ahorn.persistence["placements_placements_trigger"] = selected
        end

    elseif layerName == "fgDecals" || layerName == "bgDecals"
        global material = selected
        Ahorn.persistence["placements_placements_decal"] = selected
    end
end

function materialFiltered(list::Ahorn.ListContainer)
    layerName = Ahorn.layerName(targetLayer)

    if layerName == "entities" && isa(clonedEntity, Maple.Entity)
        return true
    end

    if layerName == "triggers" && isa(clonedEntity, Maple.Trigger)
        return true
    end

    Ahorn.selectRow!(list, row -> row[1] == materialName)
end

function updatePreviewGhost(x::Number, y::Number)
    targetX, targetY = x, y

    if !Ahorn.modifierControl()
        targetX = x * 8 - 8
        targetY = y * 8 - 8
    end

    prevGhost = deepcopy(previewGhost)
    newGhost = generatePreview!(targetLayer, material, targetX, targetY, sx=scaleX, sy=scaleY)

    if newGhost != prevGhost
        # No need to redraw if the target is on the same tile
        if isa(prevGhost, Maple.Entity) || isa(prevGhost, Maple.Trigger)
            if newGhost.data["x"] == prevGhost.data["x"] && newGhost.data["y"] == prevGhost.data["y"]
                return false
            end
        end

        if isa(prevGhost, Maple.Decal)
            if newGhost.x == prevGhost.x && newGhost.y == prevGhost.y
                return false
            end
        end

        global previewGhost = newGhost
        Ahorn.redrawLayer!(toolsLayer)

        return true
    end
end

function mouseMotion(x::Number, y::Number)
    if !Ahorn.modifierControl() && selectionRect === nothing
        updatePreviewGhost(x, y)
    end
end

function mouseMotionAbs(x::Number, y::Number)
    if Ahorn.modifierControl() && selectionRect === nothing
        updatePreviewGhost(x, y)
    end
end

function placementFunc(target::Union{Maple.Entity, Maple.Trigger})
    constructor = isa(target, Maple.Entity) ? Maple.Entity : Maple.Trigger

    return (x::Number, y::Number) -> constructor(target.name, x=x, y=y)
end

function middleClickAbs(x::Number, y::Number)
    layerName = Ahorn.layerName(targetLayer)
    selections = Ahorn.getSelected(Ahorn.loadedState.room, layerName, Ahorn.Rectangle(x, y, 1, 1))
    best = Ahorn.bestSelection(selections)

    if best !== nothing
        name, rect, target = best.layerName, best.rectangle, best.target

        if name == "fgDecals" || name == "bgDecals"
            global material = Ahorn.fixTexturePath(target.texture)
            Ahorn.selectMaterialList!(material)
            Ahorn.persistence["placements_placements_decal"] = material

        elseif name == "entities" || name == "triggers"
            horizontal, vertical = Ahorn.canResizeWrapper(target)
            placementType = horizontal || vertical ? "rectangle" : "point"

            # Fake a entity/trigger of the target
            # Copy all allowed attributes and store a fake (x, y) position for offseting nodes properly
            # If the entity/trigger is resizeable use "rectangle" placement rather than "point"

            global materialName = nothing
            global material = Ahorn.EntityPlacement(
                placementFunc(target),
                placementType,
                merge(
                    Dict{String, Any}((k, v) for (k, v) in deepcopy(target.data) if !(k in blacklistedCloneAttrs)),
                    Dict{String, Any}("__x" => target.data["x"], "__y" => target.data["y"])
                ),
                function(entity::Union{Maple.Entity, Maple.Trigger})
                    nodes = get(entity.data, "nodes", Tuple{Integer, Integer}[])
                    if length(nodes) > 0
                        x, y = entity.data["x"], entity.data["y"]
                        origx, origy = entity.data["__x"], entity.data["__y"]
                        newNodes = Tuple{Integer, Integer}[]

                        for node in nodes
                            nx, ny = node
                            push!(newNodes, (x + nx - origx, y + ny - origy))
                        end

                        entity.data["nodes"] = newNodes
                    end

                    delete!(entity.data, "__x")
                    delete!(entity.data, "__y")
                end
            )

            global clonedEntity = nothing
            updatePreviewGhost(x, y)
        end

        Ahorn.redrawLayer!(toolsLayer)
    end
end

function leftClick(x::Number, y::Number)
    if previewGhost !== nothing
        pushPreview!(targetLayer, Ahorn.loadedState.room, previewGhost)

        Ahorn.redrawLayer!(toolsLayer)
        Ahorn.redrawLayer!(targetLayer)
    end
end

function rightClickAbs(x::Number, y::Number)
    Ahorn.displayProperties(x, y, Ahorn.loadedState.room, targetLayer, toolsLayer)

    Ahorn.redrawLayer!(toolsLayer)
end

function selectionMotionAbs(x1::Number, y1::Number, x2::Number, y2::Number)
    if Ahorn.layerName(targetLayer) == "entities" || Ahorn.layerName(targetLayer) == "triggers"
        if material.placement == "line"
            if !Ahorn.modifierControl()
                x1, y1, x2, y2 = floor.(Int, [x1, y1, x2 + 8 * sign(x2 - x1), y2 + 8 * sign(y2 - y1)] ./ 8) .* 8
            end

            prevGhost = deepcopy(previewGhost)
            newGhost = generatePreview!(targetLayer, material, x1, y1, nx=x2, ny=y2)

            if newGhost != prevGhost
                global previewGhost = newGhost

                Ahorn.redrawLayer!(toolsLayer)
            end
        end
    end
end

function selectionMotionAbs(rect::Ahorn.Rectangle)
    if Ahorn.modifierControl()
        global selectionRect = rect

        if Ahorn.layerName(targetLayer) == "entities" || Ahorn.layerName(targetLayer) == "triggers"
            if material.placement == "rectangle"
                prevGhost = deepcopy(previewGhost)
                newGhost = generatePreview!(targetLayer, material, rect.x, rect.y, nx=rect.x + rect.w, ny=rect.y + rect.h)

                if newGhost != prevGhost
                    global previewGhost = newGhost

                    Ahorn.redrawLayer!(toolsLayer)
                end
            end
        end
    end
end

function selectionMotion(rect::Ahorn.Rectangle)
    if !Ahorn.modifierControl()
        ax1, ay1, ax2, ay2 = rect.x * 8 - 8, rect.y * 8 - 8, (rect.x + rect.w) * 8 - 8, (rect.y + rect.h) * 8 - 8

        global selectionRect = Ahorn.Rectangle(ax1, ay1, rect.w * 8, rect.h * 8)

        if Ahorn.layerName(targetLayer) == "entities" || Ahorn.layerName(targetLayer) == "triggers"
            if material.placement == "rectangle"
                prevGhost = deepcopy(previewGhost)
                newGhost = generatePreview!(targetLayer, material, ax1, ay1, nx=ax2, ny=ay2)

                if newGhost != prevGhost
                    global previewGhost = newGhost

                    Ahorn.redrawLayer!(toolsLayer)
                end
            end
        end
    end
end

# Doesn't matter if this is grid/abs, we only need to know the selection is done
function selectionFinish(rect::Ahorn.Rectangle)
    if previewGhost !== nothing
        pushPreview!(targetLayer, Ahorn.loadedState.room, previewGhost)

        global previewGhost = nothing
        global selectionRect = nothing

        Ahorn.redrawLayer!(toolsLayer)
        Ahorn.redrawLayer!(targetLayer)
    end
end

function layersChanged(layers::Array{Ahorn.Layer, 1})
    wantedLayer = get(Ahorn.persistence, "placements_layer", "entities")

    global drawingLayers = layers
    global toolsLayer = Ahorn.getLayerByName(layers, "tools")
    global targetLayer = Ahorn.selectLayer!(layers, wantedLayer, "entities")
end

decalScaleVals = (1.0, 2.0^4)

resizeModifiers = Dict{Integer, Tuple{Number, Number}}(
    # w, h
    # Decrease / Increase width
    Int('q') => (1, 0),
    Int('e') => (-1, 0),

    # Decrease / Increase height
    Int('a') => (0, 1),
    Int('d') => (0, -1)
)

scaleMultipliers = Dict{Integer, Tuple{Number, Number}}(
    # Vertical Flip
    Int('v') => (1, -1),

    # Horizontal Flip
    Int('h') => (-1, 1),
)

function keyboard(event::Ahorn.eventKey)
    redraw = false

    if haskey(scaleMultipliers, event.keyval)
        msx, msy = scaleMultipliers[event.keyval]

        global scaleX *= msx
        global scaleY *= msy

        redraw = true
    end

    if haskey(resizeModifiers, event.keyval)
        extraW, extraH = resizeModifiers[event.keyval]
        minVal, maxVal = decalScaleVals

        global scaleX = floor(Int, sign(scaleX) * clamp(abs(scaleX) * 2.0^extraW, minVal, maxVal))
        global scaleY = floor(Int, sign(scaleY) * clamp(abs(scaleY) * 2.0^extraH, minVal, maxVal))

        redraw = true
    end

    if redraw
        name = Ahorn.layerName(targetLayer)

        if name == "fgDecals" || name == "bgDecals"
            global previewGhost = generatePreview!(targetLayer, material, previewGhost.x, previewGhost.y, sx=scaleX, sy=scaleY)

            Ahorn.redrawLayer!(toolsLayer)
        end
    end

    return true
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
using SolverTools

"""
Return the status in GenericStats from a Stopping.
"""
function status_stopping_to_stats(stp :: AbstractStopping)
 stp_status = status(stp)
 convert = Dict([(:Optimal, :first_order),
                 (:SubProblemFailure, :unknown),
                 (:SubOptimal, :acceptable),
                 (:Unbounded, :unbounded),
                 (:UnboundedPb, :unbounded),
                 (:Stalled, :stalled),
                 (:IterationLimit, :max_iter),
                 (:Tired, :max_time), #disapear from Stopping.jl#v0.2.5
                 (:TimeLimit, :max_time),
                 (:ResourcesExhausted, :max_eval), #disapear from Stopping.jl#v0.2.5
                 (:EvaluationLimit, :max_eval),
                 (:ResourcesOfMainProblemExhausted, :max_eval),
                 (:Infeasible, :infeasible),
                 (:DomainError, :exception),
                 (:StopByUser, :user),
                 (:Exception, :exception),
                 (:Unknown, :unknown)
                 ])
 return convert[stp_status]
end

"""
Initialize a GenericStats from Stopping
"""
function stopping_to_stats(stp :: NLPStopping)
    return GenericExecutionStats(status_stopping_to_stats(stp),
                                 stp.pb,
                                 solution     = stp.current_state.x,
                                 objective    = stp.current_state.fx,
                                 primal_feas  = norm(stp.current_state.cx, Inf),
                                 dual_feas    = norm(stp.current_state.res, Inf),
                                 multipliers  = stp.current_state.lambda,
                                 iter         = stp.meta.nb_of_stop,
                                 elapsed_time = stp.current_state.current_time - stp.meta.start_time)
end

function stats_status_to_meta!(stp   :: AbstractStopping,
                               stats :: GenericExecutionStats)
  return stats_status_to_meta!(stp, stats.status)
end

function stats_status_to_meta!(stp    :: AbstractStopping,
                               status :: Symbol)
  #Update the meta boolean with the output message
  if status == :first_order stp.meta.optimal         = true end
  if status == :acceptable  stp.meta.suboptimal      = true end
  if status == :infeasible  stp.meta.infeasible      = true end
  if status == :small_step  stp.meta.stalled         = true end
  if status == :max_eval    stp.meta.resources       = true end
  if status == :max_iter    stp.meta.iteration_limit = true end
  if status == :max_time    stp.meta.tired           = true end
  if status ∈ [:neg_pred,
               :not_desc]   stp.meta.fail_sub_pb     = true end
  if status == :unbounded   stp.meta.unbounded       = true end
  if status == :user        stp.meta.stopbyuser      = true end
  if status ∈ [:stalled,
               :small_residual,
               :small_step]   stp.meta.stalled       = true end
  #if stats.status == :exception   stp.meta.exception       = true end #available ≥ 0.2.6

  return stp
end


text:
# The bearded devil
# Is forced
# To dwell
# In the only place
# Where they don't sell
# BURMA SHAVE

# http://burma-shave.org/jingles/

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
    # default = 0
    "--chain"
    help = "Write the chain to ~/web/ directory?"
    action = :store_true
    "config"
    help = "a YAML configuration file"
    required = true
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

outfmt(run_index::Int) = config["out_base"] * @sprintf("run%02d/", run_index)

# This code is necessary for multiple simultaneous runs on odyssey
# so that different runs do not write into the same output directory
if parsed_args["run_index"] == nothing
    run_index = 0
    outdir = outfmt(run_index)
    while ispath(outdir)
        println(outdir, " exists")
        run_index += 1
        outdir = outfmt(run_index)
    end
else
    run_index = parsed_args["run_index"]
    outdir = outfmt(run_index)
    println("Deleting old $outdir")
    run(`rm -rf $outdir`)
end

# make the output directory
println("Creating ", outdir)
mkdir(outdir)

using visibilities
# load data and figure out how many channels
dvarr = DataVis(config["data_file"])
nchan = length(dvarr)

if haskey(config, "exclude")
    exclude = config["exclude"]
    const global keylist = filter(x->(!in(x, exclude)), Int[i for i=1:nchan]) # which channels of the dset to fit
else
    const global keylist = Int[i for i=1:nchan]
end

println("keylist is ", keylist)

# go through any previously created directories and remove them
function cleardirs!(keylist::Vector{Int})
    println("Removing old directories")
    for key in keylist
        keydir = basedir * "jud$key"
        run(`rm -rf $keydir`)
    end
    println("Removed directories")
end

# Since we are starting with the -p option, it is no longer necessary to
# programatically add processes.
nchild = length(keylist)
# addprocs(nchild)
println("Workers allocated ", workers())

# make the values of run_index and config available on all processes
for process in procs()
    @spawnat process global run_id=run_index
    @spawnat process global cfg=config
end

@everywhere basefmt(id::Int) = cfg["base_dir"] * @sprintf("run%02d/", id)

@everywhere const global basedir = basefmt(run_id)

# make the internal Judith directory, if it doesn't exist
if !ispath(basedir)
    println("Creating ", basedir)
    mkdir(basedir)
end

# Clear all directories
cleardirs!(keylist)

@everywhere using constants
@everywhere using parallel
@everywhere using visibilities
@everywhere using image
@everywhere using gridding
@everywhere using model

# Delete the old log file (if it exists)
const logfile = outdir * "log.log"
if isfile(logfile)
    rm(logfile)
end

using Logging
# change the default logger
Logging.configure(filename=logfile, level=DEBUG)

debug("Created logfile.")

@everywhere function initfunc(key)

    # Load the relevant chunk of the dataset
    dset = DataVis(cfg["data_file"], key)
    # Conjugation is necessary for the SMA, methinks
    visibilities.conj!(dset) # Swap UV convention

    # Create a directory where all RADMC files will reside and be driven from
    keydir = basedir * "jud$key"
    mkdir(keydir)

    # Copy all relevant configuration scripts to this subdirectory
    # these are mainly setup files which will not change throughout the run
    run(`cp radmc3d.inp $keydir`)
    ag = basedir * "amr_grid.inp"
    run(`cp $ag $keydir`)
    run(`cp lines.inp $keydir`)
    run(`cp molecule_co.inp $keydir`)
    run(`cp wavelength_micron.inp $keydir`)

    # change the subprocess to reside in this directory for the remainder of the run
    # where it will drive its own independent RADMC3D process
    cd(keydir)

    # return the dataset
    return dset
end

# This is the likelihood function called by each individual process
@everywhere function f(dv::DataVis, key::Int, p::Parameters)

    # We are using the following convention: inclination ranges from
    # 0 to 180 degrees. 0 means face on, angular momentum vector pointing
    # at observer; 90 means edge on; and 180 means face on, angular momentum
    # vector pointing away from observer.
    # These are also the RADMC conventions.
    incl = p.incl # [deg]

    # We also adopt the RADMC convention for position angle, which defines position angle
    # by the angular momentum vector.
    # A positive PA angle means the disk angular momentum vector will be
    # rotated counter clockwise (from North towards East).
    PA = p.PA # [deg]

    vel = p.vel # [km/s]
    npix = cfg["npix"] # number of pixels, can alternatively specify x and y separately

    # Doppler shift the dataset wavelength to rest-frame wavelength
    beta = vel/c_kms # relativistic Doppler formula
    lam0 =  dv.lam * sqrt((1. - beta) / (1. + beta)) # [microns]

    # Run RADMC3D, redirect output to /dev/null
    run(`radmc3d image incl $incl posang $PA npix $npix lambda $lam0` |> DevNull)

    # Read the RADMC3D image from disk (we should already be in sub-directory)
    im = imread()

    # Convert raw image to the appropriate distance
    skim = imToSky(im, p.dpc)

    # Apply the gridding correction function before doing the FFT
    # shifts necessary as if the image were already offset
    corrfun!(skim, 1.0, p.mu_RA, p.mu_DEC) # alpha = 1.0 (relevant for spherical gridding function)
    # corrfun!(skim, 1.0, 0.0, 0.0) # alpha = 1.0 (relevant for spherical gridding function)

    # FFT the appropriate image channel
    vis_fft = transform(skim)

    # Apply the phase correction here, before we do the interpolation
    phase_shift!(vis_fft, p.mu_RA, p.mu_DEC)

    # Interpolate the `vis_fft` to the same locations as the DataSet
    mvis = ModelVis(dv, vis_fft)

    # Apply the phase correction here, since there are fewer data points
    # phase_shift!(mvis, p.mu_RA, p.mu_DEC)

    # Calculate chi^2 between these two
    return lnprob(dv, mvis)
end

# Create the model grid
grd = config["grid"]
global const grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)

# Regenerate all of the static files (e.g., amr_grid.inp)
# so that they may be later copied
debug("Writing grid")
write_grid(basedir, grid)
debug("Wrote grid")

debug("Initializing processes")
pipes = initialize(nchild, keylist, initfunc, f)
gather!(pipes)
debug("Initialized processes")

# Calculate the lnprior based upon the current parameter values
function lnprior(pars::Parameters)
    mu_d = 141. # [pc]
    sig_d = 7. # [pc]
    return -0.5 * (pars.dpc - mu_d)^2 / sig_d^2
end

# this function is called only on the main process, which proposes MCMC jumps
# to this function, and farms out the likelihood evaluation to all of the child
# processes
function fprob(p::Vector{Float64})

    # Here is where we make the distinction between a proposed vector of floats
    # (i.e., the parameters), and the object which defines all of the disk parameters
    # which every single subprocess will use

    # Parameters has the following definition (in model.jl)
    # M_star::Float64 # [g] stellar mass
    # r_c::Float64 # [cm] characteristic radius
    # T_10::Float64 # [K] temperature at 10 AU
    # q::Float64 # temperature gradient exponent
    # gamma::Float64 # surface temperature gradient exponent
    # M_CO::Float64 # [g] disk mass of CO
    # ksi::Float64 # [cm s^{-1}] microturbulence
    # dpc::Float64 # [pc] distance to system
    # incl::Float64 # [degrees] inclination 0 deg = face on, 90 = edge on.
    # PA::Float64 # [degrees] position angle (East of North)
    # vel::Float64 # [km/s] systemic velocity (positive is redshift/receeding)
    # mu_RA::Float64 # [arcsec] central offset in RA
    # mu_DEC::Float64 # [arcsec] central offset in DEC

    # Fix the following arguments: gamma, dpc
    gamma = 1.0 # surface temperature gradient exponent
    dpc = cfg["parameters"]["dpc"][1] # [pc] distance
    M_star = cfg["parameters"]["M_star"][1] # M_sun

    r_c, T_10, q, logM_CO, ksi, incl, PA, vel, mu_RA, mu_DEC = p

    # Enforce hard priors on physical parameters
    # Short circuit evaluation if we know the RADMC won't be valid.
    if ksi <= 0. || T_10 <= 0. || r_c <= 0.0 || M_star <= 0.0
        return -Inf
    end

    if incl < 0. || incl > 180.
        return -Inf
    end

    M_CO = 10^logM_CO

    # If we are going to fit with some parameters dropped out, here's the place to do it
    pars = Parameters(M_star, r_c, T_10, q, gamma, M_CO, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC)

    # Compute parameter file using model.jl, write to disk
    write_model(pars, basedir, grid)

    nd = basedir * "numberdens_co.inp"
    gv = basedir * "gas_velocity.inp"
    gt = basedir * "gas_temperature.inp"
    mt = basedir * "microturbulence.inp"

    # Copy new parameter files to all subdirectories
    for key in keylist
        keydir = basedir * "jud$key"
        run(`cp $nd $keydir`)
        run(`cp $gv $keydir`)
        run(`cp $gt $keydir`)
        run(`cp $mt $keydir`)
    end

    distribute!(pipes, pars)
    return gather!(pipes) # the summed lnprob
end

# wrapper for NLopt requires gradient as an argument (even if it's not used)
function fgrad(p::Vector, grad::Vector)
    val = fprob(p)
    debug(p, " : ", val)
    return val
end

function fp(p::Vector)
    val = fprob(p)
    debug(p, " : ", val)
    return val
end

debug("Initializing MCMC")
using Distributions
using PDMats

pp = config["parameters"]
# The parameters we'll be using
params = ["r_c", "T_10", "q", "logM_CO", "ksi", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
nparam = length(params)
starting_param = Array(Float64, nparam)
jumps = Array(Float64, nparam)

for i=1:nparam
    starting_param[i], jumps[i] = pp[params[i]]
end

jump_param = PDiagMat(jumps.^2)
jump_param = full(jump_param)

# Perturb the starting parameters
proposal = MvNormal(jump_param)
starting_param = starting_param .+ 7. * rand(proposal)

# If we've provided an empirically measured covariance matrix for the MCMC
# jump proposals, use that instead of our jumps without covariance
if haskey(config, "opt_jump")
    using NPZ
    jump_param = config["jump_scale"]^2 * npzread(config["opt_jump"])
end

# Now try optimizing the function using NLopt
# using NLopt
#
# nparam = length(starting_param)
# opt = Opt(:LN_COBYLA, nparam)
#
# max_objective!(opt, fgrad)
# ftol_abs!(opt, 0.05) # the precision we want lnprob to
#
# lower_bounds!(opt, lower)
# upper_bounds!(opt, upper)
#
# (optf,optx,ret) = optimize(opt, starting_param)
# println(optf, " ", optx, " ", ret)

using LittleMC

# Code to hook into the chain.js plot generator
if parsed_args["chain"]
    csv = open("/n/home07/iczekala/web/chain.js/mc.csv", "w")
else
    csv = open(outdir * "mc.csv", "w")
end

#Write the parameter header
writecsv(csv, params')

mc = MC(fp, config["samples"], starting_param, jump_param, csv)
debug("Initialized MCMC")

start(mc)

println(mean(mc.samples, 2))
println(std(mc.samples, 2))

runstats(mc)

debug("Acceptance: ", LittleMC.acceptance(mc))

LittleMC.write(mc, outdir * "mc.hdf5")
close(csv)

quit!(pipes)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using Printf

function showvampire{T<:Integer}(i::T, n::T, fangs::Array{T,2})
    s = @sprintf "%6d  %14d %s\n" i n join(fangs[1,:], "\u00d7")
    for i in 2:size(fangs)[1]
        s *= " "^23*join(fangs[i,:], "\u00d7")*"\n"
    end
    return s
end

vgoal = 25
vcnt = 0
dcnt = 0
println("Finding the first ", vgoal, " vampire numbers.")
println("     N         Vampire Fangs")
while vcnt < vgoal
    dcnt += 2
    for i in (10^(dcnt-1)):(10^dcnt-1)
        (isvampire, fangs) = vampirefangs(i)
        isvampire || continue
        vcnt += 1
        print(showvampire(vcnt, i, fangs))
        vcnt < vgoal || break
    end
end

test = [16758243290880, 24959017348650, 14593825548650]
println()
println("Checking a few numbers.")
println("     N         Vampire Fangs")
for (i, v) in enumerate(test)
    (isvampire, fangs) = vampirefangs(v)
    if isvampire
        print(showvampire(i, v, fangs))
    else
        println(@sprintf "%6d  %14d is not a vampire" i v)
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using LocalCoverage, Test

lockfile = "/tmp/testingLocalCoverage" # prevent infinite recursion when testing

covdir = normpath(joinpath(@__DIR__, "..", "coverage"))

if !isfile(lockfile)
    tracefile = joinpath(covdir, "lcov.info")
    @test !isfile(tracefile)
    touch(lockfile)
    generate_coverage("LocalCoverage"; genhtml = false)
    rm(lockfile)
    @test isfile(tracefile)
    clean_coverage("LocalCoverage")
    @test !isdir(covdir)
end


text:
#############################################################################
# variable.jl
# Defines Variable, which is a subtype of AbstractExpr
#############################################################################

"""
    VarType

Describe the type of a [`AbstractVariable`](@ref): either continuous (`ContVar`),
integer-valued (`IntVar`), or binary (`BinVar`).
"""
@enum VarType BinVar IntVar ContVar

@doc "Indicates a [`AbstractVariable`](@ref) is a binary variable." BinVar
@doc "Indicates a [`AbstractVariable`](@ref) is integer-valued." IntVar
@doc "Indicates a [`AbstractVariable`](@ref) is continuous." ContVar

"""
    abstract type AbstractVariable <: AbstractExpr end

An `AbstractVariable` should have `head` field, an `id_hash` field
and a `size` field to conform to the `AbstractExpr` interface, and
implement methods (or use the field-access fallbacks) for

* [`_value`](@ref), [`set_value!`](@ref): get or set the numeric value of the variable.
    `_value` should return `nothing` when no numeric value is set. Note:
    [`evaluate`](@ref) is the user-facing method to access the value of `x`.
* [`vexity`](@ref), [`vexity!`](@ref): get or set the `vexity` of the variable. The
    `vexity` should be `AffineVexity()` unless the variable has been
    [`fix!`](@ref)'d, in which case it is `ConstVexity()`.
* [`sign`](@ref), [`vartype`](@ref), and [`constraints`](@ref): get the `Sign`, `VarType`,
    numeric type, and a (possibly empty) vector of constraints which are
    to be applied to any problem in which the variable is used.

Optionally, also implement [`sign!`](@ref), [`vartype!`](@ref), and [`add_constraint!`](@ref)
to allow users to modify those values or add a constraint.

"""
abstract type AbstractVariable <: AbstractExpr end

# Default implementation of `AbstractVariable` interface
"""
    constraints(x::AbstractVariable)

Returns the current constraints carried by `x`.
"""
constraints(x::AbstractVariable) = x.constraints

"""
    add_constraint!(x::AbstractVariable, C::Constraint)

Adds an constraint to those carried by `x`.
"""
add_constraint!(x::AbstractVariable, C::Constraint) = push!(x.constraints, C)

"""
    vartype(x::AbstractVariable)

Returns the current [`VarType`](@ref) of `x`.
"""
vartype(x::AbstractVariable) = x.vartype

"""
    vartype!(x::AbstractVariable, vt::VarType)

Sets the current [`VarType`](@ref) of `x` to `vt`.
"""
vartype!(x::AbstractVariable, vt::VarType) = x.vartype = vt

"""
    vexity(x::AbstractVariable)

Returns the current vexity of `x`.
"""
vexity(x::AbstractVariable) = x.vexity

"""
    vexity!(x::AbstractVariable, v::Vexity)

Sets the current vexity of `x` to `v`. Should only be called
by [`fix!`](@ref) and [`free!`](@ref).
"""
vexity!(x::AbstractVariable, v::Vexity) = x.vexity = v

"""
    sign(x::AbstractVariable)

Returns the current sign of `x`.
"""
sign(x::AbstractVariable) = x.sign

"""
    sign!(x::AbstractVariable, s::Sign)

Sets the current sign of `x` to `s`.
"""
sign!(x::AbstractVariable, s::Sign) = x.sign = s

"""
    evaluate(x::AbstractExpr)

Returns the current value of `x` if assigned; errors otherwise.
"""
evaluate

function evaluate(x::AbstractVariable)
    return _value(x) === nothing ?
           error("Value of the variable is yet to be calculated") :
           output(_value(x))
end

"""
    _value(x::AbstractVariable)

Raw access to the current value of `x`; used internally by Convex.jl.
"""
_value(x::AbstractVariable) = x.value

"""
    set_value!(x::AbstractVariable, v)

Sets the current value of `x` to `v`.
"""
function set_value! end

set_value!(x::AbstractVariable, ::Nothing) = x.value = nothing

function set_value!(x::AbstractVariable, v::AbstractArray)
    size(x) == size(v) ||
        throw(DimensionMismatch("Variable and value sizes do not match!"))
    if iscomplex(x) && !(eltype(v) <: Complex)
        v = complex.(v)
    end
    return x.value = v
end

function set_value!(x::AbstractVariable, v::AbstractVector)
    size(x, 2) == 1 || throw(
        DimensionMismatch(
            "Cannot set value of a variable of size $(size(x)) to a vector",
        ),
    )
    size(x, 1) == length(v) ||
        throw(DimensionMismatch("Variable and value sizes do not match!"))
    if iscomplex(x) && !(eltype(v) <: Complex)
        v = complex.(v)
    end
    return x.value = v
end

function set_value!(x::AbstractVariable, v::Number)
    size(x) == (1, 1) ||
        throw(DimensionMismatch("Variable and value sizes do not match!"))
    if sign(x) == ComplexSign() && !(v isa Complex)
        v = complex(v)
    end
    return x.value = v
end
# End of `AbstractVariable` interface

# Some helpers for the constructors for `Variable`
iscomplex(sign::Sign) = sign == ComplexSign()
iscomplex(x::AbstractVariable) = iscomplex(sign(x))

# A concrete implementation of the `AbstractVariable` interface
mutable struct Variable <: AbstractVariable
    """
    Every `AbstractExpr` has a `head`; for a Variable it is set to `:variable`.
    """
    head::Symbol
    """
    A unique identifying hash used for caching.
    """
    id_hash::UInt64
    """
    The current value of the variable. Defaults to `nothing` until the
    variable has been [`fix!`](@ref)'d to a particular value, or the
    variable has been used in a problem which has been solved, at which
    point the optimal value is populated into this field.
    """
    value::ValueOrNothing
    """
    The size of the variable. Scalar variables have size `(1,1)`;
    `d`-dimensional vectors have size `(d, 1)`, and `n` by `m` matrices
    have size `(n,m)`.
    """
    size::Tuple{Int,Int}
    """
    `AffineVexity()` unless the variable is `fix!`'d, in which case it is
    `ConstVexity()`. Accessed by `vexity(v::Variable)`. To check if a
    `Variable` is fixed, use `vexity(v) == ConstVexity()`.
    """
    vexity::Vexity
    """
    The sign of the variable. Can be  `Positive()`, `Negative()`, `NoSign()`
    (i.e. real), or `ComplexSign()`. Accessed by `sign(v::Variable)`.
    """
    sign::Sign
    """
    Vector of constraints which are enforced whenever the variable is used
    in a problem.
    """
    constraints::Vector{Constraint}
    """
    The `VarType` of the variable (binary, integer, or continuous).
    """
    vartype::VarType
    function Variable(size::Tuple{Int,Int}, sign::Sign, vartype::VarType)
        if iscomplex(sign) && vartype != ContVar
            throw(
                ArgumentError(
                    "vartype must be `ContVar` for complex variables; got vartype = $vartype",
                ),
            )
        end
        this = new(
            :variable,
            0,
            nothing,
            size,
            AffineVexity(),
            sign,
            Constraint[],
            vartype,
        )
        this.id_hash = objectid(this)
        return this
    end
end
Variable(size::Tuple{Int,Int}, sign::Sign) = Variable(size, sign, ContVar)
function Variable(size::Tuple{Int,Int}, vartype::VarType)
    return Variable(size, NoSign(), vartype)
end
Variable(size::Tuple{Int,Int}) = Variable(size, NoSign(), ContVar)
ComplexVariable(size::Tuple{Int,Int}) = Variable(size, ComplexSign(), ContVar)

function Variable(m::Int, n::Int, sign::Sign, vartype::VarType)
    return Variable((m, n), sign, vartype)
end
Variable(m::Int, n::Int, sign::Sign) = Variable((m, n), sign, ContVar)
Variable(m::Int, n::Int, vartype::VarType) = Variable((m, n), NoSign(), vartype)
Variable(m::Int, n::Int) = Variable((m, n), NoSign(), ContVar)
ComplexVariable(m::Int, n::Int) = Variable((m, n), ComplexSign(), ContVar)

Variable(m::Int, sign::Sign, vartype::VarType) = Variable((m, 1), sign, vartype)
Variable(m::Int, sign::Sign) = Variable((m, 1), sign, ContVar)
Variable(m::Int, vartype::VarType) = Variable((m, 1), NoSign(), vartype)
Variable(m::Int) = Variable((m, 1), NoSign(), ContVar)
ComplexVariable(m::Int) = Variable((m, 1), ComplexSign(), ContVar)

Variable(sign::Sign, vartype::VarType) = Variable((1, 1), sign, vartype)
Variable(sign::Sign) = Variable((1, 1), sign, ContVar)
Variable(vartype::VarType) = Variable((1, 1), NoSign(), vartype)
Variable() = Variable((1, 1), NoSign(), ContVar)
ComplexVariable() = Variable((1, 1), ComplexSign(), ContVar)

###
# Compatability with old `sets` model
# Only dispatch to these methods when at least one set is given;
# otherwise dispatch to the inner constructors.
# There are choices of sign, no sign, or `ComplexVariable` and various
# ways to enter the size for each. We will start with:
# Case 1: with sign
function Variable(
    size::Tuple{Int,Int},
    sign::Sign,
    first_set::Symbol,
    more_sets::Symbol...,
)
    sets = [first_set]
    append!(sets, more_sets)
    if :Bin in sets
        vartype = BinVar
    elseif :Int in sets
        vartype = IntVar
    else
        vartype = ContVar
    end
    x = Variable(size, sign, vartype)
    if :Semidefinite in sets
        add_constraint!(x, x ⪰ 0)
    end
    return x
end
function Variable(
    m::Int,
    n::Int,
    sign::Sign,
    first_set::Symbol,
    more_sets::Symbol...,
)
    return Variable((m, n), sign, first_set, more_sets...)
end
function Variable(m::Int, sign::Sign, first_set::Symbol, more_sets::Symbol...)
    return Variable((m, 1), sign, first_set, more_sets...)
end
function Variable(sign::Sign, first_set::Symbol, more_sets::Symbol...)
    return Variable((1, 1), sign, first_set, more_sets...)
end

# Case 2: no sign
function Variable(size::Tuple{Int,Int}, first_set::Symbol, more_sets::Symbol...)
    return Variable(size, NoSign(), first_set, more_sets...)
end
function Variable(m::Int, n::Int, first_set::Symbol, more_sets::Symbol...)
    return Variable((m, n), NoSign(), first_set, more_sets...)
end
function Variable(m::Int, first_set::Symbol, more_sets::Symbol...)
    return Variable((m, 1), NoSign(), first_set, more_sets...)
end
function Variable(first_set::Symbol, more_sets::Symbol...)
    return Variable((1, 1), NoSign(), first_set, more_sets...)
end

# Case 3: ComplexVariable
function ComplexVariable(
    size::Tuple{Int,Int},
    first_set::Symbol,
    more_sets::Symbol...,
)
    return Variable(size, ComplexSign(), first_set, more_sets...)
end
function ComplexVariable(
    m::Int,
    n::Int,
    first_set::Symbol,
    more_sets::Symbol...,
)
    return Variable((m, n), ComplexSign(), first_set, more_sets...)
end
function ComplexVariable(m::Int, first_set::Symbol, more_sets::Symbol...)
    return Variable((m, 1), ComplexSign(), first_set, more_sets...)
end
function ComplexVariable(first_set::Symbol, more_sets::Symbol...)
    return Variable((1, 1), ComplexSign(), first_set, more_sets...)
end

Semidefinite(m::Integer) = Semidefinite(m, m)

function Semidefinite(m::Integer, n::Integer)
    if m == n
        x = Variable((m, m), NoSign(), ContVar)
        add_constraint!(x, x ⪰ 0)
        return x
    else
        error("Semidefinite matrices must be square")
    end
end

Semidefinite((m, n)::Tuple{Integer,Integer}) = Semidefinite(m, n)

function HermitianSemidefinite(m::Integer, n::Integer = m)
    if m == n
        x = Variable((m, n), ComplexSign(), ContVar)
        add_constraint!(x, x ⪰ 0)
        return x
    else
        error("`HermitianSemidefinite` matrices must be square")
    end
end

function HermitianSemidefinite((m, n)::Tuple{Integer,Integer})
    return HermitianSemidefinite(m, n)
end

function real_conic_form(x::AbstractVariable)
    vec_size = length(x)
    return sparse(1.0I, vec_size, vec_size)
end

function imag_conic_form(x::AbstractVariable)
    vec_size = length(x)
    if sign(x) == ComplexSign()
        return im * sparse(1.0I, vec_size, vec_size)
    else
        return spzeros(vec_size, vec_size)
    end
end

"""
    fix!(x::AbstractVariable, v = value(x))

Fixes `x` to `v`. It is subsequently treated as a constant in future
optimization problems. See also [`free!`](@ref).
"""
function fix! end

function fix!(x::AbstractVariable)
    _value(x) === nothing &&
        error("This variable has no value yet; cannot fix value to nothing!")
    vexity!(x, ConstVexity())
    return x
end

function fix!(x::AbstractVariable, v)
    set_value!(x, v)
    return fix!(x)
end

"""
    free!(x::AbstractVariable)

Frees a previously [`fix!`](@ref)'d variable `x`, to treat it once again as a
variable to optimize over.
"""
function free! end

function free!(x::AbstractVariable)
    vexity!(x, AffineVexity())
    return x
end


text:
using Clang
import SPIRV_Cross_jll

include_dir = joinpath(SPIRV_Cross_jll.artifact_dir, "include", "spirv_cross")
headers = joinpath.(Ref(include_dir), ["spirv_cross_c.h"])

# Set up include paths
clang_includes = [include_dir]

# Clang arguments
clang_extraargs = ["-v"]

# Test if a header should be wrapped
function wrap_header(top_hdr, cursor_header)
    any(startswith.(Ref(dirname(cursor_header)), clang_includes))
end

wc = init(;
                        headers,
                        output_file=joinpath(@__DIR__, "spirv_cross_api.jl"),
                        common_file=joinpath(@__DIR__, "spirv_cross_common.jl"),
                        clang_includes,
                        clang_args=clang_extraargs,
                        header_wrapped=wrap_header,
                        header_library=x -> "libspirv_cross",
                        clang_diagnostics=true,
                        )
run(wc)


text:
if !@isdefined sve_logistic
    const sve_logistic = Dict(
        10 => SparseIR.compute_sve(LogisticKernel(10.0)),
        42 => SparseIR.compute_sve(LogisticKernel(42.0)),
        10_000 => SparseIR.compute_sve(LogisticKernel(10_000.0)),
        (10_000, 1e-12) => SparseIR.compute_sve(LogisticKernel(10_000.0); ε=1e-12),
    )
end


text:
using StableRNGs
using VisualRegressionTests
using LinearAlgebra
using SparseArrays
using ImageMagick
using LightGraphs
using AbstractTrees

cd(@__DIR__)
include("functions.jl")
include("parse_readme.jl")

isci = "CI" ∈ keys(ENV)

default(show=false, reuse=true)

cd("../assets")

@testset "FIGURES" begin
    @plottest random_labelled_graph() "random_labelled_graph.png" popup=!isci tol=0.03

    @plottest random_3d_graph() "random_3d_graph.png" popup=!isci tol=0.02

    @plottest light_graphs() "light_graphs.png" popup=!isci tol=0.02

    @plottest directed() "directed.png" popup=!isci tol=0.02

    @plottest marker_properties() "marker_properties.png" popup=!isci tol=0.04

    @plottest edgelabel() "edgelabel.png" popup=!isci tol=0.1

    @plottest selfedges() "selfedges.png" popup=!isci tol=0.02

    @plottest multigraphs() "multigraphs.png" popup=!isci tol=0.1

    @plottest arc_chord_diagrams() "arc_chord_diagrams.png" popup=!isci tol=0.02

    @plottest ast_example() "ast_example.png" popup=!isci tol=0.02

    @plottest julia_type_tree() "julia_type_tree.png" popup=!isci tol=0.2

    @plottest julia_dict_tree() "julia_dict_tree.png" popup=!isci tol=0.2

    @plottest funky_edge_and_marker_args() "funky_edge_and_marker_args.png" popup=!isci tol=0.2

    @plottest custom_nodeshapes_single() "custom_nodeshapes_single.png" popup=!isci tol=0.2

    @plottest custom_nodeshapes_various() "custom_nodeshapes_various.png" popup=!isci tol=0.2
end

@testset "README" begin
    @plottest julia_logo_pun() "readme_julia_logo_pun.png" popup=!isci tol=0.2
end


text:
#Contains the functions for layer normalization! https://arxiv.org/abs/1607.06450

@everywhere using Statistics
if norm_active == true

    """
    Normalize
    """
    #Struct Normalized Activation
    #-------------------------------------------------------------------------------
    #Init the Net
    function init_θ_single_norm(net_s::Array{Int,1})
        x      = net_s[1]
        hidden = net_s[2:end-1]
        ysize  = net_s[end]

        w = []
        for y in [hidden...]
            push!(w, xavier(y, x))
            push!(w, ones(y))
            push!(w, zeros(y))
            x = y
        end
        push!(w, xavier(ysize, x)) # Prob actions
        push!(w, xavier(ysize, 1))

        return w
    end

    @everywhere function get_norm(x)
        x_mean = mean(x)
        x_std  = sqrt((x .- x_mean)' * (x .- x_mean) / length(x) + 1e-5)
        return x_mean, x_std
    end
    #@everywhere @zerograd get_norm(x)

    #Evaluate the Policy-Network
    @everywhere function mlp_nlayer_policy_norm(w::Any,x::AbstractVector,params::Params)
        σ       = params.σ
        u_range = params.u_range

        for i = 1:3:length(w)-2
            #1. Get the new summed input
            x = w[i] * x

            #2. Calculate the mean and std of the layer
            x_mean, x_std = get_norm(x)

            #3. Normalize activation by applying the norm and mean
            x_norm = (x .- x_mean) / x_std

            #4. Reparametrize the output
            x = x_norm .* w[i+1] + w[i+2]

            #5. Apply the non-linearity
            x = relu.(x)
        end

        mu_act   = (1 .- σ) .* u_range .* tanh.(w[end-1] * x .+ w[end])
        return vec(mu_act)
    end

    #Evaluate the Value-Network
    @everywhere function mlp_nlayer_value_norm(w::Any,x::AbstractVector,params::Params)
        eps = 1e-5

        for i = 1:3:length(w)-2
            #1. Get the new summed input
            x = w[i] * x

            #2. Calculate the mean and std of the layer
            x_mean, x_std = get_norm(x)

            #3. Normalize activation by applying the norm and mean
            x_norm = (x .- x_mean) / x_std

            #4. Reparametrize the output
            x = x_norm .* w[i+1] + w[i+2]

            #5. Apply the non-linearity
            x = relu.(x)
        end

        value    = (w[end-1] * x .+ w[end])[1]
        return value
    end

else

    """
    No normalization
    """

    init_θ_single_norm     = init_θ_single
    @everywhere mlp_nlayer_value_norm  = mlp_nlayer_value
    @everywhere mlp_nlayer_policy_norm = mlp_nlayer_policy

end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
#=
test_operations:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-06
=#

module TestOperations
    # packages
    using Test

    # extarnal modules
    include("add.jl")
    include("divide.jl")
    include("multiply.jl")
    include("subtract.jl")

    # method
    function test()
        @testset "Operations" begin
            @testset "Add" begin
                @test_nowarn Add.main()
            end
            @testset "Division" begin
                @test_nowarn Division.main()
            end
            @testset "Multiplication" begin
                @test_nowarn Multiplication.main()
            end
            @testset "Subtraction" begin
                @test_nowarn Subtraction.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using. TestOperations
    TestOperations.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui

function naive_listboxes(store::AbstractStore, get_state=Redux.get_state)
    if ReduxImGui.ListBox(store, s->get_state(s).listboxes["basic_listbox"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        item = ReduxImGui.get_item(get_state(store).listboxes["basic_listbox"])
        @info "Current item: $item"
    end
end


text:
module TestBetaDiv
  using Test
  using EcologicalNetworks

  A = convert(BinaryNetwork, web_of_life("A_HP_001"))
  B = convert(BinaryNetwork, web_of_life("A_HP_002"))

  @test βs(A, A).a == richness(A)
  @test βs(A, A).b == 0
  @test βs(A, A).c == 0
  @test βs(B, B).a == richness(B)
  @test βs(B, B).b == 0
  @test βs(B, B).c == 0

  @test βs(A,B).a+βs(A,B).b == richness(B)
  @test βs(A,B).a+βs(A,B).c == richness(A)

  @test βos(A, A).a == links(A)
  @test βos(A, A).b == 0
  @test βos(A, A).c == 0

  @test βwn(A, A).a == links(A)
  @test βwn(A, A).b == 0
  @test βwn(A, A).c == 0

  X = first(nz_stream_foodweb())

  @test βos(X, X).a == links(X)
  @test βos(X, X).b == 0
  @test βos(X, X).c == 0

end


text:
#### loads GTEX RSEM isoform level results, filtered for just HK1

function loadgtexgct(metafile="gtex_meta.tsv.gz", tpmfile="gtex_hk1.tsv.gz"; projdir=getprojectdir())

    meta = CSV.read(joinpath(projdir, "data", metafile), DataFrame)

    io = open(joinpath(projdir, "data", tpmfile)) |> GzipDecompressorStream

    verline = readline(io)

    rows, cols = parse.(Int, split(readline(io)))
    header = split(readline(io))

    M = zeros(rows, cols)

    genes = String[]
    trans = String[]

    for i = 1:rows
        line = readline(io)
        fields = split(line)
        push!(trans, fields[1])
        push!(genes, fields[2])

        for j = 1:cols
           M[i, j] = parse(Float64, fields[j+2])
        end

    end

    close(io)

    df = [DataFrame([trans, genes], header[1:2]) DataFrame(M, header[3:end])]
    sdf = stack(df, names(df)[3:end], variable_name=:SampleID, value_name=:TPM)

    tpm = innerjoin(meta, sdf, on=:SampleID)  ### isoform level
    genetpm = combine(groupby(tpm, [:SampleID, :Tissue, :gene_id]), :TPM => sum => :TPM) ## gene level

    tissues = sort(unique(genetpm.Tissue))

    tissueorder = Dict(t => i for (i, t) in enumerate(tissues));
    mediantissues = sort(combine(groupby(genetpm, :Tissue), :TPM => median => :Median), :Median, rev=true).Tissue
    mediantissueorder = Dict(t => i for (i, t) in enumerate(mediantissues))
    genetpm.TissueOrder = getindex.(Ref(tissueorder), genetpm.Tissue);
    genetpm.MedianOrder = getindex.(Ref(mediantissueorder), genetpm.Tissue);

    tpm.TissueOrder = getindex.(Ref(tissueorder), tpm.Tissue);
    tpm.MedianOrder = getindex.(Ref(mediantissueorder), tpm.Tissue);

    (tpm=tpm, genetpm=genetpm)
end

function gtexisoformheatmap(tpm; c=:blues, num_iso=3, kwargs...)

    tmm = @chain tpm begin
        groupby([:Tissue, :transcript_id])
        combine(:TPM => median => :TPM)
        unstack(:Tissue, :transcript_id, :TPM)
        coalesce.(_, 0)
    end

    M = Matrix(tmm[!, 2:end])
    sir = sortperm(maximum(M, dims=2) |> vec, rev=true)
    sic = sortperm(maximum(M, dims=1) |> vec, rev=true)

    ni = 1:num_iso
    heatmap(identity.(M[sir, sic[ni]]), yflip=true, c=c, yticks=(1:size(M, 1), tmm.Tissue[sir]), xticks=(ni, names(tmm)[2:end][sic[ni]]), xrotation=45,
        size=(500, 800), grid=false, left_margin=40mm, fontfamily="Helvetica", ylabel="Tissue", xlabel="HK1 Isoform", title="Top $(num_iso) most expressed HK1 Isoforms" ; kwargs...)

end

function gtexboxplot(genetpm; rev=false, orderfield=:MedianOrder)
    plot(grid=false)
    tissues = sort(unique(gtex.genetpm[!, [:Tissue, orderfield]]), orderfield, rev=rev).Tissue

    cc = cgrad(:viridis, length(tissues), categorical=true)
    order = genetpm[!, orderfield]
    order = ifelse(rev, maximum(order) .- order .+ 1, order)

    boxplot!(order, genetpm.TPM, group=genetpm[!, orderfield], c=hcat(cc...), marker=(stroke(0), 0.5, 2), size=(1000,600),  lab="")
    plot!(xticks=(1:length(tissues), tissues), xrotation=-45, bottom_margin=50mm, right_margin=10mm, fontfamily="helvetica", ylabel="TPM", title="GTEx HK1 Expression", left_margin=5mm, xlabel="Tissue")
end

function plot_hk1_isoforms(tpm, genecoords; num_iso=3, fs=10)
    tmm = @chain tpm begin
        groupby([:Tissue, :transcript_id])
        combine(:TPM => median => :TPM)
        unstack(:Tissue, :transcript_id, :TPM)
        coalesce.(_, 0)
    end

    M = Matrix(tmm[!, 2:end])
    sic = sortperm(maximum(M, dims=1) |> vec, rev=true)
    topisos = names(tmm)[sic[1:num_iso] .+ 1]

    topiso_ind = Int[]
    for t in topisos
        tn = replace(t, r"\.[0-9]$" => "")
        fi = findall(occursin.(tn, genecoords.TranscriptID))
        @assert length(fi) == 1
        push!(topiso_ind, first(fi))
    end

    p = plot()

    si = sortperm(topiso_ind, by = x -> genecoords.stop[x] - genecoords.start[x], rev=true)

    yos = 0:(num_iso-1)

    cc = cgrad(:viridis, num_iso+1, categorical=true)
    cc = [:steelblue, :orange, :black]
    for k in 1:num_iso
        yo = yos[k]
        i = topiso_ind[si[k]]
        plotgene!(genecoords.TranscriptID[i], genecoords, y=1 - yo , off=1, label=genecoords.TranscriptID[i], fs=fs, pos=:left, c=cc[k])
    end

    xl = Int.(round.(xlims()))
    xlt = string("chr10:", first(xl), "-", last(xl))
    plot!(size=(975, 150), bottom_margin=5mm, yticks=false, yaxis=false, left_margin=18mm, ylims=(-3, 3), xlabel=xlt, title="Top Expressed HK1 isoforms", fontfamily="helvetica")

end


text:
using DSP, AbstractFFTs

function frameSignal(z::Vector{<:Number}, windowFunction::Vector{<:Number}, frameAdvance::Int64)
    local zPad = copy(z)
    if iseven(length(windowFunction)); windowFunction=windowFunction[1:end-1]; end
    frameLength=length(windowFunction);
    padLen::Int64 = (frameLength-1)/2
    prepend!(zPad, zeros(ComplexF64,padLen))
    append!(zPad, zeros(ComplexF64,padLen))
    return [ zPad[i:frameLength+i-1].*windowFunction  for i in 1:frameAdvance:length(z)]
end

function STFT(z::Vector{<:Number}, windowFunction::Vector{<:Number}, frameAdvance::Int64; fs=1)
#References: Lim, Jae S., and Alan V. Oppenheim. Advanced topics in signal processing. Prentice-Hall, Inc., 1987.
    return ( fftshift.(fft.(frameSignal(z, windowFunction, frameAdvance))), fftshift(collect(AbstractFFTs.fftfreq(length(windowFunction), fs))) )
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
function predict_mu(gp::GPE, X::AbstractMatrix, cK::AbstractMatrix)
    n = size(X, 2)
    mu = mean(gp.mean, X) + cK*gp.alpha        # Predictive mean
    return mu
end

function _predict_raw(gp::GPE, X::AbstractMatrix)
    crossdata = KernelData(gp.kernel, gp.x, X)
    priordata = KernelData(gp.kernel, X, X)
    cK = cov(gp.kernel, gp.x, X, crossdata)
    mu = mean(gp.mean, X) + cK'*gp.alpha        # Predictive mean
    Lck = whiten!(gp.cK, cK)
    Sigma_raw = cov(gp.kernel, X, X, priordata)
    # Sigma_raw = Sigma_raw - Lck'Lck
    LinearAlgebra.BLAS.syrk!('U', 'T', -1.0, Lck, 1.0, Sigma_raw)
    LinearAlgebra.copytri!(Sigma_raw, 'U')
    return mu, Sigma_raw
end

""" Copy of a Gaussian process that shares X and cK but allows
    modification of the outcomes Y. This is useful for bootstrapping.
"""
function modifiable(gp::GPE)
    gp_copy = GPE(gp.mean, gp.kernel, gp.logNoise, gp.nobs,
        gp.x, copy(gp.y), gp.data,
        gp.dim, gp.cK, copy(gp.alpha),
        gp.mll, gp.mll, Float64[], Float64[])
    return gp_copy
end
function update_alpha!(gp::GPE)
    m = mean(gp.mean,gp.x)
    gp.alpha = gp.cK \ (gp.y - m)
end
function mll(gp::GPE)
    μ = mean(gp.mean,gp.x)
    return -dot((gp.y - μ),gp.alpha)/2.0 - logdet(gp.cK)/2.0 - gp.nobs*log(2π)/2.0 # Marginal log-likelihood
end

function addcov!(cK::AbstractMatrix{Float64}, kernel::Kernel, X::AbstractMatrix{Float64}, data::KernelData)
    dim, nobs = size(X)
    @assert size(cK, 1) == nobs
    @assert size(cK, 2) == nobs
    @inbounds for j in 1:nobs
        for i in 1:j
            cK[i,j] += cov_ij(kernel, X, X, data, i, j, dim)
            cK[j,i] = cK[i,j]
        end
    end
end
function add_diag!(mat::AbstractMatrix, d::Real)
    n = size(mat, 1)
    if n!=size(mat, 2)
        throw("Matrix is not square")
    end
    @inbounds for i in 1:n
        mat[i,i] += d
    end
    return mat
end
function update_chol!(pd::PDMats.PDMat)
    Σbuffer = mat(pd)
    Σbuffer, chol = make_posdef!(Σbuffer, cholfactors(pd))
    new_pd = wrap_cK(pd, Σbuffer, chol)
    return new_pd
end


text:
#! /usr/bin/julia

# Rosetta Code, Symmetric difference

a = Set(["John", "Bob", "Mary", "Serena"])
b = Set(["Jim", "Mary", "John", "Bob"])
c = symdiff(a, b)

println("Set a is: ", a)
println("Set b is: ", b)
println("Their symmetric difference is: ", c)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Test
using Dorothy

@DorothyAll()

@testset "Hierarchies" begin

	@testset "Chains" begin
		m1 = MolecularModel(10)
		m1.chainids =
				[fill("A", 3)..., fill("B", 2)..., fill("C", 4)..., "D"]
		chainiter = chains(m1)
		@test [length(chain) for chain in chainiter] == [3,2,4,1]
		@test parentindices(chainiter[1]) == 1:3
		@test parentindices(chainiter[2]) == 4:5
		@test parentindices(chainiter[3]) == 6:9
		@test parentindices(chainiter[4]) == 10:10
		@test chainat(m1, 1) == view(m1, 1:3)
		@test chainat(m1, 2) == view(m1, 1:3)
		@test chainat(m1, 3) == view(m1, 1:3)
		@test chainat(m1, 4) == view(m1, 4:5)
		@test chainat(m1, 5) == view(m1, 4:5)
		@test chainat(m1, 6) == view(m1, 6:9)
		@test chainat(m1, 7) == view(m1, 6:9)
		@test chainat(m1, 8) == view(m1, 6:9)
		@test chainat(m1, 9) == view(m1, 6:9)
		@test chainat(m1, 10) == view(m1, 10:10)
	end

	@testset "Residues" begin
		m1 = MolecularModel(20)
		m1.resids = [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5]
		m1.resnames = fill("WAT", 20)
		@test [length(res) for res in residues(m1)] == [4,4,4,4,4]
		m1.chainids = [fill("A", 10)..., fill("B", 10)...]
		resiter = residues(m1)
		@test [length(res) for res in resiter] == [4,4,2,2,4,4]
		@test parentindices(resiter[1]) == 1:4
		@test parentindices(resiter[2]) == 5:8
		@test parentindices(resiter[3]) == 9:10
		@test parentindices(resiter[4]) == 11:12
		@test parentindices(resiter[5]) == 13:16
		@test parentindices(resiter[6]) == 17:20
		@test residueat(m1, 1) == view(m1, 1:4)
		@test residueat(m1, 2) == view(m1, 1:4)
		@test residueat(m1, 3) == view(m1, 1:4)
		@test residueat(m1, 4) == view(m1, 1:4)
		@test residueat(m1, 5) == view(m1, 5:8)
		@test residueat(m1, 16) == view(m1, 13:16)
		@test residueat(m1, 17) == view(m1, 17:20)
		@test residueat(m1, 18) == view(m1, 17:20)
		@test residueat(m1, 19) == view(m1, 17:20)
		@test residueat(m1, 20) == view(m1, 17:20)
	end

	@testset "Fragments" begin
		m1 = MolecularModel(10)
		m1.topology = []
		pair!(m1.topology, (1, 2), (2, 3), (4, 5), (4, 6), (4, 7), (4, 8))
		fragiter = fragments(m1)
		@test [length(frag) for frag in fragiter] == [3,5,1,1]
		@test parentindices(fragiter[1]) == 1:3
		@test parentindices(fragiter[2]) == 4:8
		@test parentindices(fragiter[3]) == 9:9
		@test parentindices(fragiter[4]) == 10:10
		@test fragmentat(m1, 1) == view(m1, 1:3)
		@test fragmentat(m1, 2) == view(m1, 1:3)
		@test fragmentat(m1, 3) == view(m1, 1:3)
		@test fragmentat(m1, 4) == view(m1, 4:8)
		@test fragmentat(m1, 5) == view(m1, 4:8)
		@test fragmentat(m1, 6) == view(m1, 4:8)
		@test fragmentat(m1, 7) == view(m1, 4:8)
		@test fragmentat(m1, 8) == view(m1, 4:8)
		@test fragmentat(m1, 9) == view(m1, 9:9)
		@test fragmentat(m1, 10) == view(m1, 10:10)
	end

end # @testset