text:
#===============================================================================
# Variable scope in julia
#
# Author: Carl Cortright
# Date: 12/28/2016
===============================================================================#

module bar
x = 1
f(x) = x
end

import bar
# Now we can use the variable x
println(bar.x)
for i in 1:10
  local i # This forces i in to local scope
  println(i)
end
println(i)

# In a soft local scope, all variables are inherited from the parent scope
y = 2
for i in 1:7
  println(y)
end

# In hard local scopes, like functions and macros, an explicit global is needed
function hardLocal()
  global y = 7
end

hardLocal()
println(y)

# Nested functions can modify their parent functions local variables
function hello()
  x = 10
  function world()
    x = 20
  end
  world()
  println(x)
end
hello()

# Let blocks allocate new memory every time they run
let y = 10
  println(y)
end
println(y)

# we can also define constants
const somenum = 0.1234567891

max_issues_repo_name: carlcortright/Playground


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_status(path::AbstractString, inputs::Dict, EP::Model)

Function for writing the final solve status of the optimization problem solved.
"""
function write_status(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

	# https://jump.dev/MathOptInterface.jl/v0.9.10/apireference/#MathOptInterface.TerminationStatusCode
	status = termination_status(EP)

	# Note: Gurobi excludes constants from solver reported objective function value - MIPGap calculated may be erroneous
	if (setup["UCommit"] == 0 || setup["UCommit"] == 2)
		dfStatus = DataFrame(Status = status, Solve = inputs["solve_time"],
			Objval = objective_value(EP))
	else
		dfStatus = DataFrame(Status = status, Solve = inputs["solve_time"],
			Objval = objective_value(EP), Objbound= objective_bound(EP),FinalMIPGap =(objective_value(EP) -objective_bound(EP))/objective_value(EP) )
	end
	CSV.write(joinpath(path, "status.csv"),dfStatus)
end

ext: jl

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
const lib = joinpath(pwd(), "async.so")

done = true

function job(async)
  global done
  sleep(2)  # heavy task
  println("job is done")
  done = true
end

cond = Base.AsyncCondition(job)

function CB(handle)
  global done
  done = false
  ccall(:uv_async_send, Cint, (Ptr{Cvoid},), handle)
  while !done end
  return
end

cb_ptr = @cfunction(CB, Cvoid, (Ptr{Cvoid},))

ccall((:set,  lib), Cvoid, (Ptr{Cvoid},), cb_ptr)
ccall((:call, lib), Cvoid, (Ptr{Cvoid},), cond.handle)

max_stars_repo_name: Astellon/julia-callback-from-c

max_issues_repo_name: Astellon/julia-callback-from-c

max_forks_repo_name: Astellon/julia-callback-from-c


text:
HEIGHT = 650
WIDTH = 400
SPEED = 2
GAP = 100

bird = Actor("bird1")
bird.pos = (75, 200)
pipe_top = Actor("top")
pipe_bottom = Actor("bottom")

bird_dead=false
bird_vy=0.0
GRAVITY=0.3
FLAP_STRENGTH=5.0

function reset_pipes()
    pipe_gap_y = rand(200:(HEIGHT-200))
    pipe_top.bottomleft = (WIDTH, pipe_gap_y - GAP / 2)
    pipe_bottom.topleft = (WIDTH, pipe_gap_y + GAP / 2)
end

reset_pipes()

function draw()
    draw(bird)
    draw(pipe_top)
    draw(pipe_bottom)
end

function update()
    update_bird()
    update_pipes()
end

function update_bird()
    global bird_vy
    global bird_dead
    uy=bird_vy
    bird_vy=bird_vy+GRAVITY
    bird.x=75
    bird.y+=Int(round((uy+bird_vy)/2))
    if collide(bird,pipe_top) || collide(bird,pipe_bottom)
        bird_dead=true
        bird.image="birddead"
    end
end

function update_pipes()
    pipe_top.left = pipe_top.left - SPEED
    pipe_bottom.left = pipe_bottom.left - SPEED
    if pipe_top.right<0
        reset_pipes()
    end
end

function on_key_down(g::Game)
    global bird_vy
    if !bird_dead
        bird_vy = -FLAP_STRENGTH
    end
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


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

mutable struct IoK8sApiFlowcontrolV1alpha1QueuingConfiguration <: SwaggerModel
    handSize::Any # spec type: Union{ Nothing, Int32 } # spec name: handSize
    queueLengthLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: queueLengthLimit
    queues::Any # spec type: Union{ Nothing, Int32 } # spec name: queues

    function IoK8sApiFlowcontrolV1alpha1QueuingConfiguration(;handSize=nothing, queueLengthLimit=nothing, queues=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1QueuingConfiguration, Symbol("handSize"), handSize)
        setfield!(o, Symbol("handSize"), handSize)
        validate_property(IoK8sApiFlowcontrolV1alpha1QueuingConfiguration, Symbol("queueLengthLimit"), queueLengthLimit)
        setfield!(o, Symbol("queueLengthLimit"), queueLengthLimit)
        validate_property(IoK8sApiFlowcontrolV1alpha1QueuingConfiguration, Symbol("queues"), queues)
        setfield!(o, Symbol("queues"), queues)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1QueuingConfiguration

const _property_map_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration = Dict{Symbol,Symbol}(Symbol("handSize")=>Symbol("handSize"), Symbol("queueLengthLimit")=>Symbol("queueLengthLimit"), Symbol("queues")=>Symbol("queues"))
const _property_types_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration = Dict{Symbol,String}(Symbol("handSize")=>"Int32", Symbol("queueLengthLimit")=>"Int32", Symbol("queues")=>"Int32")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1QueuingConfiguration }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1QueuingConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1QueuingConfiguration }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1QueuingConfiguration[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1QueuingConfiguration)
    true
end

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


text:
## -------------------------------------------------------------------
# EP biomass corr
let
    objider = iJR.BIOMASS_IDER
    ps = Plots.Plot[]
    for method in ALL_METHODS
        p = plot(title = string(iJR.PROJ_IDER, " method: ", method),
            xlabel = "exp biom", ylabel = "model biom")
        ep_vals = DAT[method, :ep, :flx, objider, EXPS]
        eperr_vals = DAT[method, :eperr, :flx, objider, EXPS]
        Fd_vals = DAT[method, :Fd, :flx, objider, EXPS]
        Fderr_vals = DAT[method, :Fderr, :flx, objider, EXPS]
        color = [exp_colors[exp] for exp in EXPS]
        m, M = myminmax([Fd_vals; ep_vals])
        margin = abs(M - m) * 0.1
        scatter!(p, Fd_vals, ep_vals;
            xerr = Fderr_vals,
            yerr = eperr_vals,
            label = "", color,
            alpha = 0.7, ms = 7,
            xlim = [m - margin, M + margin],
            ylim = [m - margin, M + margin],
        )
        push!(ps, p)
    end
    layout = (1, length(ps))
    mysavefig(ps, "obj_val_ep_corr"; layout)
end

## -------------------------------------------------------------------
# correlations
let

    tot_ps = Plots.Plot[]
    for method in ALL_METHODS
        # total corr
        let
            ep_vals = DAT[method, :ep, :flx, EXCH_FLX_IDERS, EXPS] .|> abs
            ep_errs = DAT[method, :eperr, :flx, EXCH_FLX_IDERS, EXPS] .|> abs
            Fd_vals = DAT[method, :Fd, :flx, EXCH_FLX_IDERS, EXPS] .|> abs

            color = [ider_colors[ider] for ider in EXCH_FLX_IDERS, exp in EXPS]
            scatter_params = (;label = "", color, ms = 7, alpha = 0.7)
            # ep corr
            p = plot(title = "$(iJR.PROJ_IDER) (EP) $method",
                ylabel = "model abs flx",
                xlabel = "exp abs flx",
            )
            scatter!(p, Fd_vals, ep_vals; yerr = ep_errs, scatter_params...)
            all_vals = [ep_vals; Fd_vals] |> sort!
            plot!(p, all_vals, all_vals; ls = :dash, color = :black, label = "")
            push!(tot_ps, deepcopy(p))
        end

        # per ider
        let
            for ider in EXCH_FLX_IDERS
                ep_vals = DAT[method, :ep, :flx, ider, EXPS] .|> abs
                ep_errs = DAT[method, :eperr, :flx, ider, EXPS] .|> abs
                Fd_vals = DAT[method, :Fd, :flx, ider, EXPS] .|> abs

                color = ider_colors[ider]
                # ep corr
                p = plot(title = "$(iJR.PROJ_IDER) (EP) $method",
                    ylabel = "model abs flx",
                    xlabel = "exp abs flx",
                )
                scatter!(p, Fd_vals, ep_vals; yerr = ep_errs,
                    label = "", color, ms = 7, alpha = 0.7
                )
                bounds = DAT[method, :bounds, :flx, ider, EXPS]
                all_vals = [ep_vals; Fd_vals] |> sort!
                plot!(p, all_vals, all_vals;
                    ls = :dash, color = :black, label = "",
                    # xlim = [lb, ub], ylim = [lb, ub]
                )
                mysavefig(p, "corr"; ider, method)
            end
        end
    end

    layout = (1, length(tot_ps))
    mysavefig(tot_ps, "flx_tot_corr"; layout)

end


text:
### Skipmer Iterators
###
###
### Iterator over all Skipmers in a biological sequence.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

struct SkipmerFactoryResult{T<:AbstractMer}
    position::Int
    fw::T
    bw::T
end

mutable struct SkipmerFactory{S<:LongNucleotideSequence,U<:Unsigned,K}
    seq::S
    cycle_pos::Vector{UInt8}
    last_unknown::Vector{Int64}
    fkmer::Vector{U}
    rkmer::Vector{U}
    position::Int
    finished::Int
    n::Int
    cycle_len::Int
    bases_per_cycle::Int
    span::Int

    function SkipmerFactory(::Type{T},
                            seq::LongSequence{A},
                            bases_per_cycle::Int = 2,
                            cycle_len::Int = 3) where {A<:NucleicAcidAlphabet,T<:AbstractMer}

        checkmer(T)

        span = UInt(ceil(cycle_len * (ksize(T) / bases_per_cycle - 1) + bases_per_cycle))

        if eltype(seq) ∉ (DNA, RNA)
            throw(ArgumentError("element type must be either DNA or RNA nucleotide"))
        elseif bases_per_cycle > cycle_len
            throw(ArgumentError("bases per cycle must not be greater than the cycle length"))
        elseif span > length(seq)
            throw(ArgumentError(string("span of each skipmer (", span, ") is greater than the input sequence length (", length(seq), ").")))
        elseif eltype(seq) != eltype(T)
            throw(ArgumentError(string("skipmer type chosen must have same element type as chosen sequence")))
        end

        # For each of the next N skipmers being build simultaneously by the iterator,
        # keep track of a cycle_pos variable: it keeps track of which nucleotides
        # in the sequence are appended to said skipmer, and is used by the
        # _consider_position! method.
        cycle_pos = Vector{UInt8}(undef, cycle_len)

        # Vector to keep track of position of last ambiguous nucleotide we saw,
        # for each of the N skipmers being built simultaneously by the iterator.
        last_unknown = Vector{Int64}(undef, cycle_len)

        # Storage for the next N skipers being built simultaneously by the iterator.
        U = encoded_data_type(T)
        fkmer = Vector{U}(undef, cycle_len)
        rkmer = Vector{U}(undef, cycle_len)

        gen = new{LongSequence{A},U,ksize(T)}(seq, cycle_pos, last_unknown, fkmer, rkmer, 0, 0, 0, cycle_len, bases_per_cycle, span)
        _init_generator!(gen)

        return gen
    end
end

@inline function mertype(::Type{SkipmerFactory{S,UInt64,K}}) where {S<:LongNucleotideSequence,K}
    return Mer{minimal_alphabet(typeof(Alphabet(S))),K}
end

@inline function mertype(::Type{SkipmerFactory{S,UInt128,K}}) where {S<:LongNucleotideSequence,K}
    return BigMer{minimal_alphabet(typeof(Alphabet(S))),K}
end

@inline mertype(gen::SkipmerFactory) = mertype(typeof(gen))

@inline function Base.eltype(::Type{S}) where {S<:SkipmerFactory}
    return MerIterResult{mertype(S)}
end

@inline Base.IteratorSize(::Type{T}) where T <: SkipmerFactory = Base.HasLength()

@inline Base.IteratorEltype(::Type{T}) where T <: SkipmerFactory = Base.HasEltype()

@inline Base.length(gen::SkipmerFactory)::Int = length(gen.seq) - gen.span + 1

function set_sequence!(gen::SkipmerFactory{S}, seq::S) where {S<:LongNucleotideSequence}
    gen.seq = seq
    _init_generator!(gen)
end

function _init_generator!(gen::SkipmerFactory)
    N = gen.cycle_len

    @inbounds for i in 1:N
        gen.cycle_pos[i] = N - i
        gen.last_unknown[i] = -1
        gen.fkmer[i] = 0
        gen.rkmer[i] = 0
    end

    gen.finished = 0
    gen.position = 0

    for i in 1:gen.span-1
        _consider_next_position!(gen)
    end
end

@inline function _get_base_bits(seq::LongSequence{<:NucleicAcidAlphabet{2}}, position)
    return extract_encoded_element(bitindex(seq, position), encoded_data(seq))
end

@inline function _get_base_bits(seq::LongSequence{<:NucleicAcidAlphabet{4}}, position)
    return twobitnucs[extract_encoded_element(bitindex(seq, position), encoded_data(seq)) + 0x01]
end

@inline function _append_mers!(gen::SkipmerFactory{LongSequence{A},U,K}, fbits::U, rbits::U) where {A<:NucleicAcidAlphabet{2},U,K}
    # For each of the N next skipmers we are building...
    N = gen.cycle_len
    M = gen.bases_per_cycle
    for ni in 1:N
        # Update the ni'th skipmer's cycle_pos value.
        # This is used to decide if the ni'th skipmer will get the base at the
        # position currently being considered.
        nextcycle = gen.cycle_pos[ni] + 1
        gen.cycle_pos[ni] = ifelse(nextcycle == N, 0, nextcycle)
        # If the ni'th skipmer should get the base at the position currently being
        # considered, then append said base to said skipmer. We do this using some
        # bit-flipping. We also build both the forward and reverse forms of the N
        # skipmers, at the same time for efficiency.
        if gen.cycle_pos[ni] < M
            gen.fkmer[ni] = ((gen.fkmer[ni] << 0x02) | fbits)
            gen.rkmer[ni] = (gen.rkmer[ni] >> 0x02) | (rbits << 2(K - 1))
        end
    end
end

@inline function _append_mers!(gen::SkipmerFactory{LongSequence{A},U,K}, fbits::U, rbits::U) where {A<:NucleicAcidAlphabet{4},U,K}
    N = gen.cycle_len
    M = gen.bases_per_cycle
    if fbits == 0xFF
        for ni in 1:N
            nextcycle = gen.cycle_pos[ni] + 1
            gen.cycle_pos[ni] = ifelse(nextcycle == N, 0, nextcycle)
            if gen.cycle_pos[ni] < M
                gen.fkmer[ni] = (gen.fkmer[ni] << 0x02)
                gen.rkmer[ni] = (gen.rkmer[ni] >> 0x02)
            end
        end
    else
        for ni in 1:N
            nextcycle = gen.cycle_pos[ni] + 1
            gen.cycle_pos[ni] = ifelse(nextcycle == N, 0, nextcycle)
            if gen.cycle_pos[ni] < M
                gen.fkmer[ni] = ((gen.fkmer[ni] << 0x02) | fbits)
                gen.rkmer[ni] = (gen.rkmer[ni] >> 0x02) | (rbits << 2(K - 1))
            end
        end
    end
end

@inline function _consider_next_position!(gen::SkipmerFactory{S,U,K}) where {S,U,K}
    nextposition = gen.position + 1
    fbits = U(_get_base_bits(gen.seq, nextposition))
    gen.position = nextposition
    rbits = ~fbits & U(0x03)
    _append_mers!(gen, fbits, rbits)
end

function nextmer(gen::SkipmerFactory{LongSequence{A},U,K}) where {A<:NucleicAcidAlphabet{2},U,K}
    # If you've hit the end of the sequence, we're done.
    if gen.position == lastindex(gen.seq)
        return nothing
    end
    # _consider_next_position! takes the base at pos, it decides
    # which of the N skipmers being built should get that base, and then it
    # appends the base to those skipmers.
    _consider_next_position!(gen)
    # Now we take the next skipmer that will now be complete. Get its canonical
    # form, and then return it.
    nextfinished = gen.finished + 1
    nextfinished = ifelse(nextfinished > gen.cycle_len, 1, nextfinished)
    gen.finished = nextfinished
    @inbounds fkmer = gen.fkmer[nextfinished]
    @inbounds rkmer = gen.rkmer[nextfinished]
    newn = gen.n + 1
    gen.n = newn
    return MerIterResult(newn, mertype(gen)(fkmer), mertype(gen)(rkmer))
end

function nextmer(gen::SkipmerFactory{LongSequence{A},U,K}) where {A<:NucleicAcidAlphabet{4},U,K}
    lastpos = lastindex(gen.seq)
    span = gen.span
    # For every position in the sequence...
    while gen.position < lastpos
        # Consider the base at the next position, and append it to the relevant
        # skipmers of the N next skipers being built by the iterator.
        _consider_next_position!(gen)

        # If we are at pos, the skip-mer that started at pos-S is now done...
        if gen.position >= span
            nextfinished = gen.finished + 1
            gen.finished = ifelse(nextfinished > gen.cycle_len, 1, nextfinished)
            # If the currently finished skipmer, does not contain an ambiguous
            # nucleotide, then get the canonical form of the skipmer and return it.
            newn = gen.n + 1
            gen.n = newn
            if gen.last_unknown[gen.finished] + span <= gen.position
                fkmer = gen.fkmer[gen.finished]
                rkmer = gen.rkmer[gen.finished]
                return MerIterResult(newn, mertype(gen)(fkmer), mertype(gen)(rkmer))
            end
        end
    end
    # pos => lastpos, so we are done iterating through the sequence.
    return nothing
end

@inline function Base.iterate(gen::SkipmerFactory)
    _init_generator!(gen)
    return iterate(gen, nothing)
end

@inline function Base.iterate(gen::SkipmerFactory, ::Nothing)
    x = nextmer(gen)
    if isnothing(x)
        return nothing
    else
        return x, nothing
    end
end

"""
    each(::Type{T}, seq::BioSequence, bases_per_cycle::Int, cycle_len::Int) where {T<:AbstractMer}

Initialize an iterator over all overlapping skipmers in a sequence `seq`,
skipping ambiguous nucleotides without changing the reading frame.

!!! note
    Please see the BioSequences user manual for more details of how a skipmer is
    constructed.
"""
@inline function each(::Type{T}, seq::BioSequence, bases_per_cycle::Int, cycle_len::Int) where {T<:AbstractMer}
    return SkipmerFactory(T, seq, bases_per_cycle, cycle_len)
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Parallel_Computing")

# # Multithreading
#
# Now we're finally ready to start talking about running things on multiple
# processors! Most computers (even cell phones) these days have multiple cores
# or processors — so the obvious place to start working with parallelism is
# making use of those from within our Julia process.
#
# The first challenge, though, is knowing precisely how many "processors" you have.
# "Processors" is in scare quotes because, well, it's complicated.

versioninfo(verbose = true)

#-

#nb ;cat /proc/cpuinfo # on Linux machines

using Hwloc
Hwloc.num_physical_cores()

# What your computer reports as the number of processors might not be the same
# as the total number of "cores". While sometimes virtual processors can add
# performance, parallelizing a typical numerical computation over these virtual
# processors will lead to significantly worse performance because they still
# have to share much of the nuts and bolts of the computation hardware.

#-

# Julia is somewhat multithreaded by default! BLAS calls (like matrix multiplication) are
# already threaded:

using BenchmarkTools
A = rand(2000, 2000);
B = rand(2000, 2000);
@btime $A*$B;

# This is — by default — already using all your CPU cores! You can see the effect
# by changing the number of threads (which BLAS supports doing dynamically):

using LinearAlgebra
BLAS.set_num_threads(1)
@btime $A*$B
BLAS.set_num_threads(4)
@btime $A*$B

# ## What does it look like to implement your _own_ threaded algorithm?

#-

# Multithreading support is marked as "experimental" so far, but has been
# significantly improved in version 1.3.
# In the future, the core ideas will be the same, but it should become
# easier to use efficiently.

using .Threads

nthreads()

# Julia currently needs to start up knowing that it has threading support enabled.
#
# You do that with a environment variable. To get four threads, start Julia with:
#
# ```
# JULIA_NUM_THREADS=4 julia
# ```

#-

#nb # On JuliaBox, this is a challenge — we don't have access to the launching process!

#nb ;env JULIA_NUM_THREADS=4 julia -E 'using .Threads; nthreads()'

#-

threadid()

# So we're currently on thread 1. Of course a loop like this will
# just set the first element to one a number of times:

A = Array{Union{Int,Missing}}(missing, nthreads())
for i in 1:nthreads()
    A[threadid()] = threadid()
end
A

# But if we prefix it with `@threads` then the loop body runs on all threads!

@threads for i in 1:nthreads()
    A[threadid()] = threadid()
end
A

# So let's try implementing our first simple threaded algorithm — `sum`:

function threaded_sum1(A)
    r = zero(eltype(A))
    @threads for i in eachindex(A)
        @inbounds r += A[i]
    end
    return r
end

A = rand(10_000_000)
threaded_sum1(A)
@time threaded_sum1(A)

#-

sum(A)
@time sum(A)

# Whoa! What happened? Not only did we get the wrong answer, it was _slow_ to get it!

function threaded_sum2(A)
    r = Atomic{eltype(A)}(zero(eltype(A)))
    @threads for i in eachindex(A)
        @inbounds atomic_add!(r, A[i])
    end
    return r[]
end

threaded_sum2(A)
@time threaded_sum2(A)

# Alright! Now we got the correct answer (modulo some floating point associativity),
# but it's still slower than just doing the simple thing on 1 core.

threaded_sum2(A) ≈ sum(A)

# But it's still slow! Using atomics is much slower than just adding integers
# because we constantly have to go and check _which_ processor has the latest
# work! Also remember that each thread is running on its own processor — and
# that processor also supports SIMD!  Well, that is if it didn't need to worry
# about syncing up with the other processors...

function threaded_sum3(A)
    r = Atomic{eltype(A)}(zero(eltype(A)))
    len, rem = divrem(length(A), nthreads())
    @threads for t in 1:nthreads()
        rₜ = zero(eltype(A))
        @simd for i in (1:len) .+ (t-1)*len
            @inbounds rₜ += A[i]
        end
        atomic_add!(r, rₜ)
    end
    # catch up any stragglers
    result = r[]
    @simd for i in length(A)-rem+1:length(A)
        @inbounds result += A[i]
    end
    return result
end

threaded_sum3(A)
@time threaded_sum3(A)

# Dang, that's complicated. There's also a problem:

threaded_sum3(rand(10) .+ rand(10)im) # try an array of complex numbers!

# Isn't there an easier way?

R = zeros(eltype(A), nthreads())

#-

function threaded_sum4(A)
    R = zeros(eltype(A), nthreads())
    @threads for i in eachindex(A)
        @inbounds R[threadid()] += A[i]
    end
    r = zero(eltype(A))
    # sum the partial results from each thread
    for i in eachindex(R)
        @inbounds r += R[i]
    end
    return r
end

threaded_sum4(A)
@time threaded_sum4(A)

# This sacrifices our ability to `@simd` so it's a little slower, but at least we don't need to worry
# about all those indices! And we also don't need to worry about atomics and
# can again support arrays of any elements:

threaded_sum4(rand(10) .+ rand(10)im)

# ## Key takeaways from `threaded_sum`:
#
# * Beware shared state across threads — it may lead to wrong answers!
#     * Protect yourself by using atomics (or [locks/mutexes](https://docs.julialang.org/en/v1/base/multi-threading/#Synchronization-Primitives-1))
#     * Better yet: divide up the work manually such that the inner loops don't
#       share state. `@threads for i in 1:nthreads()` is a handy idiom.
#     * Alternatively, just use an array and only access a single thread's elements

#-

# # Beware of global state (even if it's not obvious!)
#
# Another class of algorithm that you may want to parallelize is a monte-carlo
# problem. Since each iteration is a new random draw, and since you're interested
# in looking at the aggregate result, this seems like it should lend itself to
# parallelism quite nicely!

using BenchmarkTools

#-

function serialpi(n)
    inside = 0
    for i in 1:n
        x, y = rand(), rand()
        inside += (x^2 + y^2 <= 1)
    end
    return 4 * inside / n
end
serialpi(1)
@time serialpi(100_000_000)

#-

using .Threads
nthreads()

# Let's use the techniques we learned above to make a fast threaded implementation:

function threadedpi(n)
    inside = zeros(Int, nthreads())
    @threads for i in 1:n
        x, y = rand(), rand()
        @inbounds inside[threadid()] += (x^2 + y^2 <= 1)
    end
    return 4 * sum(inside) / n
end
threadedpi(100_000_000)
@time threadedpi(100_000_000)

# Ok, now why didn't that work?  It's slow! Let's look at the sequence of random
# numbers that we generate:

import Random
Random.seed!(0)
N = 20000
Rserial = zeros(N)
for i in 1:N
    Rserial[i] = rand()
end
Rserial

#-

Random.seed!(0)
Rthreaded = zeros(N)
@threads for i in 1:N
    Rthreaded[i] = rand()
end
Rthreaded

#-

Set(Rserial) == Set(Rthreaded)

#-

indexin(Rserial, Rthreaded)

# Aha, `rand()` isn't threadsafe! It's mutating (and reading) some global each
# time to figure out what to get next. This leads to slowdowns — and worse — it
# skews the generated distribution of random numbers since some are repeated!!

const ThreadRNG = Vector{Random.MersenneTwister}(undef, nthreads())
@threads for i in 1:nthreads()
    ThreadRNG[Threads.threadid()] = Random.MersenneTwister()
end
function threadedpi2(n)
    inside = zeros(Int, nthreads())
    len, rem = divrem(n, nthreads())
    rem == 0 || error("use a multiple of $(nthreads()), please!")
    @threads for i in 1:nthreads()
        rng = ThreadRNG[threadid()]
        v = 0
        for j in 1:len
            x, y = rand(rng), rand(rng)
            v += (x^2 + y^2 <= 1)
        end
        inside[threadid()] = v
    end
    return 4 * sum(inside) / n
end
threadedpi2(10)
@time threadedpi2(100_000_000)

# As an aside, be careful about initializing many `MersenneTwister`s with
# different states. Better to use [`randjump`](https://docs.julialang.org/en/v1/manual/parallel-computing/#Side-effects-and-mutable-function-arguments-1) to skip ahead for a single state.

#-

# # Beware oversubscription
#
# Remember how BLAS is threaded by default? What happens if we try to `@threads`
# something that uses BLAS?

Ms = [rand(1000, 1000) for _ in 1:100]
function serial_matmul(As)
    first_idxs = zeros(length(As))
    for i in eachindex(As)
        @inbounds first_idxs[i] = (As[i]'*As[i])[1]
    end
    first_idxs
end
serial_matmul(Ms);
@time serial_matmul(Ms);

#-

using LinearAlgebra
BLAS.set_num_threads(nthreads()) # Explicitly tell BLAS to use the same number of threads
function threaded_matmul(As)
    first_idxs = zeros(length(As))
    @threads for i in eachindex(As)
        @inbounds first_idxs[i] = (As[i]'*As[i])[1]
    end
    first_idxs
end
threaded_matmul(Ms)
@time threaded_matmul(Ms);

#-

BLAS.set_num_threads(1)
@time threaded_matmul(Ms);

#-

@time serial_matmul(Ms) # Again, now that BLAS has just 1 thread

# # Beware "false sharing"

#-

# Remember the memory latency table?
#
#
# | System Event                   | Actual Latency | Scaled Latency |                          |
# | ------------------------------ | -------------- | -------------- | ------------------------ |
# | One CPU cycle                  |     0.4 ns     |     1 s        | ← work happens here     |
# | Level 1 cache access           |     0.9 ns     |     2 s        |                          |
# | Level 2 cache access           |     2.8 ns     |     7 s        |                          |
# | Level 3 cache access           |      28 ns     |     1 min      |                          |
# | Main memory access (DDR DIMM)  |    ~100 ns     |     4 min      | ← we have control here  |
#
# This is what a typical modern cpu looks like:
#
# ![Intel Core i7](https://raw.githubusercontent.com/JuliaComputing/JuliaAcademyData.jl/master/courses/Parallel_Computing/images/i7.jpg)
#
# Multiple cores on the same processor share the L3 cache, but do not share L1 and L2 caches! So what happens if we're accessing and mutating data from the same array across multiple cores?
#
# ![Cache coherency](https://raw.githubusercontent.com/JuliaComputing/JuliaAcademyData.jl/master/courses/Parallel_Computing/images/false-sharing.gif)
#
# Unlike "true" sharing — which we saw above — false sharing will still return the correct answer! But it does so at the cost of performance. The cores recognize they don't have exclusive access to the cache line and so upon modification they alert all other cores to invalidate and re-fetch the data.
#
# ```julia
# function threaded_sum4(A)
#     R = zeros(eltype(A), nthreads())
#     @threads for i in eachindex(A)
#         @inbounds R[threadid()] += A[i]
#     end
#     r = zero(eltype(A))
#     # sum the partial results from each thread
#     for i in eachindex(R)
#         @inbounds r += R[i]
#     end
#     return r
# end
# ```

#-

# ## New in 1.3: parallel @spawn

# Julia 1.3 adds Threads.@spawn, which runs a "task" in parallel
# on any available thread, using a work queue.
# This addresses the problem of having library functions implemented with
# `@threads` and then having callers call them with `@threads`.

#-

# # Threading takeaways:
#
# * It's easy! Just start Julia with `JULIA_NUM_THREADS` and tack a `@threads` on your loop
# * Well, not so fast
#     * Be aware of your hardware to set `JULIA_NUM_THREADS` appropiately
#     * Beware shared state (for both performance and correctness)
#     * Beware global state (even if it's not obvious)
#     * Beware false sharing (if Julia/LLVM don't handle it for you)
# * We need to think carefully about how to design parallel algorithms!

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
@testset "Secondary code or bit detector" begin
    detector = SecondaryCodeOrBitDetector()

    @test Tracking.get_buffer(detector) == 0
    @test Tracking.length(detector) == 0
    @test Tracking.found(detector) == false

    next_detector = Tracking.find(GPSL1, detector, -1.0 + 0.0im)
    @test Tracking.get_buffer(next_detector) == 0
    @test Tracking.length(next_detector) == 1
    @test Tracking.found(next_detector) == false

    next2_detector = Tracking.find(GPSL1, next_detector, 1.0 + 0.0im)
    @test Tracking.get_buffer(next2_detector) == 1
    @test Tracking.length(next2_detector) == 2
    @test Tracking.found(next2_detector) == false

    next3_detector = Tracking.find(GPSL1, next2_detector, 1.0 + 0.0im)
    @test Tracking.get_buffer(next3_detector) == 3
    @test Tracking.length(next3_detector) == 3
    @test Tracking.found(next3_detector) == false

    detector = SecondaryCodeOrBitDetector(1, 1, true)

    @test Tracking.get_buffer(detector) == 1
    @test Tracking.length(detector) == 1
    @test Tracking.found(detector) == true

    next_detector = Tracking.find(GPSL1, detector, 1.0 + 0.0im)
    @test Tracking.get_buffer(next_detector) == 1
    @test Tracking.length(next_detector) == 1
    @test Tracking.found(next_detector) == true

    detector = SecondaryCodeOrBitDetector(0x7ffff80000, 39, false)
    next_detector = Tracking.find(GPSL1, detector, -1.0 + 0.0im)
    @test Tracking.get_buffer(next_detector) == 0x000000fffff00000
    @test Tracking.length(next_detector) == 40
    @test Tracking.found(next_detector) == true
end


text:
function make_directory(path; delete_contents_if_existing = true)

    """makes a directory and goes inside, removes the previous contents of a
    similarly named directory if existing"""

    try
        if delete_contents_if_existing
            rm(path, recursive = true)
        else
            pass()
        end
    catch err

    finally
        mkpath(path)
        cd(path)
    end
end

function print_error(err)

	"""prints the error in a try/catch"""
	buf = IOBuffer()
    showerror(buf, err)
    message = String(take!(buf))
end

# to save julia variables : using JLD2

function save_matrix(mat, filename = "matrix.csv")

	CSV.write(filename,  DataFrame(mat), header=false)

end


text:
using BEAST

Δt, Nt = 1.0, 200
T0 = timebasisshiftedlagrange(Δt, Nt, 0)
T1 = timebasisshiftedlagrange(Δt, Nt, 1)
T2 = timebasisshiftedlagrange(Δt, Nt, 2)
T3 = timebasisshiftedlagrange(Δt, Nt, 3)
iT0 = integrate(T0)
iT1 = integrate(T1)
iT2 = integrate(T2)
iT3 = integrate(T3)
δ = timebasisdelta(Δt, Nt)

Id = BEAST.Identity()
Z = assemble(Id, δ, iT1)
Z = reshape(Z,1,1,length(Z))

taxis = range(0, step=Δt, length=Nt)
width, delay = 2.0, 6.0
f(t) = -2*(t-delay)/width^2*exp(-((t-delay)/width)^2)
F(t) = exp(-((t-delay)/width)^2)

b = F.(taxis)'
u = marchonintime(inv(Z[:,:,1]), Z, b, Nt)

plot(taxis, u[1,:])

taxis2 = range(-3Δt, 5Δt, length=200)
plot()
plot!(taxis2, iT3.(taxis2))

kmax = 2
cv = BEAST.ConvOp(
    reshape(Z[1,1,1:kmax],kmax,1,1),
    fill(1,(1,1)),
    fill(kmax,(1,1)),
    fill(1.0,(1,1)),
    Nt)

w, Q = BEAST.polyeig(cv)

u2 = marchonintime(inv(Z[:,:,1]), cv, b, Nt)


text:
# Stalin Sort

"""
A Julia implementation of the Stalin Sort algorithm.
Inspired by the Python implementation.

Elements that are out of order are eliminated.

Args:
    arr: An AbstractVector

Returns:
    Stalin sorted array.
"""

function stalinsort(arr::AbstractVector)

    if isempty(arr)
        return []
    end

    max_val = arr[1]

    function add_val(num)
        max_val = num
        return num
    end

    return [add_val(x) for x in arr if x >= max_val]
end

# Test examples

#  unordered: [1, 2, 3, 4, 5]
# -> ordered: [1, 2, 3, 4, 5]
test1 = [1, 2, 3, 4, 5]
println("# unordered: $test1\n -> ordered: ", stalinsort(test1))

#  unordered: [10, 9, 8, 7, 2]
# -> ordered: [10]
test2 = [10, 9, 8, 7, 2]
println("# unordered: $test2\n -> ordered: ", stalinsort(test2))

# unordered: [1, 2, 3, 5, 4, 6]
#-> ordered: [1, 2, 3, 5, 6]
test3 = [1, 2, 3, 5, 4, 6]
println("# unordered: $test3\n -> ordered: ", stalinsort(test3))

#  unordered: [1.0, 2.4, 2.2, 3.0, 5.0, 4.0, 6.0]
# -> ordered: [1.0, 2.4, 3.0, 5.0, 6.0]
test4 = [1, 2.4, 2.2, 3, 5, 4, 6]
println("# unordered: $test4\n -> ordered: ", stalinsort(test4))

#  unordered: [-9, 6, 4, 10, -10, -5, 1, 8, -5, 10]
# -> ordered: [-9, 6, 10, 10]
test5 = rand(-10:10, 10)
println("# unordered: $test5\n -> ordered: ", stalinsort(test5))

#  unordered: [1, 2, 2, 5, 4, 5, 3, 7]
# -> ordered: [1, 2, 2, 5, 5, 7]
test6 = [1, 2, 2, 5, 4, 5, 3, 7]
println("# unordered: $test6\n -> ordered: ", stalinsort(test6))

max_stars_repo_name: twofist/stalin-sort

max_issues_repo_name: twofist/stalin-sort

max_forks_repo_name: twofist/stalin-sort


text:
# Headcount.jl

###### headcount #####
"""
    headcount(v, z)

Compute the Headcount Ratio of a vector `v` at a specified poverty threshold `z`.

# Examples
```julia
julia> using Inequality
julia> headcount([8, 5, 1, 3, 5, 6, 7, 6, 3], 4)
0.3333333333333333
```
"""
headcount(v::AbstractVector{<:Real}, z::Real)::Float64 = length(v[v .< z]) / length(v)

###### weighted headcount #####
"""
    headcount(v, w, z)

Compute the Headcount Ratio of a vector `v` at a specified poverty threshold `z`,
using weights given by a weight vector `w`.

Weights must not be negative, missing or NaN. The weights and data vectors must have the same length.

# Examples
```julia
julia> using Inequality
julia> headcount([8, 5, 1, 3, 5, 6, 7, 6, 3], [0.1,0.5,0.3,0.8,0.1,0.5,0.3,0.8,0.2], 4)
0.36111111111111116
```
"""
function headcount(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}, z::Real)::Float64

    checks_weights(v, w)

    return sum(w[v .< z]) / sum(w)
end

function headcount(v::AbstractVector{<:Real}, w::AbstractWeights, z::Real)::Float64

    checks_weights(v, w)

    return sum(w[v .< z]) / w.sum

end

wheadcount(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}, z::Real)::Float64 = headcount(v, w, z)

wheadcount(v::AbstractVector{<:Real}, w::AbstractWeights, z::Real)::Float64 = headcount(v, w, z)


text:
module StructMapping

export convertdict

## Helper functions
_convertdict(::Type, x) = x
_convertdict(::Type{Union{T, Nothing}}, x) where T = _convertdict(T, x)
_convertdict(::Type{Union{T, Nothing}}, d::AbstractDict) where T = _convertdict(T, d)
_convertdict(::Type{T}, v::AbstractVector) where T <: AbstractVector = _convertdict.(eltype(T), v)
_convertdict(::Type{T}, d::AbstractDict) where T <: AbstractDict = d

function _convertdict(T::Type, d::AbstractDict)
    kwargs = Dict{Symbol, Any}()
    for (k, v) in pairs(d)
        symk = Symbol(k)
        kwargs[symk] = _convertdict(fieldtype(T, symk), v)
    end
    return T(;kwargs...)
end

"""
    convertdict(T::Type, d::AbstractDict)

Convert the given dictionary to a object of `T`.
`T` must be decorated with `@with_kw` or `@with_kw_noshow` of Parameters.jl.
"""
convertdict(T::Type, d::AbstractDict) = _convertdict(T, d)

end # module


text:
using NCDatasets
using Formatting

using Statistics

if length(ARGS) >= 1
    fn = ARGS[1]
end

Re  = 6371000.0 #  m
cp  = 1004.0    #  J / kg K
g   = 9.81      #  m / s^2
Lw  = 2.26e6    #  J / kg

ds = Dataset(fn, "r")

lat = ds["lat"][:]
lon = ds["lon"][:]
lev = ds["lev"][:]
ilev = ds["ilev"][:]
rec = ds["time"][:]

deg2rad = π/180.0

#cross_section = (lon[2] - lon[1]) * deg2rad * Re * cos.(lat * deg2rad)

# mass
Δp = (ilev[2:end] - ilev[1:end-1]) * 100.0

function calMeanTransport(v)

    # Get the zonal mean
    v = mean(v, dims=(1,))[1, :, :, :]

    # Weight with mass
    for j=1:length(lat), k=1:length(lev), l=1:length(rec)
        v[j, k, l] *= Δp[k] * cos(lat[j] * deg2rad)
    end

    v .*= 2 * π * Re / g

    # Get the vertical integration and time mean
    v = sum(v, dims=(2,))[:, 1, :]
    v = mean(v, dims=(2,))[:, 1]

    return v
end

# Meridional Static Energy Transport

MHT = cp * ds["VT"][:]  |> calMeanTransport
MGT = ds["VZ"][:]       |> calMeanTransport
MLT = Lw * ds["VQ"][:]  |> calMeanTransport

MMSET = MHT + MGT + MLT
close(ds)

print("Import Pyplot... ")
using PyPlot
println("Done. ")

fig, ax = plt[:subplots](1, 1, figsize=(12,8))

ax[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc", dashes=(5,3))
ax[:plot](lat, MMSET / 1e15, linewidth=2, "k-",  label="Moist static energy transport")
ax[:plot](lat, MHT   / 1e15, linewidth=2, "k--", label="Heat energy transport")
ax[:plot](lat, MGT   / 1e15, linewidth=2, "r--", label="Geopotential energy transport")
ax[:plot](lat, MLT   / 1e15, linewidth=2, "b--", label="Liquid energy transport")

ax[:legend]()
ax[:set_xlabel]("Latitude [deg]")
ax[:set_ylabel]("[PW]")

ax[:set_title](fn)
ax[:set_xlim](-90, 90)
ax[:set_ylim](-15, 15)

fig[:savefig](fn * ".png", dpi=200)
plt[:close]()

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
using Arrow
using DataFrames
using DrWatson
using MAT

"""

`fdt` files are a type of binary file format that stores "float data".

"""
fdt_parser(path, mat) =
    open(path) do io
        read!(
            io,
            Array{Float32}(undef, (Int(mat["EEG"]["nbchan"]), Int(mat["EEG"]["pnts"]))),
        )
    end

# Reference: https://stackoverflow.com/questions/20484581/search-for-files-in-a-folder
searchdir(path, key) = filter(x -> occursin(key, x), readdir(path))

function fdt_to_arrow(file_path)

    subjects = searchdir(file_path, "sub")

    for sub_num in subjects
        sessions = readdir(file_path * "/" * sub_num)
        for ses_num in sessions
            eeg_sessions = readdir(file_path * "/" * sub_num * "/" * ses_num)
            for eeg_num in eeg_sessions
                files = searchdir(
                    file_path * "/" * sub_num * "/" * ses_num * "/" * eeg_num,
                    ".fdt",
                )
                headers = searchdir(
                    file_path * "/" * sub_num * "/" * ses_num * "/" * eeg_num,
                    ".set",
                )
                for (file, header) in zip(files, headers)
                    fdt_file =
                        file_path *
                        "/" *
                        sub_num *
                        "/" *
                        ses_num *
                        "/" *
                        eeg_num *
                        "/" *
                        file
                    header_file =
                        file_path *
                        "/" *
                        sub_num *
                        "/" *
                        ses_num *
                        "/" *
                        eeg_num *
                        "/" *
                        header
                    header_data = matread(header_file)
                    fdt_data = transpose(fdt_parser(fdt_file, header_data))
                    Arrow.write(
                        datadir(
                            "exp_pro",
                            "$sub_num",
                            "$ses_num",
                            "$eeg_num",
                            "$(file[1:end - 3])arrow",
                        ),
                        DataFrame(fdt_data),
                    )
                end
            end
        end
    end
end

function move_files(file_path, file_type)
    subjects = searchdir(file_path, "sub")
    for sub_num in subjects
        sessions = readdir(file_path * "/" * sub_num)
        for ses_num in sessions
            eeg_sessions = readdir(file_path * "/" * sub_num * "/" * ses_num)
            for eeg_num in eeg_sessions
                files = searchdir(
                    file_path * "/" * sub_num * "/" * ses_num * "/" * eeg_num,
                    "$file_type",
                )
                for file in files
                    mv_file =
                        file_path *
                        "/" *
                        sub_num *
                        "/" *
                        ses_num *
                        "/" *
                        eeg_num *
                        "/" *
                        file
                    dest_path = "data/exp_pro/" * sub_num * "/" * ses_num * "/" *  eeg_num * "/" * file
                    cp(mv_file, dest_path)
                end
            end
        end
    end
end

@quickactivate "NeuriViz"

# fdt_to_arrow(datadir("exp_raw"))

max_stars_repo_name: TheCedarPrince/NeuriViz

max_issues_repo_name: TheCedarPrince/NeuriViz

max_forks_repo_name: TheCedarPrince/NeuriViz


text:
abstract type GebAction end
struct CrossOver <: GebAction end
struct Fight <: GebAction end
struct TurnClockwise <: GebAction end
struct TurnAntiClockwise <: GebAction end
struct MoveForward <: GebAction end
struct NoAction <: GebAction end

_get_prefix(::CrossOver) = "01"
_get_prefix(::Fight) = "100"
_get_prefix(::TurnAntiClockwise) = "101"
_get_prefix(::TurnClockwise) = "110"
_get_prefix(::MoveForward) = "111"

function perform!(model::GebModel, organism::GebOrganism)
    perform!(CrossOver(), model, organism)
    perform!(Fight(), model, organism)
    perform!(MoveForward(), model, organism)
    perform!(TurnAntiClockwise(), model, organism)
    perform!(TurnClockwise(), model, organism)
end

function perform!(action::GebAction, model::GebModel, organism::GebOrganism)
    prefix = _get_prefix(action)

    excitatory_sum = 0
    has_match = false

    for external in organism.network.outputs
        if _is_match(prefix, external.string)
            excitatory_sum += external.io_value
            has_match = true
        end
    end

    if has_match
        excitatory_sum = min(excitatory_sum, 1)
    else
        excitatory_sum = NOISE_LEVEL_MISSING_OUTPUT*rand()
    end

    perform!(action, model, organism, excitatory_sum)
end

function perform!(::CrossOver, model::GebModel, parent_1::GebOrganism, activation)
    if activation < CROSS_OVER_THRESHOLD return end

    parent_2 = _get_organism_in_front(model, parent_1)

    if parent_2 === nothing return end

    genotype_1 = parent_1.genotype
    genotype_2 = parent_2.genotype

    min_size = min(length(genotype_1), length(genotype_2))

    cut_point_1 = rand(1:max(1, min_size - 1)) # there is never a cut off at the end

    if length(genotype_2) == 1
        cut_point_2 = 1
    elseif cut_point_1 == 1
        cut_point_2 = 1
    elseif cut_point_1 == length(genotype_2)
        cut_point_2 = cut_point_1
    else
        cut_point_2 = cut_point_1 + rand([0, 2])
    end

    child_genotype = genotype_1[1:cut_point_1] * genotype_2[cut_point_2:end]
    child_genotype = get_mutated_genotype(child_genotype)

    if _get_organism_behind(model, parent_1) === nothing
        x_child_grid, y_child_grid = _get_cell_in_direction(model, parent_1.coordinates, 360 - parent_1.direction)

        child = GebOrganism(child_genotype,
                            (x_child_grid - 0.5, y_child_grid - 0.5),
                            parent_1.genotype,
                            parent_2.genotype,
                            model.time)

        push!(parent_1.daughters, child)
        push!(parent_2.daughters, child)

        add_organism!(model, child)

        log_birth(model.logger, model, child, [parent_1, parent_2])
    else
        child = GebOrganism(child_genotype,
        parent_2.coordinates,
        parent_1.genotype,
        parent_2.genotype,
        model.time)

        push!(parent_1.daughters, child)
        push!(parent_2.daughters, child)

        kill!(model, parent_2)

        add_organism!(model, child)

        log_birth(model.logger, model, child, [parent_1, parent_2])
    end
end

function perform!(::Fight, model::GebModel, fighter::GebOrganism, activation)
    if activation < FIGHT_THRESHOLD return end

    target = _get_organism_in_front(model, fighter)

    if target !== nothing
        kill!(model, target)
    end
end

function perform!(::TurnAntiClockwise, ::GebModel, organism::GebOrganism, angle)
    organism.direction = mod(organism.direction - min(MAX_TURN, Int(round(MAX_TURN*angle))), 360)
end
function perform!(::TurnClockwise, ::GebModel, organism::GebOrganism, angle)
    organism.direction = mod(organism.direction + min(MAX_TURN, Int(round(MAX_TURN*angle))), 360)
end

function perform!(::MoveForward, model::GebModel, organism::GebOrganism, distance)
    x_start, y_start = organism.coordinates
    x_start_grid, y_start_grid = _get_grid_coordinates((x_start, y_start))
    dir = organism.direction

    x_end, y_end = x_start + cosd(dir)*distance, y_start + sind(dir)*distance
    x_end, y_end = mod(x_end, model.size), mod(y_end, model.size)

    x_end_grid, y_end_grid = _get_grid_coordinates((x_end, y_end))

    # has the cell even changed?
    if (x_start_grid, y_start_grid) == (x_end_grid, y_end_grid)
        organism.coordinates = x_end, y_end
        return
    end

    # is there already an organism on the new cell?
    if model.grid[x_end_grid, y_end_grid] !== nothing
        return
    end

    # move to new cell
    organism.coordinates = x_end, y_end
    model.grid[x_start_grid, y_start_grid] = nothing
    model.grid[x_end_grid, y_end_grid] = organism
end

function perform!(::NoAction, ::GebModel, ::GebOrganism, _) end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# Integration test combining AtomsBase, NBodySimulator, and DFTK

@testset "Integration: AB & NBS & DFTK" begin
    # TODO
end


text:
struct TerminalPresenter <: Presenter
end
struct BrowserPresenter <: Presenter
end
@enum YesNoCancel yes=1 no=0 cancel=-1

# TODO: use https://juliagizmos.github.io/Interact.jl/latest/widgets.html#InteractBase.checkboxes
# to display ingredients and let contents of barcabinet be indicated by checked checkboxes. Possibly use for selecting ingredients to include in recipe
# TODO: On checkbox, update barcabinet
# TODO: Button for Add: Recipe, Ingredient
# TODO: Interface for updateing all fields of Recipe
# TODO: Menu items for save and load files

# Presenter interface specification =====================
error(p::Presenter, msg) = error("error not defined for ::$(typeof(p))")
info(p::Presenter, msg) = error("info not defined for ::$(typeof(p))")
prompt(p::Presenter, msg)::String = error("prompt not defined for ::$(typeof(p))")
prompt(p::Presenter, ::Type{YesNoCancel}, msg)::YesNoCancel = error("yesnocancel not defined for ::$(typeof(p))")

# TerminalPresenter implementation
error(p::TerminalPresenter, msg) = @error(msg)
info(p::TerminalPresenter, msg)  = @info(msg)

function prompt(p::TerminalPresenter, msg)::String
    println(msg)
    readline(stdin)
end

function prompt(p::TerminalPresenter, ::Type{YesNoCancel}, msg)::YesNoCancel
    println(msg)
    print("yes/no/cancel (y/n/c)? ")
    r = lowercase(readline(stdin)[1])
    r == 'y' ? yes : r == 'n' ? no : cancel
end


text:
# If not already installed, first hit "]" within a Julia REPL. Then type:
# add DifferentialEquations DiffEqBiological PyPlot Plots Latexify

using DifferentialEquations, DiffEqBiological, Plots, Latexify
pyplot(fmt=:svg);

repressilator = @reaction_network begin
    hillr(P₃,α,K,n), ∅ --> m₁
    hillr(P₁,α,K,n), ∅ --> m₂
    hillr(P₂,α,K,n), ∅ --> m₃
    (δ,γ), m₁ ↔ ∅
    (δ,γ), m₂ ↔ ∅
    (δ,γ), m₃ ↔ ∅
    β, m₁ --> m₁ + P₁
    β, m₂ --> m₂ + P₂
    β, m₃ --> m₃ + P₃
    μ, P₁ --> ∅
    μ, P₂ --> ∅
    μ, P₃ --> ∅
end α K n δ γ β μ;

latexify(repressilator; env=:chemical)

mathjax = WEAVE_ARGS[:doctype] == "pdf" ? false : true
x = latexify(repressilator; env=:chemical, starred=true, mathjax=mathjax);
display("text/latex", "$x");

latexify(repressilator, cdot=false)

x = latexify(repressilator, cdot=false, starred=true);
display("text/latex", "$x");

speciesmap(repressilator)

paramsmap(repressilator)

# parameters [α,K,n,δ,γ,β,μ]
p = (.5, 40, 2, log(2)/120, 5e-3, 20*log(2)/120, log(2)/60)

# initial condition [m₁,m₂,m₃,P₁,P₂,P₃]
u₀ = [0.,0.,0.,20.,0.,0.]

# time interval to solve on
tspan = (0., 10000.)

# create the ODEProblem we want to solve
oprob = ODEProblem(repressilator, u₀, tspan, p)

sol = solve(oprob, saveat=10.)
plot(sol, fmt=:svg)

# first we redefine the initial condition to be integer valued
u₀ = [0,0,0,20,0,0]

# next we create a discrete problem to encode that our species are integer valued:
dprob = DiscreteProblem(repressilator, u₀, tspan, p)

# now we create a JumpProblem, and specify Gillespie's Direct Method as the solver:
jprob = JumpProblem(dprob, Direct(), repressilator, save_positions=(false,false))

# now let's solve and plot the jump process:
sol = solve(jprob, SSAStepper(), saveat=10.)
plot(sol, fmt=:svg)

rjs = regularjumps(repressilator)
lprob = JumpProblem(dprob, Direct(), rjs)
lsol = solve(lprob, SimpleTauLeaping(), dt=.1)
plot(lsol, plotdensity=1000, fmt=:svg)

bdp = @reaction_network begin
  c₁, X --> 2X
  c₂, X --> 0
  c₃, 0 --> X
end c₁ c₂ c₃
p = (1.0,2.0,50.)
u₀ = [5.]
tspan = (0.,4.);

latexify(bdp, noise=true, cdot=false)

x = latexify(bdp, noise=true, cdot=false, starred=true);
display("text/latex", "$x");

# SDEProblem for CLE
sprob = SDEProblem(bdp, u₀, tspan, p)

# solve and plot, tstops is used to specify enough points
# that the plot looks well-resolved
sol = solve(sprob, tstops=range(0., step=4e-3, length=1001))
plot(sol, fmt=:svg)

latexify(jacobianexprs(repressilator), cdot=false)

x = latexify(jacobianexprs(repressilator), cdot=false, starred=true);
display("text/latex", "$x");

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file], remove_homedir=true)

max_stars_repo_name: SciML/SciMLTutorialsOutput

max_issues_repo_name: SciML/SciMLTutorialsOutput

max_forks_repo_name: SciML/SciMLTutorialsOutput


text:
function fibonacci_sum(limit)
    sequence = [1, 2]; total = 0
    while sequence[end] <= limit
        if iseven(sequence[end])
            total += sequence[end]
        end
        push!(sequence, sequence[end] + sequence[end-1])
    end
    println(total)
end

fibonacci_sum(4000000)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
using Random
include("tools.jl")
include("pam250scoring.jl")

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

# creates a number of gaps to fill a certain string length randomly
# the gaps are index positions in sorted order as the function it
# serves is to strings with gaps to a certain length depending on missing
# items, like 'aa' and 'bbb' going to size 4 would need 2 gaps for 'aa' and
# 1 gap for 'bbb', but without duplication it is fine to assume insertions
# because for 'aa', we could only at worst get (3, 4) for size 4, so putting
# incrementally adding the gaps works by placing at end or inserting into
# existing index, benefit of this one is it ends sorted
function create_random_more_gaps(number_of_gaps, string_length)
    # fill [1, 2, 3, ..., string_length]
    my_range = collect(Int64, 1:string_length)

    # delete items randomly until number of gaps remain
    for i = 1:string_length - number_of_gaps
        deleteat!(my_range, Int(rand(1:string_length)))
        string_length -= 1
    end

    return my_range
end

# same as above, but puts gaps into a new list, instead of deleting from
# until gaps remain, above is bad when list is large and gaps are small
# this scales better for very large strings, assuming portion of gaps is
# not large as well, but this ends unsorted which is a downside - added sort
function create_random_less_gaps(number_of_gaps, string_length)
    # fill [1, 2, 3, ..., string_length]
    my_range = collect(Int64, 1:string_length)
    gap_list = []

    # delete items randomly until number of gaps remain
    for i = 1:number_of_gaps
        index = Int(rand(1:string_length))
        push!(gap_list, my_range[index])
        deleteat!(my_range, index)
        string_length -= 1
    end

    return sort!(gap_list)
end

# score using the gap list and the sequence
function score_gap_list(input_sequence_set, gap_list, sequence_length, score_list)
    score = 0

    num_items = length(input_sequence_set)

    # start offsets for accessing with respect to current index,
    # ex: [AAA] w/ gap [2,3] = [A__AA] - a 5 item string
    # but instead of building string, we want to look at A, then gap of 2, 3, then AAA
    # so index will be 1-5, when needing second A, which should be index 4, offsets will be
    # +2 gap_list, -2 sequence_list, so index of 4 will give 2 into sequence_list and show
    # that we've looked past the end of the gap list (because we already had 2 gaps)
    gap_offset_list = fill(0, num_items)

    # make the last item an overflow position
    for i = 1:num_items
        push!(gap_list[i], 0)
    end

    # save variable allocation later
    score_buffer = 0
    item_one = nothing
    item_two = nothing

    # go through sequence length, to index by column, all input sequences + gaps should = sequence length
    for i = 1:sequence_length
        score_buffer = 0

        # go through each input sequence, which is a single line of characters
        # last one will just pass through, either not doing a score or it was a gap
        # so that gap count and offset still needs to increase
        for j = 1:num_items
            # get either a gap or the character
            # if it is a gap, this layer is all gap gap_penalty
            # so save a lot of calculation by mult and continue
            if i == gap_list[j][1 + gap_offset_list[j]]
                gap_offset_list[j] += 1
                item_one = '_'
            else
                item_one = input_sequence_set[j][i - gap_offset_list[j]]
            end

            # go through each input seqence from the next item up to the last
            # final pattern of access = [1, 2, 3] => 12, 13, 23
            for k = j + 1:num_items
                # get either a blank space or the character
                if i == gap_list[k][1 + gap_offset_list[k]]
                    item_two = '_'
                else
                    item_two = input_sequence_set[k][i - gap_offset_list[k]]
                end

                if item_one == '_' || item_two == '_'
                    if item_one != item_two
                        score_buffer += score_list[3]
                    end
                    continue
                end

                if (item_one == item_two)
                    score_buffer += score_list[1]
                else
                    score_buffer += score_list[2]
                end
            end
        end

        score += score_buffer
    end

    # remove the added overflow position
    for i = 1:num_items
        pop!(gap_list[i])
    end

    return score
end

# "save" some of the population for next generation, optionally return a modified input list
function elitism(chrom_pop_fitness, proportion, children_cap)
    num_items = Int(ceil(length(chrom_pop_fitness) * proportion))

    # return the leftover as in the paper
    #return [chrom_pop_fitness[1:num_items], chrom_pop_fitness[num_items + 1:length(chrom_pop_fitness)]]

    # returns the original list back, better results
    return [chrom_pop_fitness[1:min(children_cap, num_items)], chrom_pop_fitness]
end

# take the population and create children off of random points left and right combined
# against the rest of the population by an amount num_best
# ex: 5 would take the first 5 and make children with the rest
# if a cap is specified, stops when enough accepted children are added
# children only accepted if higher score than highest or lowest parent (design decision)
# original population dies off, only children returned
function crossover(input_sequence, population, num_best, cap, len, score_list, score_type, criteria, return_early, gap_ext_penalty)
    ret_children = []
    input_length = length(input_sequence)

    for i = 1:num_best
        for j = i + 1:length(population)
            cut = Int(rand(1:input_length - 1))

            child1 = []
            child2 = []
            append!(child1, population[i][2][1:cut])
            append!(child1, population[j][2][cut + 1:input_length])
            append!(child2, population[j][2][1:cut])
            append!(child2, population[i][2][cut + 1:input_length])

            if score_type == 2
                child1_w_fitness = [ score_gap_list_pam250(input_sequence, child1, len, score_list, gap_ext_penalty), child1 ]
                child2_w_fitness = [ score_gap_list_pam250(input_sequence, child2, len, score_list, gap_ext_penalty), child2 ]
            else
                child1_w_fitness = [ score_gap_list(input_sequence, child1, len, score_list), child1 ]
                child2_w_fitness = [ score_gap_list(input_sequence, child2, len, score_list), child2 ]
            end

            if criteria == 1
                # average of the 2
                fitness_criteria = (population[i][1] + population[j][1]) / 2
            elseif criteria == 2
                # score of worst of 2 parents
                fitness_criteria = min(population[i][1], population[j][1])
            else
                # score of best of 2 parents
                fitness_criteria = max(population[i][1], population[j][1])
            end

            if (child1_w_fitness[1] > fitness_criteria)
                push!(ret_children, child1_w_fitness)
            end

            if (child2_w_fitness[1] > fitness_criteria)
                push!(ret_children, child2_w_fitness)
            end

            # not exact, because 2 children are made but this check is done once
            if !isnothing(cap) && length(ret_children) > cap
                if return_early
                    return ret_children
                end

                sort!(ret_children, by = x -> (x[1]), rev = true)
                pop!(ret_children)
                pop!(ret_children)
            end
        end
    end
    return ret_children
end

# the difference to crossover 1 is this picks one or the other randomly
# for each cut, instead of left right on all
function crossover2(input_sequence, population, num_best, cap, len, score_list, score_type, criteria, return_early, gap_ext_penalty)
    ret_children = []
    input_length = length(input_sequence)

    for i = 1:num_best
        for j = i + 1:length(population)
            child1 = []
            child2 = []

            current = [ i, j ]

            for k = 1:input_length
                shuffle!(current)
                push!(child1, population[current[1]][2][k])
                push!(child2, population[current[2]][2][k])
            end

            if score_type == 2
                child1_w_fitness = [ score_gap_list_pam250(input_sequence, child1, len, score_list, gap_ext_penalty), child1 ]
                child2_w_fitness = [ score_gap_list_pam250(input_sequence, child2, len, score_list, gap_ext_penalty), child2 ]
            else
                child1_w_fitness = [ score_gap_list(input_sequence, child1, len, score_list), child1 ]
                child2_w_fitness = [ score_gap_list(input_sequence, child2, len, score_list), child2 ]
            end

            if criteria == 1
                # average of the 2
                fitness_criteria = (population[i][1] + population[j][1]) / 2
            elseif criteria == 2
                # score of worst of 2 parents
                fitness_criteria = min(population[i][1], population[j][1])
            else
                # score of best of 2 parents
                fitness_criteria = max(population[i][1], population[j][1])
            end

            if (child1_w_fitness[1] > fitness_criteria)
                push!(ret_children, child1_w_fitness)
            end

            if (child2_w_fitness[1] > fitness_criteria)
                push!(ret_children, child2_w_fitness)
            end

            # not exact, because 2 children are made but this check is done once
            if !isnothing(cap) && length(ret_children) > cap
                if return_early
                    return ret_children
                end

                sort!(ret_children, by = x -> (x[1]), rev = true)
                pop!(ret_children)
                pop!(ret_children)
            end
        end
    end
    return ret_children
end

# needs input sequence for rescoring, children population, and a mutation string_length
# as well as segment length to pick the random number for each new mutation item
function mutate(input_sequence, children_population, mutation_strength, mutation_chance, segment_length, score_list, score_type, gap_ext_penalty)
    index_list_mutation = []
    rand_index = 0

    # for every child
    for i = 1:length(children_population)
        # only half the time do a small mutation
        # otherwise keep it
        if rand() < mutation_chance
            continue
        end

        # for every gap list in every child
        for j = 1:length(children_population[i])
            randsubseq!(index_list_mutation, 1:length(children_population[i][j]), mutation_strength)

            for k = 1:length(index_list_mutation)
                rand_index = Int(rand(1:segment_length))
                while rand_index in children_population[i][j]
                    rand_index = Int(rand(1:segment_length))
                end

                children_population[i][j][index_list_mutation[k]] = rand_index
            end
            sort!(children_population[i][j])
        end
    end

    children_population_scored = get_fitness_population(input_sequence, children_population, segment_length, score_list, score_type, gap_ext_penalty)

    return children_population_scored
end

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

# create a number of gap lists representing a chromosome population
# that will work by taking the list of number of gaps to needed
# and the length of the final string
# ex: [ AAA, BBBB, C ] = lengths [ 3, 4, 1 ]
# ex: length = 5, then we would get counts of [2, 1, 4]
# as many times as population size
function create_gap_population(gap_counts, final_length, population_size)
    population = []

    for i = 1:population_size
        push!(population, create_gap_list(gap_counts, final_length))
    end

    return population
end

# take a list of gaps ex: [2, 1, 4]
# give a final length ex: 5
# create a number of gaps randomly in the range up to length
# that would pad the current length to that amount
# ex: [ [1, 4], [3], [1, 2, 4, 5] ]
function create_gap_list(gap_counts, final_length)
    median = final_length / 2

    gap_list = []
    # go through each number of gaps to create
    for i = 1:length(gap_counts)
        # when large number of gaps compared to final size, do first, otherwise do second
        if gap_counts[i] > median
            # deletes from index list of final_length until number of gaps remain
            # ends sorted for us
            push!(gap_list, create_random_more_gaps(gap_counts[i], final_length))
        else
            # create new list from index list by selecting randomly until gap count is reached
            # needs to sort at end
            push!(gap_list, create_random_less_gaps(gap_counts[i], final_length))
        end
    end

    return gap_list
end

function get_fitness_population(input_sequence, gap_population, len, score_list, score_type, gap_ext_penalty)
    gap_pop_w_fitness = []

    for i = 1:length(gap_population)
        if score_type == 2
            push!(gap_pop_w_fitness, [score_gap_list_pam250(input_sequence, gap_population[i], len, score_list, gap_ext_penalty), gap_population[i]])
        else
            push!(gap_pop_w_fitness, [score_gap_list(input_sequence, gap_population[i], len, score_list), gap_population[i]])
        end
    end

    # needs to be sorted, for easy checking of best scores
    #sort!(gap_pop_w_fitness, by = x -> (x[1]), rev = true)

    return gap_pop_w_fitness
end

function print_fitness_population(input_sequence, fitness_population)
    out = ""
    for i = 1:length(fitness_population)
        out = string(out, fitness_population[i][1], ", ")
    end
    out = string(out, "\n")

    strings_list = []

    for i = 1:length(fitness_population)
        push!(strings_list, build_strings(input_sequence, fitness_population[i][2]))
    end

    for i = 1:length(input_sequence)
        for j = 1:length(strings_list)
            out = string(out, strings_list[j][i], "   ")
        end
        out = string(out, "\n")
    end
    print(out)
    return out
end

function MSA(input_sequence, score_list, init_pop_size, gap_growth, elitism_proportion, num_crossover, children_cap, generation_count,
             mutation_strength, crossover_criteria, crossover_version, return_early, mutation_chance, printout, score_type, gap_ext_penalty)

    sequence = [
        "___________________________________MEVKKTSWTEEEDRILYQAHKRLG_NRWAEIAKLLP_____GRTDNAIKNHWNSTMRRKV____________________",
        "______________________________________SHPTYSEMIAAAIRAEKSRGG_SSRQSIQKYIKSHYKVGHNADLQIKLSIRRLLAAGVLKQTKGVGASGSFRLAK__",
        "______________________________________PRGSALSDTERAQLDVMKLLN_VSLHEMSRKIS______RSRHCIRVYLKDPVSYGTS___________________",
        "_____________________________________MRSSAKQEELVKAFKALLKEEKFSSQGEIVAALQEQGFDNINQSKVSRMLTKFGAVRTRNAKMEMVYCLPAELGVPTT",
        "SAAMAEQRHQEWLRFVDLLKNAYQNDLHLPLLNLMLTPDEREALGTRVRIIEELLRGE__MSQRELKNELG______AGIATITRGSNSLKAAPVELRQWLEEVLLKSD______"
    ]

    # get the gaps needed to increase each string to a calculated length
    gap_counts, len = get_gap_count(input_sequence, gap_growth)

    # print the gap count and length
    if printout
        # show gap counts and len
        print("Gaps needed to extend sequences to a given calculated length.", '\n')
        print(gap_counts, '\n')
        print(len, '\n')
        print('\n')
    end

    # initialize the first population, return it's chromosome representation
    @time chrom_rep_population = create_gap_population(gap_counts, len, initial_population_size)

    # print the randomized chromosome representation population
    if printout
        print("The population of random gap positions that make new sequence alignments.", '\n')
        # show the randomized population
        for i = 1:length(chrom_rep_population)
            println(chrom_rep_population[i], '\n')
        end
        println()
    end

    # get list of the chromosome representations sorted by fitness
    fitness_population = get_fitness_population(input_sequence, chrom_rep_population, len, score_list, score_type, gap_ext_penalty)

    # it is unsorted at this point
    sort!(fitness_population, by = x -> (x[1]), rev = true)

    # print the score of the populations and the reconstructed strings
    if printout
        print("The first populations score and associated strings.", '\n')
        print_fitness_population(input_sequence, fitness_population)
        print('\n')
    end

    # now we have a list with fitness values, start building generations
    # and operating on them
    while generation_count >= 1
        if printout
            print("Generation: ", generation_count)
            print("-----------------------------------------------------------", '\n')
        end
        # elitism, select a proportion of these
        fitness_population, original_population = elitism(fitness_population, elitism_prop, children_cap)

        if printout
            print("The elitism selection.", '\n')
            print_fitness_population(input_sequence, fitness_population)
            print("The population subject to crossover.", '\n')
            print_fitness_population(input_sequence, original_population)
            print('\n')
        end

        # do crossover on best, plus some of rest depending on the number of crossover
        if crossover_version == 1
            children = crossover(input_sequence, original_population, num_crossover, children_cap, len, score_list, score_type, crossover_criteria, return_early, gap_ext_penalty)
        else
            children = crossover2(input_sequence, original_population, num_crossover, children_cap, len, score_list, score_type, crossover_criteria, return_early, gap_ext_penalty)
        end

        if printout
            println("New children: ",  length(children), '\n')
        end

        # remove score for mutation, which will rescore
       children_unscored = []

        # TODO see all where not doing a deepcopy will effect the outcome
       for i = 1:length(children)
           push!(children_unscored, deepcopy(children[i][2]))
       end

        if printout
            print("Crossed over children.", '\n')
            print_fitness_population(input_sequence, children[1:min(end, 20)])
            print('\n')
        end

        # do a mutation on children, return scored
        children = mutate(input_sequence, children_unscored, mutation_strength, mutation_chance, len, score_list, score_type, gap_ext_penalty)

        if printout
            print("Mutated crossed over children.", '\n')
            print_fitness_population(input_sequence, children[1:min(end, 20)])
            print('\n')
        end

        # get list of the chromosome representations sorted by fitness
        if length(children) == 0
            append!(fitness_population, get_fitness_population(input_sequence, create_gap_population(gap_counts, len, 10), len, score_list, score_type, gap_ext_penalty))
        else
            # place the children into the next population set
            append!(fitness_population, children)
            append!(fitness_population, get_fitness_population(input_sequence, create_gap_population(gap_counts, len, 1), len, score_list, score_type, gap_ext_penalty))
        end

        # make sure to keep it sorted so elitism and crossover works again
        # became unsorted when adding children back to fitness (even though children are sorted atm)
        sort!(fitness_population, by = x -> (x[1]), rev = true)

        if printout
            print("The new population after crossover and mutation, max 20", '\n')
            print_fitness_population(input_sequence, fitness_population[1:min(end, 20)])
            print('\n')
        end

        if generation_count % 100 == 0
            println(string(generation_count, ' ', fitness_population[1][1], " New children: ",  length(children)))
            print_fitness_population(input_sequence, fitness_population[1:1])
            fitness_string = build_strings(input_sequence, fitness_population[1][2])
            sequence2 = [
                fitness_string[1],
                fitness_string[2],
                fitness_string[3],
                fitness_string[4],
                fitness_string[5]
            ]
            pairwise_score(sequence, sequence2)
        end

        generation_count -= 1
    end

    #print_fitness_population(input_sequence, fitness_population[1:min(end, 20)])
    #print_fitness_population(input_sequence, fitness_population[1:1])
    print_fitness_population(input_sequence, fitness_population[1:1])

    println(score_gap_list_pam250(input_sequence, fitness_population[1][2], 115, score_list, gap_ext_penalty))

    fitness_string = build_strings(input_sequence, fitness_population[1][2])
    #println(fitness_string)

    sequence = [
        "___________________________________MEVKKTSWTEEEDRILYQAHKRLG_NRWAEIAKLLP_____GRTDNAIKNHWNSTMRRKV____________________",
        "______________________________________SHPTYSEMIAAAIRAEKSRGG_SSRQSIQKYIKSHYKVGHNADLQIKLSIRRLLAAGVLKQTKGVGASGSFRLAK__",
        "______________________________________PRGSALSDTERAQLDVMKLLN_VSLHEMSRKIS______RSRHCIRVYLKDPVSYGTS___________________",
        "_____________________________________MRSSAKQEELVKAFKALLKEEKFSSQGEIVAALQEQGFDNINQSKVSRMLTKFGAVRTRNAKMEMVYCLPAELGVPTT",
        "SAAMAEQRHQEWLRFVDLLKNAYQNDLHLPLLNLMLTPDEREALGTRVRIIEELLRGE__MSQRELKNELG______AGIATITRGSNSLKAAPVELRQWLEEVLLKSD______"
    ]

    println()
    println(string(fitness_string[1], '\n', fitness_string[2], '\n', fitness_string[3], '\n', fitness_string[4], '\n', fitness_string[5]))
    println()

    sequence2 = [
        fitness_string[1],
        fitness_string[2],
        fitness_string[3],
        fitness_string[4],
        fitness_string[5]
    ]

    println(string("balibase score: ", score_sequence_pam250(sequence, 115, score_list, gap_ext_penalty)))
    println(string("our score: ", score_sequence_pam250(sequence2, 115, score_list, gap_ext_penalty)))

    pairwise_score(sequence, sequence2)
    return fitness_population[1]
end

export MSA


text:
# called directly by alg_cache for iip algorithms and indirectly by non-sc algorithms
function SaddleSystem(A,f::ConstrainedODEFunction{true},p,pold,ducache,solver)
    nully, nullz = state(ducache), constraint(ducache)
    du_aux = aux_state(ducache)
    @inline B₁ᵀ(z) = (zero_vec!(ducache);
                     _ode_neg_B1!(ducache,f,solvector(state=nully,constraint=z,aux_state=du_aux),pold,0.0);
                     state(ducache) .*= -1.0; return state(ducache))
    @inline B₂(y) = (zero_vec!(ducache);
                     _constraint_neg_B2!(ducache,f,solvector(state=y,constraint=nullz,aux_state=du_aux),p,0.0);
                     constraint(ducache) .*= -1.0; return constraint(ducache))
    SaddleSystem(A,B₂,B₁ᵀ,mainvector(ducache),solver=solver)
end

# called directly by oop algorithms
function SaddleSystem(A,f::ConstrainedODEFunction{false},p,pold,ducache,solver)

    nully, nullz = state(ducache), constraint(ducache)
    du_aux = aux_state(ducache)
    @inline B₁ᵀ(z) = (zero_vec!(ducache); ducache .= _ode_neg_B1(f,solvector(state=nully,constraint=z,aux_state=du_aux),pold,0.0);
                     state(ducache) .*= -1.0; return state(ducache))
    @inline B₂(y) = (zero_vec!(ducache);
                     ducache .= _constraint_neg_B2(f,solvector(state=y,constraint=nullz,aux_state=du_aux),p,0.0);
                     constraint(ducache) .*= -1.0; return constraint(ducache))
    SaddleSystem(A,B₂,B₁ᵀ,mainvector(ducache),solver=solver)
end

# this version is called by in-place algorithms
@inline SaddleSystem(S::SaddleSystem,A,f::ConstrainedODEFunction,p,pold,
                      cache::ConstrainedODEMutableCache{sc,solverType}) where {sc,solverType} =
          SaddleSystem(S,A,f,p,pold,cache.dutmp,solverType,Val(sc))

# non-static constraints
@inline SaddleSystem(S::SaddleSystem,A,f::ConstrainedODEFunction,p,pold,ducache,solver,
                      ::Val{false}) = SaddleSystem(A,f,p,pold,ducache,solver)

# static constraints
@inline SaddleSystem(S::SaddleSystem,A,f::ConstrainedODEFunction,p,pold,ducache,solver,
                      ::Val{true}) = S

function B1_times_z!(u,S::SaddleSystem)
    state(u) .= typeof(state(u))(S.A⁻¹B₁ᵀf)
    return u
end

function B1_times_z(u,S::SaddleSystem)
    out = zero(u)
    state(out) .= typeof(state(u))(S.A⁻¹B₁ᵀf)
    return out
end


text:
mutable struct ObjectiveContainer{MOIForm,ArrayForm} <: MOI.ModelLike
    num_variables::Int64
    is_sense_set::Bool
    sense::MOI.OptimizationSense
    is_function_set::Bool
    func::Union{Nothing,ArrayForm}
    function ObjectiveContainer{MOIForm,ArrayForm}() where {MOIForm,ArrayForm}
        o = new{MOIForm,ArrayForm}()
        MOI.empty!(o)
        return o
    end
end

function MOI.empty!(o::ObjectiveContainer)
    o.num_variables = 0
    o.is_sense_set = false
    o.sense = MOI.FEASIBILITY_SENSE
    o.is_function_set = false
    o.func = nothing
    return
end

function MOI.is_empty(o::ObjectiveContainer)
    return !o.is_sense_set && !o.is_function_set
end

# FIXME type piracy
function MOI.add_variable(model::MOI.Utilities.AbstractModel)
    x = MOI.add_variable(model.variables)
    MOI.Utilities._add_variable(model.objective)
    MOI.Utilities._add_variable(model.constraints)
    return x
end
function MOI.Utilities._add_variable(::MOI.Utilities.ObjectiveContainer) end
function MOI.Utilities._add_variable(model::ObjectiveContainer)
    if !model.is_function_set
        error("Adding variables after setting objective is not implemented yet")
    end
    model.num_variables += 1
    return
end

###
### ObjectiveSense
###

MOI.supports(::ObjectiveContainer, ::MOI.ObjectiveSense) = true

MOI.get(o::ObjectiveContainer, ::MOI.ObjectiveSense) = o.sense

function MOI.set(o::ObjectiveContainer, ::MOI.ObjectiveSense, value)
    if value == MOI.FEASIBILITY_SENSE
        MOI.empty!(o)
    end
    o.is_sense_set = true
    o.sense = value
    return
end

###
### ObjectiveFunctionType
###

function MOI.get(
    ::ObjectiveContainer{MOIForm},
    ::MOI.ObjectiveFunctionType,
) where {MOIForm}
    return MOIForm
end

###
### ObjectiveFunction
###

function MOI.supports(
    ::ObjectiveContainer{MOIForm},
    ::MOI.ObjectiveFunction{MOIForm},
) where {MOIForm}
    return true
end

function MOI.get(
    o::ObjectiveContainer{MOIForm},
    ::MOI.ObjectiveFunction{MOIForm},
) where {MOIForm}
    return convert(MOIForm, o.func)
end

function MOI.set(
    o::ObjectiveContainer,
    ::MOI.ObjectiveFunction{MOIForm},
    f::MOIForm,
) where {MOIForm}
    o.is_function_set = true
    o.func = sparse_array_representation(f, o.num_variables)
    return
end

###
### MOI.ListOfModelAttributesSet
###

function MOI.get(o::ObjectiveContainer, ::MOI.ListOfModelAttributesSet)
    ret = MOI.AbstractModelAttribute[]
    if o.is_sense_set
        push!(ret, MOI.ObjectiveSense())
    end
    if o.is_function_set
        F = MOI.get(o, MOI.ObjectiveFunctionType())
        push!(ret, MOI.ObjectiveFunction{MOIForm}())
    end
    return ret
end


text:
using .Turing: DynamicPPL, Prior

function (alg::ThermInt)(model::DynamicPPL.Model; progress=true, kwargs...)
    nsteps = length(alg.schedule)
    p = Progress(nsteps; enabled=progress, desc="TI sampling")
    ΔlogZ = [
        begin
            ProgressMeter.next!(p)
            evaluate_loglikelihood(model, alg, β)
        end for β in alg.schedule
    ]
    return trapz(alg.schedule, ΔlogZ)
end

function (alg::ThermInt)(
    model::DynamicPPL.Model, ::TIParallelThreads; progress=true, kwargs...
)
    nsteps = length(alg.schedule)
    nthreads = min(Threads.nthreads(), nsteps)
    ΔlogZ = zeros(Float64, nsteps)
    algs = [deepcopy(alg) for _ in 1:nthreads]
    p = Progress(nsteps; enabled=progress, desc="TI sampling")
    Threads.@threads for i in 1:nsteps
        id = Threads.threadid()
        ΔlogZ[i] = evaluate_loglikelihood(model, algs[id], alg.schedule[i]; kwargs...)
        ProgressMeter.next!(p)
    end
    return trapz(alg.schedule, ΔlogZ)
end

function evaluate_loglikelihood(model::DynamicPPL.Model, alg::ThermInt, β::Real)
    powerlogπ = power_logjoint(model, β)
    loglikelihood = get_loglikelihood(model)
    x_init = vec(Array(sample(model, Prior(), 1))) # Bad ugly hack cause I don't know how to sample from the prior
    samples = sample_powerlogπ(powerlogπ, alg, x_init)
    return mean(loglikelihood, samples)
end

function power_logjoint(model, β)
    ctx = DynamicPPL.MiniBatchContext(DynamicPPL.DefaultContext(), β)
    spl = DynamicPPL.SampleFromPrior()
    vi = DynamicPPL.VarInfo(model)
    return function f(z)
        varinfo = DynamicPPL.VarInfo(vi, spl, z)
        model(varinfo, spl, ctx)
        return DynamicPPL.getlogp(varinfo)
    end
end

function get_loglikelihood(model)
    ctx = DynamicPPL.LikelihoodContext()
    spl = DynamicPPL.SampleFromPrior()
    vi = DynamicPPL.VarInfo(model)
    return function f(z)
        varinfo = DynamicPPL.VarInfo(vi, spl, z)
        model(varinfo, spl, ctx)
        return DynamicPPL.getlogp(varinfo)
    end
end