text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0101 - Basic Syntax Elements

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
ENV["PYTHON"] = ""  # This will be needed in a further segment
ENV["R_HOME"] = "*" # This wil lbe needed in a further segment
## If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will hower be lost)
## Pkg.resolve()
Pkg.instantiate()
using Random
Random.seed!(123)

# ## Comments

## This is a comment
a = 1 # also this one
a = #= also this one =# 1
#= also
#= this =#
one
=#

# ## Code organisation

## Semicolon:
a = 1
a = 1;
for i in 1:3
   println("i is $i")
end # Keyword `end` to finish a block
println("Hello world!")
## println("This would error!")

# ## Unicode support

# We can can use any fancy Unicode symbol, including modifiers, in the names of variables, types, functions..
using Statistics # for the `mean` function, in the Standard Library
σ²(x) = sum( (x .- mean(x)).^2 )/length(x) # The superscript `²` is just another character, it has no syntactic value.
σ²([1,2,3])
x̄ₙ = 10
vàlidVarNαme! = 2
🥞 = 8
🍴(🥫,parts=4) = 🥫/parts
🍰 = 🍴(🥞)

# ## Broadcasting

# _Broadcasting_ refers to the capacity to dispatch a function that accept a scalar as argument to any collection for that argument.
# It is implemented by postfixing a function name (or prefixing an operator) with a single dot at call time:

10 .+ [1,2,3]
add2(x) = x + 2
add2(10)
## add2([1,2,3]) # would return an error
add2.([1,2,3])  # any, including user defined functions, can be broadcasted. No need for map, for loops, etc..
add2.((1,2,3)) # not only for arrays, here is a Tuple
add2.(Set([1,2,3,2])) # and here a Set
.+([1,2,3],[10,20,30]) # fine here
## .+([1,2,3],[10,20]) # DimensionMismatch error: the input of the the broadcasted arguments must have the same size or be a scalar
# To "protect" one specific argument to be broadcasted, use `Ref(arg)`:
foo(x,y::AbstractArray) = [yi+x for yi in y] # a function that wants the first arguent as scalar and the second as vector
foo(1,[10,20])
## foo.([1,2],[10,20])  # error, as we try to broadcast also the seond element, but we can't call `foo` with two integers
foo.([1,2],Ref([10,20])) # now it is fine, we broadcast only the first argument

# ## 1 based arrays
# Arrays start indexing at 1 (like R, Fortran, Matlab,...) instead of 0 (like in C, Python, Java...)
a = [1,2,3]
a[1]
collect(1:3) # in ranges, both extremes are included

# !!! joke "0 or 1 ?"
#     Should array indices start at 0 or 1?  My compromise of 0.5 was rejected without, I thought, proper consideration. --Stan Kelly-Bootle

# ## Basic Mathematic operations

# All standard mathematical arithmetic operators (`+`,`-`,`*`,`/`) are supported in the obvious way:

a = 2^4         # rise to power
b = ℯ^2; #= or =# b = exp(2) # Exponential with base ℯ
d = log(7.3890) # base ℯ
e = log(10,100) # custom base
f = 5 ÷ 2       # integer division
e = 5 % 2       # reminder (modulo operator)
a = 2//3 + 1//3     # rational numbers
typeof(a)
π == pi         # some irrational constants
typeof(ℯ)
convert(Float64,a)

# ## Quotation

a = 'k'              # single quotation mark: a single Char
b = "k"              # double quotation mark: a (Unicode) String
#c = 'hello'         # error !
c = "hello"
d = `echo hello`     # backtick: define a command to (later) run (e.g. on the OS)
e = """a
multiline
string"""
println(c)
println(e)
using Markdown
f = md"""a
**markdown**
_string_"""

# ## Missingness implementations

# Attention to these three different implementations (of 3 different concepts) of "missingness":

a = nothing # C-style, "software engineer's null → run-time error
b = missing # Data scientist's null → silent propagation
c = NaN     # Not a number → silent propagation
typeof(a)
typeof(b)
typeof(c)
d = 0/0
## a2 = mean([1,a,3]) # would error
b2 = mean([1,b,3])
c2 = mean([1,c,3])
b3 = mean(skipmissing([1,b,3]))
b == missing # propagate
ismissing(b)
isequal(b,2)  # exception to the propagation rule, it allows comparations when one of the item may be missing
isequal(b,missing)
b4 = [1,missing,3]
typeof(b4)
eltype(b4)
nonmissingtype(eltype(b4))

# ## Random values

rand()       # [0,1] continuous
rand(30:40)  # [30,40] integer
rand(30:0.01:40) # [30,40] with precision to the second digit
using Distributions
rand(Exponential(10)) # We'll see Distributions more in detail in the Scientific Programming lesson
rand(30:40,3) # A vector of 3 random numbers.
rand(Exponential(10),3,4) # Same syntax for sampling from any distribution !
using Random
myRNG = MersenneTwister(123) # use StableRNG for a RNG guaranteed to remain stable between Julia-versions
a1 = rand(myRNG,10:1000,5)
a2 = rand(myRNG,10:1000,5)
a1 == a2
myRNG = MersenneTwister(123)
b1 = rand(myRNG,10:1000,5)
b2 = rand(myRNG,10:1000,5)
b1 == b2
a1 == b1
a2 == b2

a = rand(myRNG,Exponential(10),5)

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
module UniDerivativeTests

import ..Algopt # ensure we are using the correct Algopt

using Test

using Algopt.UniDerivative:
df,
CentralDiff, ComplexDiff

quadratic = x->x^2

@testset "Algopt.UniDerivative.df(::CentralDiff)" begin
    @test df(CentralDiff(), quadratic, 3) ≈ 2 * 3
end

@testset "Algopt.UniDerivative.df(::ComplexDiff)" begin
    @test df(ComplexDiff(), quadratic, 3) ≈ 2 * 3
end

end # module

ext: jl


text:
module Drums
using MusicManipulations
using MIDI: NOTEON, NOTEOFF
# all pitches of digitally connected instruments that support extra velocities
const DIGITAL = [0x26,0x28,0x33,0x35,0x3b]

"""
    getnotes_td50(track::MIDITrack, tpq = 960)

Get notes from midi track. Take care of Roland TD-50's ability to have velocities up to 159 for snare and ride.
"""
getnotes_td50(midi::MIDIFile, i::Int) = getnotes_td50(midi.tracks[i], midi.tpq)
function getnotes_td50(track::MIDI.MIDITrack, tpq = 960)
    notes = MoreVelNote[]
    tracktime = UInt(0)
    extravel = 0
    for (i, event) in enumerate(track.events)
        tracktime += event.dT
        # Read through events until a noteon with velocity higher tha 0 is found
        if isa(event, MIDIEvent) && event.status & 0xF0 == NOTEON && event.data[2] > 0
            duration = UInt(0)
            #Test if the next event is an extra velocity event and modify velocity if needed.
            if event.data[1] in DIGITAL && event.data[2]==0x7f && track.events[i+1].status==0xb0 && track.events[i+1].data[1]==0x58
                extravel = floor(UInt8,track.events[i+1].data[2]/2)
                if extravel > 32
                    extravel = 32
                end
            end
            #Test if the previous event is an extra velocity event and modify velocity if needed.
            if i>2 #first event is alwas METAEvent
                if event.data[1] in DIGITAL && event.data[2]==0x7f && track.events[i-1].status==0xb0 && track.events[i-1].data[1]==0x58
                    extravel = floor(UInt8,track.events[i-1].data[2]/2)
                    if extravel > 32
                        extravel = 32
                    end
                end
            end
            for event2 in track.events[i+1:length(track.events)]
                duration += event2.dT
                # If we have a MIDI event & it's a noteoff (or a note on with 0 velocity), and it's for the same note as the first event we found, make a note
                # Many MIDI files will encode note offs as note ons with velocity zero
                if isa(event2, MIDI.MIDIEvent) && (event2.status & 0xF0 == MIDI.NOTEOFF || (event2.status & 0xF0 == MIDI.NOTEON && event2.data[2] == 0)) && event.data[1] == event2.data[1]
                    push!(notes, MoreVelNote(event.data[1], event.data[2]+extravel, tracktime, duration, event.status & 0x0F))
                    break
                end
            end
            extravel = 0
        end
    end
    sort!(notes, lt=((x, y)->x.position<y.position))
    return Notes(notes, tpq)
end

"""
    rm_hihatfake(notes;
                 back = 100, forw = 100, cutoff_fc = 22,
                 cutoff_vel = 10, cut_pitches = [0x2e])

Remove fake hihat notes (typically in a drums recording of the TD50 drumset)
generated by certain actions.
Remove every pitch in `cut_pitches` (by default the Hihat Head event)
with velocity less than `cutoff_fc`
and position maximum `back` ticks before or `forw` ticks after a foot close.
Triggers to search for are pitches are given by the keyword
`triggers = (0x2c, 0x1a, 0x16)` (foot close and hi hat rim).

After that, all notes with pitches in `cut_pitches` and `velocity ≤ cutoff_vel`
are deleted as well, which accounts for spurious note triggering.
"""
function rm_hihatfake!(notes::Notes;
    back = 100, forw = 100, cutoff_fc = 22,
    cutoff_vel = 10, cut_pitches = [0x2e],
    triggers = (0x2c, 0x1a, 0x16))

    #first map special closed notes
    for note in notes
        if note.pitch == 0x16
            note.pitch = 0x1a
        elseif note.pitch == 0x2a
            note.pitch = 0x2e
        end
    end
    #then look for fake notes
   i = 1
   deleted = 0
   len = length(notes)
   while i <= len
      #find foot close
      if notes[i].pitch ∈ triggers
         #go back and remove all fake tip strokes
         j = i-1
         #search all notes in specified back region
         while j>0 && notes[i].position-notes[j].position < back
            #if they are quiet enough
            if notes[j].pitch == 0x2e && notes[j].velocity ≤ cutoff_fc
               #remove them
               deleteat!(notes.notes,j)
               deleted += 1
               i-=1
               len-=1
            else
               j-=1
            end
         end
         #go forward and remove all fake tip strokes
         j=i+1
         #search all notes in specified forw region
         while j ≤ len && notes[j].position-notes[i].position < forw
            #if they are quiet enough
            if notes[j].pitch == 0x2e && notes[j].velocity ≤ cutoff_fc
               #remove them
               deleteat!(notes.notes,j)
               deleted += 1
               len-=1
            else
               j+=1
            end
         end
      end
      i+=1
   end

   deletes = Int[]
   for i = 1:length(notes)
       if notes[i].pitch in cut_pitches && notes[i].velocity ≤ cutoff_vel
           push!(deletes, i)
       end
   end
   deleted += length(deletes)
   deleteat!(notes.notes, deletes)
   println("rm_hihatfake! deleted $(deleted) notes")
end

###############################################################################
#velocity quantization
###############################################################################

"""
    td50_velquant_interval(notes::MIDI.Notes, numintervals::Int)

Divide the velocity range in `numintervals` intervals and quantize the
velocities of each `Note` to the mean value of all notes of the corresponding
instrument in this interval.
"""
function td50_velquant_interval(notes::MIDI.Notes{N}, numintervals::Int) where {N}
    #get notes separated by pitches
    sep = separatepitches(notes)
    newnotes = MIDI.Notes{N}(Vector{N}(), notes.tpq)

    for pitch in keys(sep)
        #short acces to needed notes
        pitchnotes = sep[pitch].notes

        #take care of different maximum velocities
        maxvel = 0
        pitch in DIGITAL ? maxvel = 160 : maxvel = 128

        #do a histogram and weight it with the velocities
        hist = zeros(maxvel)
        for note in pitchnotes
            hist[note.velocity] += 1
        end
        whist = copy(hist)
        for i = 1:length(hist)
            whist[i] *= i
        end

        #create the partitioning and compute corresponding means
        intlength = ceil(Int, maxvel/numintervals)
        meanvals = zeros(Int, numintervals)
        for i = 0:numintervals-1
            start = i*intlength+1
            ende = i*intlength+intlength
            if ende > maxvel
                ende = maxvel
            end
            piece = whist[start:ende]
            hits = sum(hist[start:ende])
            if hits != 0
                piece ./= hits
            end
            meanvals[i+1] = round(Int,mean(piece)*intlength)
        end

        #quantize notes
        for note in pitchnotes
            quant = ceil(Int, note.velocity/intlength)
            note.velocity = meanvals[quant]
        end

        # append to field of quantized notes
        append!(newnotes.notes, pitchnotes)
    end

    #restore temporal order
    sort!(newnotes.notes, lt=((x, y)->x.position<y.position))
    return newnotes
end

end


text:
const RequiemFusion = Dict(
    :name => "RequiemFusion",
    :current_members => [
        "AngelaPapi",
        "Annanana",
        "Anraka",
        "Berserkresse",
        "bnv",
        "Breyn",
        "ChacoKevy",
        "deathmonk",
        "DocStrange",
        "Dom",
        "Drizzt4Ever",
        "F33T",
        "Faeyla",
        "Fisichella",
        "Foxgloves",
        "Imper0s",
        "IrishWrath",
        "KARMA",
        "Kefass",
        "Kobun",
        "Mbad",
        "Mikaeus",
        "Minimoose27",
        "MrMilanna",
        "MrMoJoRisin",
        "Psykosis",
        "SpeedyMoon",
        "uncleb",
        "VladderPfahler",
        "WhoTFCars"
        ],
    :past_members => [
        "Rhellik",
        "Tristanslav",
        "Fryseboksen",
        "sillien",
        "TRUNG",
        "Black Dahlia",
        "matrixttd",
        "Azrael",
        "D_Mac",
        "Glory",
        "Kropman",
        "Lux",
        "Mazsola",
        "skinwalker",
        "The Wiza",
        "Displace",
        ]
    )

const BoethiasChampions = Dict(
    :name => "BoethiasChampions",
    :current_members => [
        "Alnitak",
        "Antar",
        "Aras-LT",
        "Ashen",
        "Chizzy",
        "Clam",
        "Dark",
        "Defeated",
        "Eraser",
        "Flmfx",
        "GodofWar",
        "Grakas",
        "Hellsorcerer",
        "Hollaaa",
        "Hsubox",
        "Lie_to_me",
        "loki",
        "Malachai",
        "Mantis",
        "Mellow",
        "Myname",
        "Peppe",
        "Pigbag",
        "Rukes",
        "Sagitarius",
        "Saint",
        "Sicarium",
        "Skytng",
        "Spinesnap"
    ],
    :past_members => [
        "Ilya"
    ]
)


text:
"""
Return the LOCUS entry of the header.
"""
function parseheader(header::String)
    lines = split(header, "\n")
    locus = split(lines[1], r" +")[2]
    return locus
end

"""
Parse lines encoding genomic position, returning the feature as a `String`, and an instance of `Locus`.
"""
function parseposition(line::String)
    feature, posstring = split(strip(line), r" +")
    if occursin(r"(\.\.|\^)", posstring)
        position = UnitRange(parse.(Int, filter.(c->isnumeric(c), split(posstring, r"(\.\.|\^)(.*(\.\.|\^))?")))...)
        strand = occursin("complement", posstring) ? '-' : '+'
    else
        position = parse(Int, posstring)
        position = position:position
        strand = '.'
    end
    complete_left = !occursin('<', posstring)
    complete_right = !occursin('>', posstring)
    order = Vector{UnitRange{Int}}()
	join = occursin("join", posstring)
    if join || occursin("order", posstring)
        for m in eachmatch(r"\d+(\.\.|\^)\d+", posstring)
            r = Meta.parse.(split(m.match, r"(\.\.|\^)"))
            push!(order, r[1]:r[2])
        end
    end
    return feature, Locus(position, strand, complete_left, complete_right, order, join)
end

"""
Parse footer (sequence) portion of a GenBank file, returning a `String`.
"""
function filterseq(io::IOBuffer)
    line = String(take!(io))
    String(replace(filter(x -> !(isspace(x) || isnumeric(x)), line), r"[^ATGCNatgcn]" => "n"))
end

"""
    parsechromosome(lines)

Parse and return one chromosome entry, and the line number that it ends at.
"""
function parsechromosome(lines)
    genes = Gene[]
    iobuffer = IOBuffer()
	isheader = true
    isfooter = false
    spanning = false
    position_spanning = false
    qualifier = String("")
    content = String("")
    header = ""

    feature = ""
    locus = Locus()

    chromosome = Chromosome()

    linecount = 0
    for line in lines
        linecount += 1

        # Catch cases where there's no header
        if linecount == 1 && occursin(r"gene", line)
            isheader = false
        end

        ### HEADER
		if isheader && occursin(r"FEATURES", line)
            chromosome.header = String(take!(iobuffer))
			isheader = false

        elseif isheader
            linecount == 1 ? print(iobuffer, line) : print(iobuffer, '\n', line)

        # Check if the footer has been reached
		elseif !isheader && !isfooter && (occursin(r"^BASE COUNT", line) || occursin(r"ORIGIN", line))
			# Stop parsing the file when the list of genes is over
            isfooter = true
			iobuffer = IOBuffer()

        ### BODY
		elseif !isheader && !isfooter
            if position_spanning && occursin(r"  /", line)
                position_spanning = false
                spanningline = filter(x -> x != ' ', String(take!(iobuffer)))
                try
                    feature, locus = parseposition(spanningline)
                catch
                    println(spanningline)
                    println(line)
                    @error "parseposition(spanningline) failed at line $linecount"
                end
            elseif position_spanning
                print(iobuffer, line)
            end
            if occursin(r"^ {5}\S", line)
                spanning = false
                try
                    feature, locus = parseposition(line)
                catch
                    println(line)
                    @error "parseposition(line) failed at line $linecount"
                end
                addgene!(chromosome, feature, locus)
            elseif !spanning && occursin(r"^ +/", line)
                if occursin(r"=", line)
                    if occursin("=\"", line)
                        (qualifier, content) = match(r"^ +/([^=]+)=\"?([^\"]*)\"?$", line).captures
                        content = String(content)
                    else
                        (qualifier, content) = match(r"^ +/(\S+)=(\S+)$", line).captures
                        try
                            content = Meta.parse(content)
                            content isa Expr && throw(Meta.ParseError)
                        catch
                            content = Symbol(content)
                        end
                    end

                    if occursin(r"\".*[^\"]$", line)
                        spanning = true
                    end

                    pushproperty!(chromosome.genes[end], Symbol(qualifier), content)

                else
                    # Qualifiers without a value assigned to them end up here
                    qualifier = split(line, '/')[2]
                    pushproperty!(chromosome.genes[end], Symbol(qualifier), true)
                end
            elseif spanning
                try
                    content = match(r" {21}([^\"]*)\"?$", line)[1]
                catch
                    @warn "Couldn't read content (line $linecount)"
                end
                if line[end] == '"'
                    spanning = false
                end
                if eltype(chromosome.genedata[!, Symbol(qualifier)]).b <: AbstractArray
                    i = chromosome.genes[end].index
                    chromosome.genedata[!, Symbol(qualifier)][end][end] = Base.getproperty(chromosome.genes[end], Symbol(qualifier))[end] * "\n" * content
                else
                    Base.setproperty!(chromosome.genes[end], Symbol(qualifier), Base.getproperty(chromosome.genes[end], Symbol(qualifier)) * "\n" * content)
                end
            end

        ### FOOTER
        elseif isfooter
            if line == "//"
                break
            end
            if occursin(r"^ ", line)
                print(iobuffer, line)
            end
        end
    end
    chromosome.name = parseheader(chromosome.header)
    chromosome.sequence = LongDNASeq(filterseq(iobuffer))
    return linecount, chromosome
end

"""
    readgbk(filename)

Parse GenBank-formatted file `filename`, returning a `Vector{Chromosome}`.
"""
function readgbk(filename)
    gz = filename[end-2:end] == ".gz"
    finished = false
    chrs = Chromosome[]
    if gz
        f = GZip.open(filename)
    else
        f = open(filename)
    end
    lines = readlines(f)
    currentline = 1
    while !finished
        if currentline >= length(lines)
            break
        end
        i, chr = parsechromosome(lines[currentline:end])
        currentline += i
        push!(chrs, chr)
    end
    return chrs
end


text:
for (func, typ, complex_typ) in [(:clblasCher2k, cl.CL_float, CL_float2),
                                 (:clblasZher2k, cl.CL_double, CL_double2)]

    @eval @blasfun $func(order::clblasOrder, uplo::clblasUplo,
                              trans::clblasTranspose,
                              N::Csize_t, K::Csize_t,
                              alpha::$complex_typ,
                              A::cl.CL_mem, offa::Csize_t, lda::Csize_t,
                              B::cl.CL_mem, offb::Csize_t, ldb::Csize_t,
                              beta::$typ,
                              C::cl.CL_mem, offc::Csize_t, ldc::Csize_t,
                              n_queues::cl.CL_uint,
                              queues::Ptr{cl.CL_command_queue},
                              n_events_in_wait_list::cl.CL_uint,
                              event_wait_list::Ptr{cl.CL_event},
                              events::Ptr{cl.CL_event})

    @eval @blasfun2 $func(order::clblasOrder, uplo::clblasUplo,
                               trans::clblasTranspose,
                               N::Csize_t, K::Csize_t,
                               alpha::$complex_typ,
                               A::cl.CL_mem, offa::Csize_t, lda::Csize_t,
                               B::cl.CL_mem, offb::Csize_t, ldb::Csize_t,
                               beta::$typ,
                               C::cl.CL_mem, offc::Csize_t, ldc::Csize_t)

end


text:
import SpecialFunctions: gamma, lgamma, erf, erfc, beta

import Base.FastMath: sqrt_llvm_fast, neg_float_fast,
    add_float_fast, sub_float_fast, mul_float_fast, div_float_fast, rem_float_fast,
    eq_float_fast, ne_float_fast, lt_float_fast, le_float_fast,
    sign_fast, abs_fast, log_fast, log2_fast, log10_fast, sqrt_fast,
    pow_fast

# Implicitly defined:
#binary: mod
#unary: exp, abs, log1p, sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh, atanh, erf, erfc, gamma, relu, round, floor, ceil, round, sign.

# Use some fast operators from https://github.com/JuliaLang/julia/blob/81597635c4ad1e8c2e1c5753fda4ec0e7397543f/base/fastmath.jl
# Define allowed operators. Any julia operator can also be used.
plus(x::Float64, y::Float64)::Float64 = add_float_fast(x, y) #Do not change the name of this operator.
sub(x::Float64, y::Float64)::Float64 = sub_float_fast(x, y) #Do not change the name of this operator.
mult(x::Float64, y::Float64)::Float64 = mul_float_fast(x, y) #Do not change the name of this operator.
square(x::Float64)::Float64 = mul_float_fast(x, x)
cube(x::Float64)::Float64 = mul_float_fast(mul_float_fast(x, x), x)
pow(x::Float64, y::Float64)::Float64 = sign_fast(x)*pow_fast(abs(x), y)
div(x::Float64, y::Float64)::Float64 = div_float_fast(x, y)
logm(x::Float64)::Float64 = log_fast(abs_fast(x) + 1e-8)
logm2(x::Float64)::Float64 = log2_fast(abs_fast(x) + 1e-8)
logm10(x::Float64)::Float64 = log10_fast(abs_fast(x) + 1e-8)
sqrtm(x::Float64)::Float64 = sqrt_fast(abs_fast(x))
neg(x::Float64)::Float64 = neg_float_fast(x)

function greater(x::Float64, y::Float64)::Float64
    if x > y
        return 1e0
    end
    return 0e0
end

function relu(x::Float64)::Float64
    if x > 0e0
        return x
    end
    return 0e0
end

function logical_or(x::Float64, y::Float64)::Float64
    if x > 0e0 || y > 0e0
        return 1e0
    end
    return 0e0
end

# (Just use multiplication normally)
function logical_and(x::Float64, y::Float64)::Float64
    if x > 0e0 && y > 0e0
        return 1e0
    end
    return 0e0
end

max_stars_repo_name: LuchkoLab/PySR

max_issues_repo_name: LuchkoLab/PySR

max_forks_repo_name: LuchkoLab/PySR


text:
@testset "evaluate" begin
    ground_truth = vcat(fill(:outlier,30), fill(:inlier, 20))
    classification = vcat(fill(:outlier,20), fill(:inlier, 10), fill(:outlier,5), fill(:inlier,15))

    positive_class = :outlier
    negative_class = :inlier
    @testset "confusion matrix" begin
        conf_matrix = OneClassActiveLearning.ConfusionMatrix(classification, ground_truth)
        @test conf_matrix.tp == 20
        @test conf_matrix.fp == 5
        @test conf_matrix.fn == 10
        @test conf_matrix.tn == 15
        @test OneClassActiveLearning.get_positive(conf_matrix) == sum(ground_truth .== positive_class)
        @test OneClassActiveLearning.get_negative(conf_matrix) == sum(ground_truth .== negative_class)
        @test OneClassActiveLearning.get_n(conf_matrix) == length(ground_truth)
    end

    @testset "matthews corr" begin
        conf_matrix = OneClassActiveLearning.ConfusionMatrix(2, 3, 4, 5)
        numerator = (2 * 4) - (3 * 5)
        denominator = sqrt((2 + 3) * (2 + 5) * (4 + 3) * (4 + 5))
        expected = numerator/denominator
        @test OneClassActiveLearning.matthews_corr(conf_matrix) ≈ expected
        @test OneClassActiveLearning.matthews_corr(OneClassActiveLearning.ConfusionMatrix(0, 0, 4, 5)) ≈ 0
    end

    @testset "roc auc" begin
        predictions = vcat(rand(10), -rand(20))
        ground_truth = vcat(fill(:outlier, 10), fill(:inlier, 20))

        @test OneClassActiveLearning.roc_auc(predictions, ground_truth) ≈ 1.0
        @test OneClassActiveLearning.roc_auc(predictions, ground_truth, fpr = 0.1, normalize=true) ≈ 1.0
        ground_truth = vcat(fill(:inlier, 10), fill(:outlier, 20))
        @test OneClassActiveLearning.roc_auc(predictions, ground_truth) ≈ 0.0
        ground_truth = vcat(fill(:inlier, 5), fill(:outlier, 25))
        @test OneClassActiveLearning.roc_auc(predictions, ground_truth, fpr = 0.04, normalize=true) .> OneClassActiveLearning.roc_auc(predictions, ground_truth)
    end
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra
using SpecialFunctions
using SparseArrays

@testset "compute_types" begin
    A0=(1:5)*(1:5)'
    A1=(1:5)*(3:7)'+I;
    nep=DEP([A0,A1]);
    x=ones(5);
    s=compute_rf(ComplexF64,nep,x,
                 NonlinearEigenproblems.NEPSolver.ScalarNewtonInnerSolver(),λ=-0.5+3im)[1];

    @test abs(x'*compute_Mlincomb(nep,s,x)) < eps()*100;

    y=zeros(5); y[1]=3;
    s=compute_rf(ComplexF64,nep,x,
                 NonlinearEigenproblems.NEPSolver.ScalarNewtonInnerSolver(),λ=-0.5+3im,
                 y=y)[1];

    @test abs(y'*compute_Mlincomb(nep,s,x)) < eps()*100;
end


text:
module effective_interaction

# Optimized functions for calculation of Omega and Gamma

function Omega(a, θ)
    ξ = 2*pi*a
    cosθpow2 = cos(θ)^2
    cosξdivξ = cos(ξ)/ξ
    return 3.0/4*((1.0-3.0*cosθpow2) * (sin(ξ)+cosξdivξ)/ξ^2 - (1.0-cosθpow2)*cosξdivξ)
end

function Gamma(a, θ)
    ξ = 2*pi*a
    cosθpow2 = cos(θ)^2
    if ξ<1e-3
        return 3/2*(2/3 + 1/15*ξ^2 *(-2 + cosθpow2) + ξ^4*(1/140-cosθpow2/210))
    else
        sinξ = sin(ξ)/ξ
        return 3/2*((1-cosθpow2)*sinξ + (1-3*cosθpow2) * (cos(ξ)-sinξ)/ξ^2)
    end
end

function OmegaGamma(a, cosθpow2)
    #cosθpow2 = cos(θ)^2
    cosξ = cospi(2*a)
    sinξ = sinpi(2*a)
    omega = 0.75*((1.0-3.0*cosθpow2) * (sinξ+cosξ/(2*pi*a))/(2*pi*a)^2 - (1.0-cosθpow2)*cosξ/(2*pi*a))
    gamma = 1.5*((1-cosθpow2)*sinξ/(2*pi*a) + (1-3*cosθpow2) * (cosξ-sinξ/(2*pi*a))/(2*pi*a)^2)
    return omega, gamma
end

function OmegaGamma_orthogonal(a)
    cosξ = cospi(2*a)
    sinξ = sinpi(2*a)
    omega = 0.75*((sinξ+cosξ/(2*pi*a))/(2*pi*a)^2 - cosξ/(2*pi*a))
    gamma = 1.5*(sinξ/(2*pi*a) + (cosξ-sinξ/(2*pi*a))/(2*pi*a)^2)
    return omega, gamma
end

function Omega_orthogonal(a)
    ξ = 2*pi*a
    cosξdξ = cos(ξ)/ξ
    return 3.0/4.0*(-cosξdξ + (sin(ξ)+cosξdξ)/ξ^2)
end

function Gamma_orthogonal(a)
    ξ = 2*a
    sincξ = sinc(ξ)
    return 3.0/2.0*(sincξ + (cospi(ξ)-sincξ)/(pi*ξ)^2)
end

function map_effectiveinteraction(a::Vector{Float64}, f::Function)
    omega_list = Float64[]
    gamma_list = Float64[]
    for ai=a
        omega, gamma = f(ai, e_dipole)
        push!(omega_list, omega)
        push!(gamma_list, gamma)
    end
    return omega_list, gamma_list
end

# Finite symmetric systems

"""
    effective_interaction.triangle_orthogonal(a)

Effective Omega and Gamma for a equilateral triangle with edge length `a`.

The polarization axis is orthogonal to the triangle plane.
"""
function triangle_orthogonal(a::Float64)
    omega_eff = 2*Omega_orthogonal(a)
    gamma_eff = 2*Gamma_orthogonal(a)
    return omega_eff, gamma_eff
end

"""
    effective_interaction.square_orthogonal(a)

Effective Omega and Gamma for a square with edge length `a`.

The polarization axis is orthogonal to the square plane.
"""
function square_orthogonal(a::Float64)
    omega_eff = 2*Omega_orthogonal(a) + Omega_orthogonal(sqrt(2.)*a)
    gamma_eff = 2*Gamma_orthogonal(a) + Gamma_orthogonal(sqrt(2.)*a)
    return omega_eff, gamma_eff
end

"""
    effective_interaction.rectangle_orthogonal(a, b)

Effective Omega and Gamma for a rectangle with edge lengths `a` and `b`.

The polarization axis is orthogonal to the rectangle plane.
"""
function rectangle_orthogonal(a::Float64, b::Float64)
    d = sqrt(a^2 + b^2)
    omega_eff, gamma_eff = OmegaGamma_orthogonal(a)
    om, ga = OmegaGamma_orthogonal(b)
    omega_eff += om; gamma_eff+=ga
    om, ga = OmegaGamma_orthogonal(d)
    omega_eff += om; gamma_eff+=ga
    return omega_eff, gamma_eff
end

"""
    effective_interaction.cube_orthogonal(a)

Effective Omega and Gamma for a cube with edge length `a`

The polarization axis is orthogonal to the xy faces.
"""
function cube_orthogonal(a::Float64)
    omega_eff, gamma_eff = square_orthogonal(a)
    sqrt2 = sqrt(2.)
    sqrt3 = sqrt(3.)
    Θdiag = atan(sqrt2, 1.)
    omega_eff += Omega(a, 0.) + 2*Omega(sqrt(2.)*a, pi/4.) + Omega(sqrt3*a, Θdiag)
    gamma_eff += Gamma(a, 0.) + 2*Gamma(sqrt(2.)*a, pi/4.) + Gamma(sqrt3*a, Θdiag)
    return omega_eff, gamma_eff
end

"""
    effective_interaction.box_orthogonal(a, b, c)

Effective Omega and Gamma for a box with edge lengths `a`, `b` and `c`.

The polarization axis is orthogonal to the top face.
"""
function box_orthogonal(a::Float64, b::Float64, c::Float64)
    omega_eff, gamma_eff = rectangle_orthogonal(a, b)

    a2 = a^2
    b2 = b^2
    c2 = c^2

    om, ga = OmegaGamma(c, 1.)
    omega_eff += om; gamma_eff+=ga

    dxz_pow2 = a2+c2
    om, ga = OmegaGamma(sqrt(dxz_pow2), c2/dxz_pow2)
    omega_eff += om; gamma_eff+=ga

    dyz_pow2 = b2+c2
    om, ga = OmegaGamma(sqrt(dyz_pow2), c2/dyz_pow2)
    omega_eff += om; gamma_eff+=ga

    dxyz_pow2 = a2+b2+c2
    om, ga = OmegaGamma(sqrt(dxyz_pow2), c2/dxyz_pow2)
    omega_eff += om; gamma_eff+=ga
    return omega_eff, gamma_eff
end

# Infinite 1D symmetric systems
"""
    effective_interaction.chain(a, Θ, N)

Effective Omega and Gamma for an infinite chain.

The calculation is done by adding N spins left and N spins right of a
central spin.

# Arguments
* `a`: Spin-spin distance.
* `θ`: Angle between polarization axis and spin chain.
* `N`: Number of included spins.
"""
function chain(a::Float64, Θ, N::Int)
    omega_eff::Float64 = 0.
    gamma_eff::Float64 = 0.
    d = 0.
    for j=1:N
        d += a
        omega_eff += Omega(d, Θ)
        gamma_eff += Gamma(d, Θ)
    end
    return 2*omega_eff, 2*gamma_eff
end

"""
    effective_interaction.chain_orthogonal(a, N)

Effective Omega and Gamma for an infinite chain with orthogonal polarization axis.

The calculation is done by adding N spins left and N spins right of a
central spin.

# Arguments
* `a`: Spin-spin distance.
* `N`: Number of included spins.
"""
function chain_orthogonal(a::Float64, N::Int)
    omega_eff::Float64 = 0.
    gamma_eff::Float64 = 0.
    d::Float64 = 0.
    for j=1:N
        d += a
        omega_eff += Omega_orthogonal(d)
        gamma_eff += Gamma_orthogonal(d)
    end
    return 2*omega_eff, 2*gamma_eff
end

# Infinite 2D symmetric systems
"""
    effective_interaction.squarelattice_orthogonal(a, N)

Effective Omega and Gamma for a infinite square lattice.

The polarization axis is orthogonal to the square lattice plane and the
calculation is done by creating a (2N+1)*(2N+1) square lattice
and calculate the combined interaction for the central spin.

# Arguments
* `a`: Spin-spin distance.
* `N`: Number of included spins.
"""
function squarelattice_orthogonal(a::Float64, N::Int)
    omega_eff::Float64 = 0.
    gamma_eff::Float64 = 0.
    for j=1:N
        d = a*j
        omega_eff += 4*Omega_orthogonal(d)
        gamma_eff += 4*Gamma_orthogonal(d)
        dx2 = j^2
        for i=1:j-1
            dy2 = i^2
            d = a*sqrt(dx2+dy2)
            omega_eff += 8*Omega_orthogonal(d)
            gamma_eff += 8*Gamma_orthogonal(d)
        end
        d = a*sqrt(2)*j
        omega_eff += 4*Omega_orthogonal(d)
        gamma_eff += 4*Gamma_orthogonal(d)
    end
    return omega_eff, gamma_eff
end

"""
    effective_interaction.hexagonallattice_orthogonal(a, N)

Effective Omega and Gamma for a infinite hexagonal lattice.

The polarization axis is orthogonal to the square lattice plane and the
calculation is done by creating hexagonal lattice consisting of
N rings and calculate the combined interaction for the central spin.

# Arguments
* `a`: Spin-spin distance.
* `N`: Number of included spins.
"""
function hexagonallattice_orthogonal(a::Float64, N::Int)
    omega_eff::Float64 = 0.
    gamma_eff::Float64 = 0.
    for n=1:2:N
        nv1x = -0.5*n
        nv1y_square = 3.0/4*n^2
        d = a*sqrt(nv1y_square + nv1x^2)
        omega_eff += 6*Omega_orthogonal(d)
        gamma_eff += 6*Gamma_orthogonal(d)
        for i=1:div(n-1,2)
            d = a*sqrt(nv1y_square + (nv1x+1.0*i)^2)
            omega_eff += 12*Omega_orthogonal(d)
            gamma_eff += 12*Gamma_orthogonal(d)
        end
    end
    for n=2:2:N
        nv1x = -0.5*n
        nv1y_square = 3.0/4*n^2
        d = a*sqrt(nv1y_square + nv1x^2)
        omega_eff += 6*Omega_orthogonal(d)
        gamma_eff += 6*Gamma_orthogonal(d)
        for i=1:div(n-2,2)
            d = a*sqrt(nv1y_square + (nv1x+1.0*i)^2)
            omega_eff += 12*Omega_orthogonal(d)
            gamma_eff += 12*Gamma_orthogonal(d)
        end
        i = div(n-2,2)+1
        d = a*sqrt(nv1y_square + (nv1x+1.0*i)^2)
        omega_eff += 6*Omega_orthogonal(d)
        gamma_eff += 6*Gamma_orthogonal(d)
    end
    return omega_eff, gamma_eff
end

# Infinite 3D symmetric systems
"""
    effective_interaction.cubiclattice_orthogonal(a, N)

Effective Omega and Gamma for a infinite cubic lattice.

The polarization axis is orthogonal to the top face of a unit cell and the
calculation is done by creating a (2N+1)*(2N+1)*(2N+1) cubic lattice
and calculate the combined interaction for the central spin.

# Arguments
* `a`: Spin-spin distance.
* `N`: Number of included spins.
"""
function cubiclattice_orthogonal(a::Float64, N::Int)
    omega_eff, gamma_eff = squarelattice_orthogonal(a, N)
    for nx=1:N
        for ny=0:nx
            multiplicity = ((ny==0 || ny==nx) ? 8 : 16)
            nrpow2 = nx^2 + ny^2
            nr = sqrt(nrpow2)
            for nz=1:N
                nzpow2 = nz^2
                d = a*sqrt(nrpow2+nzpow2)
                cosθpow2 = nzpow2/(nrpow2+nzpow2)
                om, ga = OmegaGamma(d, cosθpow2)
                omega_eff += multiplicity*om
                gamma_eff += multiplicity*ga
            end
        end
    end
    d::Float64 = 0.
    for nz=1:N
        d += a
        om, ga = OmegaGamma(d, 1.)
        omega_eff += 2*om
        gamma_eff += 2*ga
    end
    return omega_eff, gamma_eff
end

"""
    effective_interaction.tetragonallattice_orthogonal(a, b, N)

Effective Omega and Gamma for a infinite tetragonal lattice.

The polarization axis is orthogonal to the top face of a unit cell and the
calculation is done by creating a (2N+1)*(2N+1)*(2N+1) tetragonal lattice
and calculate the combined interaction for the central spin.

# Arguments
* `a`: Spin-spin distance for bottom side square.
* `b`: Height of the unit cell.
* `N`: Number of included spins.
"""
function tetragonallattice_orthogonal(a::Float64, b::Float64, N::Int)
    omega_eff, gamma_eff = squarelattice_orthogonal(a, N)
    for nx=1:N
        for ny=0:nx
            multiplicity = ((ny==0 || ny==nx) ? 8 : 16)
            nrpow2 = nx^2 + ny^2
            nr = sqrt(nrpow2)
            for nz=1:N
                nzpow2 = nz^2
                d = sqrt(a^2*nrpow2 + b^2*nzpow2)
                θ = atan(a*nr, b*nz)
                omega_eff += multiplicity*Omega(d, θ)
                gamma_eff += multiplicity*Gamma(d, θ)
            end
        end
    end
    d::Float64 = 0.
    for nz=1:N
        d += b
        omega_eff += 2*Omega(d, 0.)
        gamma_eff += 2*Gamma(d, 0.)
    end
    return omega_eff, gamma_eff
end

"""
    effective_interaction.hexagonallattice3d_orthogonal(a, b, N)

Effective Omega and Gamma for a infinite 3D hexagonal lattice.

The lattice consists of stacked planes of hexagonal lattices where the
the polarization axis is orthogonal to the planes. The
calculation is done by creating hexagonal lattices with N rings, stacking
2N+1 lattices of this kind above each other and calculating the combined
interaction for the central spin.

# Arguments
* `a`: Spin-spin distance for hexagons.
* `b`: Distance between planes of hexagonal lattices
* `N`: Number of included spins.
"""
function hexagonallattice3d_orthogonal(a::Float64, b::Float64, N::Int)
    omega_eff, gamma_eff::Float64 = hexagonallattice_orthogonal(a, N)
    a2 = a^2
    b2 = b^2
    for n=1:2:N
        nv1x = -0.5*n
        nv1y_square = 3.0/4*n^2
        nr_square = a2*(nv1y_square + nv1x^2)
        for iz=1:N
            nz_square = iz^2*b2
            d = sqrt(nr_square+nz_square)
            cosθ_square = nz_square/(nr_square+nz_square)
            om, ga = OmegaGamma(d, cosθ_square)
            omega_eff += 12*om
            gamma_eff += 12*ga
        end
        for i=1:div(n-1,2)
            nr_square = a2*(nv1y_square + (nv1x+1.0*i)^2)
            for iz=1:N
                nz_square = iz^2*b2
                d = sqrt(nr_square+nz_square)
                cosθ_square = nz_square/(nr_square+nz_square)
                om, ga = OmegaGamma(d, cosθ_square)
                omega_eff += 24*om
                gamma_eff += 24*ga
            end
        end
    end
    for n=2:2:N
        nv1x = -0.5*n
        nv1y_square = 3.0/4*n^2
        nr_square = a2*(nv1y_square + nv1x^2)
        for iz=1:N
            nz_square = iz^2*b2
            d = sqrt(nr_square+nz_square)
            cosθ_square = nz_square/(nr_square+nz_square)
            om, ga = OmegaGamma(d, cosθ_square)
            omega_eff += 12*om
            gamma_eff += 12*ga
        end
        for i=1:div(n-2, 2)
            nr_square = a2*(nv1y_square + (nv1x+1.0*i)^2)
            for iz=1:N
                nz_square = iz^2*b2
                d = sqrt(nr_square+nz_square)
                cosθ_square = nz_square/(nr_square+nz_square)
                om, ga = OmegaGamma(d, cosθ_square)
                omega_eff += 24*om
                gamma_eff += 24*ga
            end
        end
        i = div(n-2, 2)+1
        nr_square = a2*(nv1y_square + (nv1x+1.0*i)^2)
        for iz=1:N
            nz_square = iz^2*b2
            d = sqrt(nr_square+nz_square)
            cosθ_square = nz_square/(nr_square+nz_square)
            om, ga = OmegaGamma(d, cosθ_square)
            omega_eff += 12*om
            gamma_eff += 12*ga
        end
    end
    return omega_eff, gamma_eff
end

end # module