text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Plots

"""
Plots a NeuralFMU (ME).
"""
function fmiPlot(nfmu::ME_NeuralFMU)
    FMI.fmiPlot(nfmu.fmu, nfmu.solution, true)
end

"""
Plots a NeuralFMU (CS).
"""
function fmiPlot(nfmu::CS_NeuralFMU)
    numSig = length(nfmu.valueStack[1])-1
    t = collect(data[1] for data in nfmu.valueStack)

    fig = Plots.plot(xlabel="t [s]")
    for i in 1:numSig
        values = collect(data[i+1] for data in nfmu.valueStack)
        Plots.plot!(fig, t, values, label="NeuralFMU output #$(i)")
    end

    fig
end


text:
# Takes an output name which will be made into a directory
# The dataframe and associated information exported out
# and appropriate plots made
# possibly a web happy report made
# this is then zipped up
#s

include(Pkg.dir("BismarkSummary","src","bismark-report.jl"))
test_bismark_report1 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report.txt")
test_bismark_report2 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report1.txt")
test_bismark_report3 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report2.txt")
bismark_report_filenames = [test_bismark_report1,test_bismark_report2,test_bismark_report3]
sampleinfo_path = joinpath(Pkg.dir(), "BismarkSummary","testdata","datasource.tsv")
outputdir = Pkg.dir("BismarkSummary","testdata","bismarkreport1")
groupby_fields = ["strain"]

make_bismark_summary_report(outputdir, bismark_report_filenames, sampleinfo_path, groupby_fields)

ext: jl


text:
export iauPvdpv
"""
Inner (=scalar=dot) product of two pv-vectors.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   a        double[2][3]      first pv-vector
   b        double[2][3]      second pv-vector

Returned:
   adb      double[2]         a . b (see note)

Note:

   If the position and velocity components of the two pv-vectors are
   ( ap, av ) and ( bp, bv ), the result, a . b, is the pair of
   numbers ( ap . bp , ap . bv + av . bp ).  The two numbers are the
   dot-product of the two p-vectors and its derivative.

Called:
   iauPdp       scalar product of two p-vectors

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPvdpv(double a[2][3], double b[2][3], double adb[2])
function iauPvdpv(a::AbstractMatrix{<:Real}, b::AbstractMatrix{<:Real})
   # Preallocate return values
   adb = zeros(Float64, 2)

   ccall((:iauPvdpv, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          convert(Array{Float64, 2}, a'),
          convert(Array{Float64, 2}, b'),
          adb)

   return SVector{2}(adb)
end


text:
using StatsFuns
export Normal, HalfNormal

# The `@parameterized` macro below does three things:
# 1. Defines the `Normal` struct
#
#    struct Normal{N, T} <: (ParameterizedMeasure){N}
#        par::MeasureTheory.NamedTuple{N, T}
#    end
#
#    (Note that this is the same form required by the `@kwstruct` in the
#    KeywordCalls package.)
#
# 2. Sets the base measure
#
#    basemeasure(::Normal) = (1 / sqrt2π) * Lebesgue(ℝ)
#
#    This just means that (log-)densities will be defined relative to this.
#    Including the `(1 / sqrt2π)` in the base measure allows us to avoid
#    carrying around the normalization constant.
#
# 3. Sets up a method for a call with two unnamed arguments,
#
#    Normal(μ, σ) = Normal((μ=μ, σ=σ))
#
@parameterized Normal(μ,σ) ≪ (1/sqrt2π) * Lebesgue(ℝ)

@kwstruct Normal(μ,Σ)
@kwstruct Normal(μ,Σ⁻¹)

TV.as(::Normal) = asℝ

# `@kwalias` defines some alias names, giving users flexibility in the names
# they use. For example, σ² is standard notation for the variance parameter, but
# it's a lot to type. Some users might prefer to just use `var` and have us do
# the conversion (at compile time).
@kwalias Normal [
    mean  => μ
    mu    => μ
    std   => σ
    sigma => σ
    var   => σ²
    cov   => Σ
]

# It's often useful to be able to map into the parameter space for a given
# measure. We (currently, at least) do this *independently* per parameter. This
# allows us to do things like, e.g.
#
# julia> transform(asparams(Normal(2,3)))(randn(2))
# (μ = -0.6778046205440658, σ = 3.7357686957861898)
#
# julia> transform(asparams(Normal(μ=-3,logσ=2)))(randn(2))
# (μ = 0.3995295982002209, logσ = -1.3902312393777492)
#
# Or using types:
# julia> transform(asparams(Normal{(:μ,:σ²)}))(randn(2))
# (μ = -0.4548087051528626, σ² = 11.920775478312793)
#
# And of course, you can apply `Normal` to any one of the above.
#
asparams(::Type{<:Normal}, ::Val{:σ²}) = asℝ₊
asparams(::Type{<:Normal}, ::Val{:τ}) = asℝ₊
asparams(::Type{<:Normal}, ::Val{:logτ}) = asℝ

# Rather than try to reimplement everything in Distributions, measures can have
# a `distproxy` method. This just delegates some methods to the corresponding
# Distributions.jl methods. For example,
#
#     julia> entropy(Normal(2,4))
#     2.805232894324563
#
distproxy(d::Normal{(:μ, :σ)}) = Dists.Normal(d.μ, d.σ)

###############################################################################
# Some distributions have a "standard" version that takes no parameters
@kwstruct Normal()

# Instead of setting default values, the `@kwstruct` call above makes a
# parameter-free instance available. The log-density for this is very efficient.
logdensity(d::Normal{()} , x) = - x^2 / 2

Base.rand(rng::Random.AbstractRNG, T::Type, μ::Normal{()}) = randn(rng, T)

###############################################################################
# `μ` and `σ` parameterizations are so common, we have a macro to make them easy
# to build. Note that `μ` and `σ` are *not* always mean and standard deviation,
# but instead should be considered "location" and "scale" parameters,
# respectively. Also, it's valid to have existing parameters here as a starting
# point, in particular a "shape parameter" is very common. See `StudentT` for an
# example of this.
#
# This defines methods for `logdensity` and `Base.rand`. The latter works out
# especially nicely, because it lets us do things like
#
#     julia> using Symbolics
#     [ Info: Precompiling Symbolics [0c5d862f-8b57-4792-8d23-62f2024744c7]
#
#     julia> @variables μ σ
#     2-element Vector{Num}:
#      μ
#      σ
#
#     julia> rand(Normal(μ,σ))
#     μ + 1.2517620265570832σ
#
@μσ_methods Normal()

###############################################################################
# The `@half` macro takes a symmetric univariate measure and efficiently creates
# a truncated version.
@half Normal()
@kwstruct HalfNormal()

# @μσ_methods, without the μ.
@σ_methods HalfNormal()

# A single unnamed parameter for `HalfNormal` should be interpreted as a `σ`
HalfNormal(σ) = HalfNormal(σ = σ)

###############################################################################
@kwstruct Normal(μ,σ²)

function logdensity(d::Normal{(:σ²)}, x)
    σ² = d.σ²
    -0.5 * (log(σ²) + (x^2/σ²))
end

function logdensity(d::Normal{(:μ,:σ²)}, x)
    μ = d.μ
    σ² = d.σ²
    -0.5 * (log(σ²) + ((x - μ)^2/σ²))
end

###############################################################################
@kwstruct Normal(μ,τ)

function logdensity(d::Normal{(:τ)}, x)
    τ = d.τ
    0.5 * (log(τ) - τ * x^2)
end

function logdensity(d::Normal{(:μ,:τ)}, x)
    μ = d.μ
    τ = d.τ
    0.5 * (log(τ) - τ * (x - μ)^2)
end

###############################################################################
@kwstruct Normal(μ, logσ)

function logdensity(d::Normal{(:μ,:logσ)}, x)
    μ = d.μ
    logσ = d.logσ
    -logσ - 0.5(exp(-2logσ)*((x - μ)^2))
end


text:
using HydroModels
using Base.Test


text:
module num

using DSP

export trapz, conv_same, convergent, mixing

function trapz(x,y)
    dx = x[2] - x[1]
    intgrl = 0.0
    @fastmath @inbounds for i in 1:length(y)-1
        intgrl += (y[i+1] + y[i])
    end
    return intgrl * 0.5 * dx
end

function conv_same(A,B)
    C = conv(A,B)
    lower_bound = floor(Int,length(C)/4)
    upper_bound = floor(Int,3*length(C)/4-1)
    lengthc = upper_bound-lower_bound + 1
    if lengthc != length(A)
        upper_bound = floor(Int,3*length(C)/4)
    end
    return C[lower_bound:upper_bound]
end

function convergent(value_old,value_new,ω,nω,tol)
    error = zeros(Float64,nω)
    if length(value_new) != nω
        for i = 1:nω
            error[i] = abs( sum( imag.(value_new[i,:]) .- imag.(value_old[i,:]) ) )
        end
    else
        for i = 1:nω
            error[i] = abs( imag.(value_new[i]) .- imag.(value_old[i]) )
        end
    end
    err = trapz(ω,error)
    convg = err <= tol
    return convg,err
end

function mixing(value_old,value_new,mixings)
    return (1 - mixings) .* value_old .+ mixings .* value_new
end

end

max_stars_repo_name: muhammadgaffar/SKRIPSI

max_issues_repo_name: muhammadgaffar/SKRIPSI

max_forks_repo_name: muhammadgaffar/SKRIPSI


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

module TestsRunner
    # packages
    using Test

    # external modules
    include("../HelloWorld/test_hello_world.jl")
    include("../Module/test_module.jl")
    include("../Operations/test_operations.jl")
    include("../ComplexNumber/test_complex_number.jl")
    include("../String/test_string.jl")
    include("../Control/test_control.jl")
    include("../Type/test_types.jl")
    include("../Collection/test_collection.jl")
    include("../Array/test_array.jl")
    include("../Plots/test_plots.jl")
    include("../DataFrames/test_df.jl")
    include("../StacksAndQueues/test_stack_queue.jl")
    include("../Trees/test_trees.jl")
    include("../LinkedLists/test_linked_lists.jl")
    include("../Search/test_search.jl")

    # methods
    function tests()
        @testset "JuliaPractice" begin
            TestHelloWorld.test()
            TestModule.test()
            TestOperations.test()
            TestCompNum.test()
            TestString.test()
            TestControl.test()
            TestTypes.test()
            TestCollection.test()
            TestArray.test()
            TestPlots.test()
            TestDf.test()
            TestStackQueue.test()
            TestTrees.test()
            TestLinkedLists.test()
            TestSearch.test()
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestsRunner
    @time TestsRunner.tests()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
function binary_string_to_decimal(binstring)
    total = 0
    num_bits = length(binstring)
    for i in 1:num_bits
        if binstring[i] == '1'
            total += 2^(num_bits-i)
        end
    end
    total
end

# read the input data into an array (strips newlines for us)
rows = readlines("input.txt")

# how many binary digits in each row
num_digits = length(rows[1])
# initialize arrays to count ones and zeros
one_counts = zeros(num_digits)
zero_counts = zeros(num_digits)
# count ones and zeros in each position
for row in rows
    for i in 1:num_digits
        if row[i] == '1' # has to be single quotes !?
            one_counts[i] += 1
        else
            zero_counts[i] += 1
        end
    end
end
println("zero counts", zero_counts)
println("one counts", one_counts)

# now construct binary strings depending on whether there were
# more ones or zeros in each position
gamma_rate = ""
epsilon_rate = ""

for i in 1:num_digits
    if one_counts[i] > zero_counts[i]
        global gamma_rate *= "1"
        global epsilon_rate *= "0"
    else # what if they're equal??? never mind....
        global epsilon_rate *= "1"
        global gamma_rate *= "0"
    end
end

println("gamma rate ", gamma_rate)
println("epsilon rate ", epsilon_rate)

# convert to decimal
gamma_rate_dec = binary_string_to_decimal(gamma_rate)
epsilon_rate_dec = binary_string_to_decimal(epsilon_rate)
println("gamma rate (decimal) ",gamma_rate_dec)
println("epsilon rate (decimal) ",epsilon_rate_dec)

# multiply
println("total power consumption ",gamma_rate_dec * epsilon_rate_dec)


text:
@testset "Gen&Opt" begin
    isdir("basic") && rm("basic", recursive=true, force=true)
    newsite("basic")

    global silly_call = 0
    function foo_on_write(pg, vars)
        global silly_call
        silly_call += length(pg)
        vars["fd_rpath"] # should not error
        return nothing
    end

    serve(single=true, on_write=foo_on_write)

    @test silly_call > 0

    # ---------------
    @test all(isdir, ("_assets", "_css", "_libs", "_layout", "__site"))
    @test isfile(joinpath("__site", "index.html"))
    @test isfile(joinpath("__site", "menu1", "index.html"))
    @test isfile(joinpath("__site", "menu2", "index.html"))
    @test isfile(joinpath("__site", "menu3", "index.html"))
    @test isfile(joinpath("__site", "css", "basic.css"))
    @test isfile(joinpath("__site", "css", "franklin.css"))

    # ---------------
    if Franklin.FD_CAN_MINIFY
        presize1 = stat(joinpath("__site", "css", "basic.css")).size
        presize2 = stat(joinpath("__site", "index.html")).size
        optimize(prerender=false)
        @test stat(joinpath("__site", "css", "basic.css")).size < presize1
        @test stat(joinpath("__site", "index.html")).size < presize2
    end
    # ---------------
    # verify all links
    Franklin.verify_links()

    # ---------------
    # change the prepath
    index = read(joinpath("__site","index.html"), String)
    @test occursin("=\"/css/basic.css", index)
    @test occursin("=\"/css/franklin.css", index)
    @test occursin("=\"/libs/highlight/github.min.css", index)
    @test occursin("=\"/libs/katex/katex.min.css", index)

    optimize(minify=false, prerender=false, prepath="prependme")
    index = read(joinpath("__site","index.html"), String)
    @test occursin("=\"/prependme/css/basic.css", index)
    @test occursin("=\"/prependme/css/franklin.css", index)
    @test occursin("=\"/prependme/libs/highlight/github.min.css", index)
    @test occursin("=\"/prependme/libs/katex/katex.min.css", index)
end

if F.FD_CAN_PRERENDER; @testset "prerender" begin
    @testset "katex" begin
        hs = raw"""
        <!doctype html>

        <meta charset=UTF-8>

range is \(10\sqrt{3}\)–\(20\sqrt{2}\)

Consider an invertible matrix \(M\) made of blocks \(A\), \(B\), \(C\) and \(D\) with

        \[ M \quad\!\! =\quad\!\! \begin{pmatrix} A & B \\ C & D \end{pmatrix} \]

        """

        jskx = F.js_prerender_katex(hs)
        # conversion of the non-ascii endash (inline)
        @test occursin("""–""", jskx)
        # # conversion of `\(M\)` (inline)
        # @test occursin("""<math xmlns=\"http://www.w3.org/1998/Math/MathML\"><semantics><mrow><mi>M</mi></mrow>""", jskx)
        # # conversion of the equation (display)
        # @test occursin("""<math xmlns=\"http://www.w3.org/1998/Math/MathML\"><semantics><mrow><mi>M</mi>""", jskx)
    end

    if F.FD_CAN_HIGHLIGHT; @testset "highlight" begin
        hs = raw"""
        <!doctype html>

        <meta charset=UTF-8>

Title

Blah

using Test
        # Woodbury formula
        b = 2
        println("hello $b")

        """
        jshl = F.js_prerender_highlight(hs)
        # conversion of the code
        @test occursin("""
using""", jshl)
        @test occursin(raw"""# Woodbury formula""", jshl)
    end; end # if can highlight
end; end # if can prerender


text:
"""
Package: Forecast

    london()

Return ten years of monthly data about weather and crime in Greater London from 2008 to 2018.

Data has been collected and joined from london.gov.uk and metoffice.gov.uk (Heathrow Station).

# Weather Variables

- `MaxTemp`:  Mean daily maximum temperature in C°
- `MinTemp`:  Mean daily minimum temperature in C°
- `AirFrost`: Days of air frost
- `Rain`:     Total rainfall in mm
- `Sun`:      Total sunshine durationin hours

# Crime Variables and its aggregated categories
```
┌─────────────────────────────┬────────────────────────────────────────┐
│                       Crime │                               Category │
├─────────────────────────────┼────────────────────────────────────────┤
│                    Burglary │            Burglary in Other Buildings │
│                             │                 Burglary in a Dwelling │
│                      Damage │            Criminal Damage To Dwelling │
│                             │       Criminal Damage To Motor Vehicle │
│                             │      Criminal Damage To Other Building │
│                             │                  Other Criminal Damage │
│                       Drugs │                       Drug Trafficking │
│                             │                            Other Drugs │
│                             │                    Possession Of Drugs │
│                       Fraud │                     Counted per Victim │
│                             │                  Other Fraud & Forgery │
│                       Other │                         Going Equipped │
│                             │                       Other Notifiable │
│                     Robbery │                      Business Property │
│                             │                      Personal Property │
│                      Sexual │                           Other Sexual │
│                             │                                   Rape │
│                       Theft │                  Handling Stolen Goods │
│                             │ Motor Vehicle Interference & Tampering │
│                             │                            Other Theft │
│                             │                     Other Theft Person │
│                             │               Theft From Motor Vehicle │
│                             │                       Theft From Shops │
│                             │          Theft/Taking Of Motor Vehicle │
│                             │            Theft/Taking of Pedal Cycle │
│                    Violence │                    Assault with Injury │
│                             │                         Common Assault │
│                             │                   Grievous Bodily Harm │
│                             │                             Harassment │
│                             │                                 Murder │
│                             │                       Offensive Weapon │
│                             │                         Other violence │
│                             │                           Wounding/GBH │
└─────────────────────────────┴────────────────────────────────────────┘
```

# Returns
Dataframe containing the descrived dataset.

# Examples
```julia-repl
julia> london()
132×15 DataFrame
 Row │ Date        MaxTemp  MinTemp  AirFrost
     │ Date        Float64  Float64  Int64
─────┼─────────────────────────────────────────[...]
   1 │ 2008-01-01     10.4      4.7         0
   2 │ 2008-02-01     11.0      2.0         7
   3 │ 2008-03-01     10.6      3.7         2
   [...]
```
"""
function london()

    path = art_path("london")

    df = open(joinpath(path, "london.csv")) do file
        CSV.read(file,DataFrame,
                 dateformat = "yyyy-mm-dd",
                 types = Dict(:Date => Date))
    end

end


text:
# IRI_assF_kfold_include.jl

include("../tcm_code/tcm_transactionsToCounts.jl")

include("../tcm_code/tcm_estimateRanking_CS.jl")
include("../tcm_code/tcm_predictRanking.jl")

include("../tcm_code/tcm_evaluateKL.jl")

include("../tcm_code/tcm_rankingStats.jl")

include("R2_IRI_temporalRCS_helperfn.jl")

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# This file contains the two PV cell models (SDM and DDM)
# and functions to compute the SSE as well as the RMSE.

const k = 1.3806503e-23     # Boltamann constant
const q = 1.60217646e-19    # electron charge

"""
    Single diode model with five parameters

V, I are the input voltage and current.

Ir:  photocurrent generated by the solar cell (A)
I0: the reverse saturation current of the diode (μA)
n: diode ideality factor
Rs: series resistance (Ω)
Rp: parallel resistance (Ω)
"""
function model_sdm(V, I, Iph, I0, n, Rs, Rp, T)
    tmp = V+ I*Rs
    return Iph - I0*1e-6*(exp(q*tmp/(n*k*T)) - 1) - tmp/Rp
end

"""
    Double diode model with seven parameters

V, I are the input voltage and current.

Iph:  photocurrent generated by the solar cell (A)
I01, I02: the reverse saturation current of the diode (μA)
n1, n2: diode ideality factor
Rs: series resistance (Ω)
Rp: parallel resistance (Ω)
"""
function model_ddm(V, I, Iph, I01, I02, n1, n2, Rs, Rp, T)
    tmp = V+ I*Rs
    c = k*T/q
    return Iph - I01*1e-6*(exp(tmp/(n1*c)) - 1) - I02*1e-6*(exp(tmp/(n2*c)) - 1) - tmp/Rp
end

"""
	Compute the sum of squared error (SSE)

data: a matrix of V-I data, each row denoting a (V, I) point
model: a function for the SDM or the DDM, whose input is (V, I) and output is the estimated current
"""
function calculate_sse(data, model::Function)
    sse = 0.0
    for (Vm, Im) in eachrow(data)
        Ic = model(Vm, Im)
        sse += (Im - Ic)^2
    end
    return sse
end

"""
	Compute the root mean square error (RMSE)

data: a matrix of V-I data, each row denoting a (V, I) point
model: a function for the SDM or the DDM, whose input is (V, I) and output is the estimated current
"""
function calculate_rmse(data, model::Function)
    sse = calculate_sse(data, model)
    return sqrt(sse / size(data, 1))
end

"""
    Compute the root mean square error (RMSE)

data: a matrix of V-I data
θ: parameter vector. Depending on the length of θ, either SDM or DDM will be used.
T: temperature in Calvin (K)
"""
function calculate_rmse(data, θ, T)
    model = length(θ) == 5 ? model_sdm : model_ddm
    return calculate_rmse(data, (V, I)->model(V, I, θ..., T))
end

"""
	Compute the mean bias error (MBE)

data: a matrix of V-I data, each row denoting a (V, I) point
model: a function for the SDM or the DDM, whose input is (V, I) and output is the estimated current
"""
function calculate_mbe(data, model::Function)
    mbe = 0.0
    for (Vm, Im) in eachrow(data)
        Ic = model(Vm, Im)
        mbe += (Im - Ic)
    end
    return mbe / size(data, 1)
end

"""
	Compute the mean absolute percentage error (MAPE).

See https://en.wikipedia.org/wiki/Mean_absolute_percentage_error for a definition.

data: a matrix of V-I data, each row denoting a (V, I) point
model: a function for the SDM or the DDM, whose input is (V, I) and output is the estimated current
"""
function calculate_mape(data, model::Function)
    mape = 0.0
    for (Vm, Im) in eachrow(data)
        Ic = model(Vm, Im)
        d = abs((Im - Ic) / Im)
        # @show d
        mape += d
    end
    return mape / size(data, 1) * 100
end

"""
	Compute the mean absolute error (MAE)

data: a matrix of V-I data, each row denoting a (V, I) point
model: a function for the SDM or the DDM, whose input is (V, I) and output is the estimated current
"""
function calculate_mae(data, model::Function)
    mae = 0.0
    for (Vm, Im) in eachrow(data)
        Ic = model(Vm, Im)
        mae += abs(Im - Ic)
    end
    return mae / size(data, 1)
end

max_stars_repo_name: ShuhuaGao/rePVest

max_issues_repo_name: ShuhuaGao/rePVest

max_forks_repo_name: ShuhuaGao/rePVest


text:
# built-in: true, false and missing

using Printf

const tril = (true, missing, false)

@printf("\n%8s | %8s\n", "A", "¬A")
for A in tril
    @printf("%8s | %8s\n", A, !A)
end

@printf("\n%8s | %8s | %8s\n", "A", "B", "A ∧ B")
for (A, B) in Iterators.product(tril, tril)
    @printf("%8s | %8s | %8s\n", A, B, A & B)
end

@printf("\n%8s | %8s | %8s\n", "A", "B", "A ∨ B")
for (A, B) in Iterators.product(tril, tril)
    @printf("%8s | %8s | %8s\n", A, B, A | B)
end

@printf("\n%8s | %8s | %8s\n", "A", "B", "A ≡ B")
for (A, B) in Iterators.product(tril, tril)
    @printf("%8s | %8s | %8s\n", A, B, A == B)
end

⊃(A, B) = B | !A

@printf("\n%8s | %8s | %8s\n", "A", "B", "A ⊃ B")
for (A, B) in Iterators.product(tril, tril)
    @printf("%8s | %8s | %8s\n", A, B, A ⊃ B)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Gadfly
using ANOVADDPTest
using DataFrames
using Distributions
using Random

function generate_plot_bernoulli(N, G, signigicant_groups)
    rng = MersenneTwister(1)
    m = BernoulliDDP(rng, N, G)
    x = rand(rng, 1:G, N)
    y = rand(rng, Bernoulli(0.25), N)
    for i = 1:N
        x[i] in signigicant_groups && (y[i] = rand(rng, Bernoulli(0.75)))
    end
    data = BernoulliData(x, y)
    pred = BernoulliData(expandgrid(1:G, false:true)...)
    ch = train(rng, m, data, pred)
    df = DataFrame(x = pred.x, y = pred.y, f = mean(ch.f))
    plot(df, x = :y, y = :f, color = :x, Geom.line, Scale.color_discrete())
end

generate_plot_bernoulli(1000, 3, [])
generate_plot_bernoulli(1000, 3, [2])
generate_plot_bernoulli(1000, 3, [3])
generate_plot_bernoulli(1000, 3, [2, 3])


text:
################################################################################
#
#  Basic field access
#
################################################################################

@doc Markdown.doc"""
    nf(O::NumFieldOrd) -> NumField

Returns the ambient number field of $\mathcal O$.
"""
nf(O::NumFieldOrd)

@inline nf(O::NfAbsOrd{S, T}) where {S, T} = O.nf::S

@inline nf(O::NfRelOrd{S, T, U}) where {S, T, U} = O.nf::parent_type(U)

_algebra(O::NumFieldOrd) = nf(O)

@doc Markdown.doc"""
    number_field(O::NumFieldOrd)

Return the ambient number field of $\mathcal O$.
"""
@inline function NumberField(O::NumFieldOrd)
  return nf(O)
end

is_simple(O::NumFieldOrd) = is_simple(nf(O))

is_commutative(O::NumFieldOrd) = true

@doc Markdown.doc"""
    is_equation_order(O::NumFieldOrd) -> Bool

Returns whether $\mathcal O$ is the equation order of the ambient number
field $K$.
"""
@inline is_equation_order(O::NumFieldOrd) = O.is_equation_order

@inline is_maximal_known(O::NumFieldOrd) = O.is_maximal != 0

@inline is_maximal_known_and_maximal(O::NumFieldOrd) = isone(O.is_maximal)

@doc Markdown.doc"""
    is_maximal(R::NfAbsOrd) -> Bool

Tests if the order $R$ is maximal. This might trigger the
computation of the maximal order.
"""
function is_maximal(R::NumFieldOrd)
  if R.is_maximal == 1
    return true
  end
  if R.is_maximal == 2
    return false
  end
  S = maximal_order(R)
  if discriminant(S) == discriminant(R)
    R.is_maximal = 1
  else
    R.is_maximal = 2
  end
  return R.is_maximal == 1
end

################################################################################
#
#  Degree
#
################################################################################

@doc Markdown.doc"""
    degree(O::NumFieldOrd) -> Int

Returns the degree of $\mathcal O$.
"""
degree(O::NumFieldOrd) = degree(O.nf)

function absolute_degree(O::NumFieldOrd)
  return absolute_degree(nf(O))
end

################################################################################
#
#  Signature
#
################################################################################

@doc Markdown.doc"""
    signature(O::NumFieldOrd) -> Tuple{Int, Int}

Returns the signature of the ambient number field of $\mathcal O$.
"""
function signature(x::NumFieldOrd)
  return signature(nf(x))
end

################################################################################
#
#  Discriminant
#
################################################################################

function discriminant(O::NumFieldOrd, K::NumField)
  K = nf(O)
  if K == base_field(K)
    return discriminant(O)
  else
    return norm(discriminant(O), K)*discriminant(base_ring(O), K)^degree(O)
  end
end

discriminant(O::NumFieldOrd, ::FlintRationalField) = absolute_discriminant(O)

################################################################################
#
#   Absolute basis
#
################################################################################

function absolute_basis(O::NfAbsOrd)
  return basis(O)
end

function absolute_basis(O::NfRelOrd)
  pb = pseudo_basis(O, copy = false)
  res = Vector{elem_type(O)}()
  for i = 1:degree(O)
    for b in absolute_basis(pb[i][2])
      push!(res, O(b*pb[i][1]))
    end
  end
  return res
end


text:
module test_julia_task

using Test

t = @async begin
    return true
end
@test fetch(t)

t = @async begin
    return "hello"
end
@test fetch(t) == "hello"

t = @async "world"
@test fetch(t) == "world"

@test endswith(task_local_storage(:SOURCE_PATH), "task.jl")

end # module test_julia_task


text:
using JSON
using HTTP
using Test

@testset "HTTP.URL" begin
    tests = JSON.parse(String(read("resources/cweb-urls.json")))["tests"]["group"]

    @testset " - $test - $group" for group in tests, test in group["test"]
        name = group["name"]

        url = get(test, "url", nothing)

        if url !== nothing
            uri = nothing

            try
                uri = HTTP.URIs.parse_uri_reference(url; strict=true)
            catch e
                if e isa HTTP.URIs.ParseError || e isa AssertionError
                    continue
                else
                    rethrow(e)
                end
            end

            if haskey(test, "expect_protocol")
                @test uri.scheme == test["expect_protocol"][1:end-1]
            end

            if haskey(test, "expect_hostname")
                @test uri.host == test["expect_hostname"]
            end

            if haskey(test, "expect_port")
                @test uri.port == test["expect_port"]
            end
        end
    end
end


text: # Welcome to this introduction!


text:
import Base: iterate, length, eachcol

struct JDFFileColIterator
    jdf::JDFFile
    cols::Vector{Union{String,Symbol}}
end

eachcol(jdf::JDFFile) = JDFFileColIterator(jdf, names(jdf))

Base.length(jdf::JDFFileColIterator) = length(jdf.cols)

Base.iterate(jdf::JDFFileColIterator, state = 1) = begin
    if state > length(jdf.cols)
        return nothing
    end

    res1 = sload(jdf.jdf, cols = [jdf.cols[state]])[1]
    return (res1, state + 1)
end


text:
module DuckDB
using DataFrames, Dates, DuckDB_jll
using DBInterface, Tables

export DBInterface

include("api.jl")
include("consts.jl")

const DBHandle = Ref{Ptr{Cvoid}}   # DuckDB DB connection handle

"""
    toDataFrame(connection::Ref{Ptr{Cvoid}},query::String)::DataFrame

Creates a DataFrame from a SQL query within a connection.

* `connection`: The connection to perform the query in.
* `query`: The SQL query to run.
* returns: the abstract dataframe

"""
function toDataFrame(connection::Ref{Ptr{Cvoid}}, query::String)::DataFrame
    res = execute(connection, query)::Ref{duckdb_result}
    return toDataFrame(res)
end
"""
    toDataFrame(result::Ref{duckdb_result})::DataFrame

Creates a DataFrame from the full result
* `result`: the full result from `execute`
* returns: the abstract dataframe

"""
function toDataFrame(result::Ref{duckdb_result})::DataFrame
    columns =
        unsafe_wrap(Array{duckdb_column}, result[].columns, Int64(result[].column_count))
    df = DataFrame()
    for i = 1:Int64(result[].column_count)
        rows = Int64(result[].row_count)
        name = (unsafe_string(columns[i].name))
        type = DUCKDB_TYPE(Int64(columns[i].type))
        if type == DUCKDB_TYPE_INVALID
            print("invalid type for column - \"" * name * "\"")
        else
            mask = unsafe_wrap(Array, columns[i].nullmask, rows)
            data = unsafe_wrap(Array, Ptr{DUCKDB_TYPES[type]}(columns[i].data), rows)
            bmask = reinterpret(Bool, mask)

            if 0 != sum(mask)
                data = data[.!bmask]
            end
            if type == DUCKDB_TYPE_DATE
                data = Dates.epochdays2date.(data .+ 719528)
            elseif type == DUCKDB_TYPE_TIME
                data = Dates.Time.(Dates.Nanosecond.(data .* 1000))
            elseif type == DUCKDB_TYPE_TIMESTAMP
                data = Dates.epochms2datetime.((data ./ 1000) .+ 62167219200000)
            elseif type == DUCKDB_TYPE_INTERVAL
                data = map(
                    x -> Dates.CompoundPeriod(
                        Dates.Month(x.months),
                        Dates.Day(x.days),
                        Dates.Microsecond(x.micros),
                    ),
                    data,
                )
            elseif type == DUCKDB_TYPE_HUGEINT
                data = map(x -> x.upper < 1 ? (x.lower::UInt64)%Int64 : x, data)
            elseif type == DUCKDB_TYPE_VARCHAR
                data = unsafe_string.(data)
            end

            if 0 != sum(mask)
                fulldata = Array{Union{Missing,eltype(data)}}(missing, rows)
                fulldata[.!bmask] = data
                data = fulldata
            end

            df[!, name] = data
        end
    end
    return df
end
"""
    disconnect(connection)
Closes the specified connection and de-allocates all memory allocated for that connection.
* `connection`: The connection to close.

"""
function disconnect(connection)
    return duckdb_disconnect(connection)
end

"""
    close(database)
Closes the specified database and de-allocates all memory allocated for that database.\n
This should be called after you are done with any database allocated through duckdb_open.\n
Note that failing to call duckdb_close (in case of e.g. a program crash) will not cause data corruption. Still it is recommended to always correctly close a database object after you are done with it.
*`database`: the database object to shut down.

"""
function close(database)
    return duckdb_close(database)
end

connect() = connect(":memory:")

"""
    connect(path)
Creates a new database or opens an existing database file stored at the the given path. If no path is given a new in-memory database is created instead.
* `path`: Path to the database file on disk or `:memory:` to open an in-memory database.
* returns: a connection handle

"""
function connect(path::String)::Ref{Ptr{Cvoid}}
    database = Ref{Ptr{Cvoid}}()
    connection = Ref{Ptr{Cvoid}}()
    duckdb_open(path, database)
    duckdb_connect(database, connection)
    return connection
end

"""
    execute(connection, query)

Executes a SQL query within a connection and returns the full (materialized) result. If the query fails to execute, `DuckDBError` is returned and the error message can be retrieved by calling `duckdb_result_error`.

Note that after running duckdb_query, duckdb_destroy_result must be called on the result object even if the query fails, otherwise the error stored within the result will not be freed correctly.
* `connection`: The connection to perform the query in.
* `query`: The SQL query to run.
* returns: the full result pointer

"""
function execute(connection::Ref{Ptr{Cvoid}}, query::String)::Ref{duckdb_result}
    result = Ref{duckdb_result}()
    duckdb_query(connection, query, result)
    if result[].error_message != Ptr{UInt8}(0)
        print(unsafe_string(result[].error_message))
    end
    return result
end

struct DuckDBException <: Exception
    msg::AbstractString
end

duckdb_errmsg(handle) = unsafe_string(handle[].error_message)
duckdbexception(handle::DBHandle) = DuckDBException(duckdb_errmsg(handle))
duckdberror(handle::DBHandle) = throw(duckdbexception(handle))

"""
    `DuckDB.DB()` => in-memory DuckDB database
    `DuckDB.DB(file)` => file-based DuckDB database

Constructors for a representation of a DuckDB database, either backed by an on-disk file or in-memory.

`DuckDB.DB` requires the `file` string argument in the 2nd definition
as the name of either a pre-defined DuckDB database to be opened,
or if the file doesn't exist, a database will be created.

The `DuckDB.DB` object represents a single connection to a DuckDB database.
All other DuckDB.jl functions take an `DuckDB.DB` as the first argument as context.

To create an in-memory temporary database, call `DuckDB.DB()`.

The `DuckDB.DB` will be automatically closed/shutdown when it goes out of scope
(i.e. the end of the Julia session, end of a function call wherein it was created, etc.)

NOTE: This borrows heavily from SQLite.jl, [here](https://github.com/JuliaDatabases/SQLite.jl/blob/9724a175ae30b22cb01775250ce2e87317706215/src/SQLite.jl)
"""
mutable struct DB <: DBInterface.Connection
    file::String
    handle::DBHandle

    # TODO: Support prepared statements

    function DB(f::AbstractString)
        f = String(isempty(f) ? f : expanduser(f))

        try
            handle = connect(f)

            db = new(f, handle)
            finalizer(_close, db)
            return db
        catch
            duckdberror(handle)
        end
    end
end

function _close(db::DB)
    db.handle == C_NULL || duckdb_disconnect(db.handle)
    db.handle = C_NULL
    return
end

DB() = DB(":memory:")
DBInterface.close!(db::DB) = _close(db)
Base.close(db::DB) = _close(db)

Base.show(io::IO, db::DB) = print(io, string("DB(", "\"$(db.file)\"", ")"))

"""
    DBInterface.execute(db::DB, sql::String, [params])

TODO: Support DBInterface statements

Take inputs given by `db` and `sql` execute the query and return an iterator of result rows.

Note that the returned result row iterator only supports a single-pass, forward-only iteration of the result rows.

TODO: Support Tables.jl
"""
function DBInterface.execute(db::DB, sql::String)

    result = execute(db.handle, sql)

    return result
end

end # module


text:
push!(LOAD_PATH,".")

using TestManager, TestSpecDef, JuMP, JLD, Dictionaries, NamedDims

maxSeconds = 30
maxSamples = 1000

testSpec1 = TestSpec(A,1,maxSeconds,maxSamples,"Data/arrayData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D10",1)
testSpec2 = TestSpec(A,2,maxSeconds,maxSamples,"Data/arrayData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F10",1)

testSpec3 = TestSpec(AL,1,maxSeconds,maxSamples,"Data/arrayData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D16",1)
testSpec4 = TestSpec(AL,2,maxSeconds,maxSamples,"Data/arrayData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F16",1)

testSpec5 = TestSpec(ALB,1,maxSeconds,maxSamples,"Data/arrayData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D22",1)
testSpec6 = TestSpec(ALB,2,maxSeconds,maxSamples,"Data/arrayData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F22",1)

testSpec7 = TestSpec(DAA,1,maxSeconds,maxSamples,"Data/denseAxisArrayData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D28",1)
testSpec8 = TestSpec(DAA,2,maxSeconds,maxSamples,"Data/denseAxisArrayData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F28",1)

testSpec9 = TestSpec(DAAL,1,maxSeconds,maxSamples,"Data/denseAxisArrayData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D34",1)
testSpec10 = TestSpec(DAAL,2,maxSeconds,maxSamples,"Data/denseAxisArrayData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F34",1)

testSpec11 = TestSpec(DAALS,1,maxSeconds,maxSamples,"Data/denseAxisArrayData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D40",1)
testSpec12 = TestSpec(DAALS,2,maxSeconds,maxSamples,"Data/denseAxisArrayData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F40",1)

testSpec13 = TestSpec(DAALB,1,maxSeconds,maxSamples,"Data/denseAxisArrayData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D46",1)
testSpec14 = TestSpec(DAALB,2,maxSeconds,maxSamples,"Data/denseAxisArrayData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F46",1)

testSpec15 = TestSpec(DAAS20,1,maxSeconds,maxSamples,"Data/denseAxisArrayDataStr1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D52",1)
testSpec16 = TestSpec(DAAS20,2,maxSeconds,maxSamples,"Data/denseAxisArrayDataStr1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F52",1)

testSpec17 = TestSpec(D,1,maxSeconds,maxSamples,"Data/dictData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D58",1)
testSpec18 = TestSpec(D,2,maxSeconds,maxSamples,"Data/dictData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F58",1)

testSpec19 = TestSpec(DD,1,maxSeconds,maxSamples,"Data/dictionaryData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D64",1)
testSpec20 = TestSpec(DD,2,maxSeconds,maxSamples,"Data/dictionaryData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F64",1)

testSpec21 = TestSpec(ND,1,maxSeconds,maxSamples,"Data/namedDimsData1000.jld",false,false,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","D70",1)
testSpec22 = TestSpec(ND,2,maxSeconds,maxSamples,"Data/namedDimsData1000.jld",false,true,true,true,"Results/Excel/DatatypeBenchmarks2.xlsx","F70",1)

#testSpecs = [testSpec1 testSpec2 testSpec3 testSpec4 testSpec5 testSpec6 testSpec7 testSpec8 testSpec9 testSpec10 testSpec11 testSpec12 testSpec13 testSpec14 testSpec15 testSpec16 testSpec17 testSpec18 testSpec19 testSpec20 testSpec21 testSpec22]
testSpecs = [testSpec11 testSpec12 testSpec13 testSpec14 testSpec15 testSpec16 testSpec17 testSpec18 testSpec19 testSpec20 testSpec21 testSpec22]
for testSpec in testSpecs
   Test(testSpec)
end

max_stars_repo_name: SonntagK/JuMPContainerBenchmarks

max_issues_repo_name: SonntagK/JuMPContainerBenchmarks

max_forks_repo_name: SonntagK/JuMPContainerBenchmarks


text:
"""
$(SIGNATURES)

Given the raw data from [`Simulation.sequencing`](@ref) returns two DataFrames

1. `guide_data`: This DataFrame contains the per-guide level data including the
    log2 fold change in the normalized frequencies of each guide between each
    pairwise combination of bins. Thus, if there are `n` bins, then it computes
    the log2 fold changes for the ``\\frac{n!}{2(n-2)!}`` combinations

2. `gene_data`: This DataFrame contains the same information but grouped by
    gene. The log2 fold change data from the first DataFrame is used to calculate
    the average log2 fold change per gene and a pvalue computed using a
    [Mann-Whitney U-test](https://en.wikipedia.org/wiki/Mann-Whitney_U_test) as
    measure of how consistently shifted the guides are of this gene versus the
    population of negative control guides. (see below for more info)

A typical 2 bin `guide_data` DataFrame contains the following columns:

| Column Name  | Meaning  |
|:--|:--|
| `gene` | the gene ID of that this guide targets|
| `knockdown` | activity of the guide on 0 to 1 scale, where 1 is complete knockout|
| `barcodeid` | the ID of this specific guide|
| `theo_phenotype` | expected phenotype of this guide, generally a -1 to 1 scale|
| `behavior` | whether the target gene displays a linear or sigmoidal response to incomplete knockdown (see [`Simulation.Library`](@ref) for more details)|
| `class` | which phenotype distribution the target gene was drawn from (see [`Simulation.Library`](@ref) for more details). Serves as the "ground truth" label against which screen performance is evaluated, e.g. with [`Simulation.auprc`](@ref) |
| `initial_freq` | frequency of guide post-transfection (see [`Simulation.transfect`](@ref))|
| `counts_bin1` | the raw number of reads for each guide in the first bin|
| `freqs_bin1` | the number of reads for each guide divided by the total number of reads in this bin|
| `rel_freqs_bin1` | the frequency of each guide divided by the median frequency of negative control guides|
| `counts_bin2` | the raw number of reads for each guide in the second bin|
| `freqs_bin2` | the number of reads for each guide divided by the total number of reads in this bin|
| `rel_freqs_bin2` | the frequency of each guide divided by the median frequency of negative control guides for this bin|
| `log2fc_bin2_div_bin1` | the log2 fold change in relative guide frequencies between `bin2` and `bin1`, equivalent to `log2(rel_freqs_bin2/rel_freqs_bin1)` |

A typical `gene_data` DataFrame contains the following data:

| Column Name  | Meaning  |
|:--|:--|
| `gene` | this gene's ID|
| `class` | see above |
| `behavior` | see above |
| `mean_bin2_div_bin1` | the mean log 2 fold change in relative frequencies between from `bin1` to `bin2` for all the guides targeting this gene. Calculated as ``\\frac{1}{k}\\sum_k \\text{log2fc_bin2_div_bin1}_k`` for the ``k`` guides targeting each gene |
| `pvalue_bin2_div_bin1` | the -log10 pvalue of the log2 fold changes of all guides targeting this gene as computed by the non-parametric Mann-Whitney U-test. A measure of the consistency of the log 2 fold changes[^1]|
| `absmean_bin2_div_bin1` | absolute value of `mean_bin2_div_bin1` per-gene |
| `pvalmeanprod_bin2_div_bin1` | `mean_bin2_div_bin1` multiplied with the `pvalue_bin2_div_bin1` per-gene|

# Further reading

[^1]: Kampmann M, Bassik MC, Weissman JS. Integrated platform for genome-wide screening
    and construction of high-density genetic interaction maps in mammalian cells.
    *Proc Natl Acad Sci U S A*. 2013;110:E2317–26.

"""
function differences_between_bins(raw_data::Associative{Symbol, DataFrame})
    for (bin, seq_data) in raw_data
        sort!(seq_data, [:barcodeid])
        # add a pseudocount of 0.5 to every value to prevent -Inf's when
        # taking the log
        seq_data[:counts] += 0.5
        seq_data[:freqs] = seq_data[:counts]./sum(seq_data[:counts])
        # normalize to median of negative controls, fixes #19
        # TODO: consider normalizing the std dev
        negcontrol_freqs = seq_data[seq_data[:class] .== :negcontrol, :freqs]
        (length(negcontrol_freqs) == 0) && error("No negative control guides found. Try increasing "*
            "the frequency of negative controls or increase the number of genes.")
        med = median(negcontrol_freqs)
        seq_data[:rel_freqs] = seq_data[:freqs] ./ med
    end

    bin1 = first(raw_data)[2]
    cols_to_copy = [col for col in names(bin1) if !(col in (:counts, :freqs, :rel_freqs))]
    guide_data = bin1[cols_to_copy]
    gene_data = DataFrame()

    # processed bins; so we don't re-add bins when doing all combos
    proc_bins = Set{Symbol}()

    # compute pairwise log 2 fold changes between every combination of bins
    for bins in combinations(collect(keys(raw_data)), 2)
        for bin in bins
            (bin in proc_bins) && continue
            push!(proc_bins, bin)
            guide_data[[Symbol("counts_", bin), Symbol("freqs_", bin), Symbol("rel_freqs_", bin)]] =
                raw_data[bin][[:counts, :freqs, :rel_freqs]]
        end
        suffix = "_$(bins[2])_div_$(bins[1])"
        log2fc_col = Symbol(:log2fc, suffix)
        guide_data[log2fc_col] = log2.(guide_data[Symbol(:rel_freqs_, bins[2])]
                                        ./guide_data[Symbol(:rel_freqs_, bins[1])])

        nonnegs = guide_data[guide_data[:class] .!= :negcontrol, :]
        negcontrols = guide_data[guide_data[:class] .== :negcontrol, log2fc_col]

        tmp = by(nonnegs, [:gene, :behavior, :class]) do barcodes
            log2fcs = barcodes[log2fc_col]
            result = MannWhitneyUTest(log2fcs, negcontrols)
            DataFrame(pvalue = -log10(pvalue(result)), mean= mean(log2fcs))
        end
        gene_data[[:gene, :behavior, :class, Symbol(:pvalue, suffix), Symbol(:mean, suffix)]] = tmp[[:gene, :behavior, :class, :pvalue, :mean]]
        gene_data[Symbol(:absmean, suffix)] = abs.(gene_data[Symbol(:mean, suffix)])
        gene_data[Symbol(:pvalmeanprod, suffix)] = gene_data[Symbol(:mean, suffix)] .* gene_data[Symbol(:pvalue, suffix)]

    end

    guide_data, gene_data
end


text:
function makewave(f::Int64, t::Float32; amplitude::Float32=1, fs::Int64=44100)
    # f : 频率，单位 HZ
    # t : 时长，单位秒
    # amplitude: 振幅，默认为 1
    # Fs: 采样频率，默认为 44100

    N = ceil(Int64, Fs*t) # 总的采样点个数
    L = floor(Int64, Fs/f) # 一个周期所含的采样数(整)
    part = rand(L, 1) - 0.5 # 循环段
    np = floor(Int64, N/L) # 总的循环次数

    r = 0.5; # 平均权值

    wave=zeros(N, 1)
    for i in 1:np
        part = part .* r + [part[end]; part[1:end-1]] .* (1 - r)   # 吉他音质算法
        pos = L * (i -1) + 1
        wave[pos:pos+L-1] = part
    end

    wave = wave .* range(1, stop=0, length=length(wave))  # 归一化

    wave * amplitude
end

# 需要用到第三方模块：WAV
using WAV

Fs = 44100
a1 = 440;  #标准音
r = 2^(1/12); #12平均律的比值
t = .5; #1拍的时长，半秒

#主旋律 〈八辈子〉
song = [3 1; 5 1; 3 1; 3 .5; 5 .5; 10 1; 10 1; 10 .5; 7 .5; 5 .5; 7 .5; 5 1; 3 1; 3 1
   3 .5; 5 .5; 5 1; 3 1; 3 1; -2 .5; -2 .5]

#伴奏    〈八辈子〉
chords = [3 -2 3 5 7 5 10 7 3 -2 3 5 7 5 10 7 -2 -2 3 5 7 5 10 7 -2 -2 3 5 7 5 7 10]'
t_chord = 0.5  # 所有伴奏音符的时长都为半拍

beats = sum(song[:,2])

l = Int(beats * 2); # 伴奏每个音符都只有半拍，所以长度乘个二
while length(chords) < l
    chords = [chords; chords]
end
chords = chords[1:l]

wave =[]
for i = 1:size(song, 1)
    f = a1 * r^song[i,1]
    f_t = t * song[i,2]

    tone = makewave(f, f_t, 1, Fs)
    wave = isempty(wave) ? tone : [wave; tone]
end

wave_c = []
for i = 1:size(chords, 1)
    f = a1 * r^chords[i]
    f_t = t * .5;  # 伴奏每个音符都只有半拍

    tone = makewave(f, f_t, .3, Fs)
    wave_c = isempty(wave_c) ? tone : [wave_c; tone]
end

ll = min(length(wave), length(wave_c))

wav = wave[1:ll] + wave_c[1:ll]
wavplay(wav, Fs)

max_stars_repo_name: yuansuye/generate-song

max_issues_repo_name: yuansuye/generate-song

max_forks_repo_name: yuansuye/generate-song