text:
@inline function get_geom_weights_outerloop(
    geom_weights::NTuple{3, <:AbstractArray{T, 2}}, i, ::Type{Cylindrical}
) where {T}
    geom_weights[1][1, i],
    geom_weights[1][2, i],
    geom_weights[1][3, i],
    geom_weights[1][4, i],
    geom_weights[1][5, i],
    geom_weights[1][6, i]
end

@inline function get_geom_weights_outerloop(
    geom_weights::NTuple{3, <:AbstractArray{T, 2}}, i, ::Type{Cartesian}
) where {T}
    geom_weights[1][1, i],
    geom_weights[1][2, i],
    geom_weights[1][3, i],
    geom_weights[1][4, i + 1],
    geom_weights[1][4, i]
end

@inline function prepare_weights_in_middleloop(
    geom_weights::NTuple{3, <:AbstractArray{T, 2}}, ::Type{Cylindrical},
    i2, in2,
    pwwrr, pwwrl, r_inv_pwΔmpr, Δr_ext_inv_r_pwmprr, Δr_ext_inv_l_pwmprl, Δmpr_squared,
    is_r0_t::Val{is_r0}
) where {T, is_r0}
    pwwφr        = geom_weights[2][1, in2]
    pwwφl        = geom_weights[2][2, in2]
    pwΔmpφ       = geom_weights[2][3, in2]
    Δφ_ext_inv_r = geom_weights[2][4,  i2]
    Δφ_ext_inv_l = geom_weights[2][4, in2]

    if is_r0
        pwwφr = T(0.5)
        pwwφl = T(0.5)
        pwΔmpφ = T(2π)
        Δφ_ext_inv_r = inv(pwΔmpφ)
        Δφ_ext_inv_l = Δφ_ext_inv_r
    end
    pwwrr_pwwφr = pwwrr * pwwφr
    pwwrl_pwwφr = pwwrl * pwwφr
    pwwrr_pwwφl = pwwrr * pwwφl
    pwwrl_pwwφl = pwwrl * pwwφl

    pwΔmpφ_Δmpr_squared = pwΔmpφ * Δmpr_squared
    Δr_ext_inv_r_pwmprr_pwΔmpφ = Δr_ext_inv_r_pwmprr * pwΔmpφ
    Δr_ext_inv_l_pwmprl_pwΔmpφ = Δr_ext_inv_l_pwmprl * pwΔmpφ
    r_inv_pwΔmpr_Δφ_ext_inv_r = r_inv_pwΔmpr * Δφ_ext_inv_r
    r_inv_pwΔmpr_Δφ_ext_inv_l = r_inv_pwΔmpr * Δφ_ext_inv_l
    return (
        pwwrr, pwwrl, pwwφr, pwwφl,
        pwwrr_pwwφr, pwwrl_pwwφr, pwwrr_pwwφl, pwwrl_pwwφl,
        pwΔmpφ_Δmpr_squared,
        Δr_ext_inv_r_pwmprr_pwΔmpφ, Δr_ext_inv_l_pwmprl_pwΔmpφ,
        r_inv_pwΔmpr_Δφ_ext_inv_r, r_inv_pwΔmpr_Δφ_ext_inv_l
    )
end

@inline function prepare_weights_in_middleloop(
    geom_weights::NTuple{3, <:AbstractArray{T, 2}}, ::Type{Cartesian},
    i2, in2,
    pww3r, pww3l, pwΔmp3, Δ3_ext_inv_r, Δ3_ext_inv_l,
    is_r0_t
) where {T}
    pww2r  = geom_weights[2][1, in2]
    pww2l  = geom_weights[2][2, in2]
    pwΔmp2 = geom_weights[2][3, in2]
    pwΔmp2_pwΔmp3 = pwΔmp2 * pwΔmp3
    Δ2_ext_inv_r_pwΔmp3  = geom_weights[2][4, in2 + 1] * pwΔmp3
    Δ2_ext_inv_l_pwΔmp3  = geom_weights[2][4, in2]     * pwΔmp3
    Δ3_ext_inv_r_pwΔmp2 = Δ3_ext_inv_r * pwΔmp2
    Δ3_ext_inv_l_pwΔmp2 = Δ3_ext_inv_l * pwΔmp2

    pww2r_pww3r = pww2r * pww3r
    pww2r_pww3l = pww2r * pww3l
    pww2l_pww3r = pww2l * pww3r
    pww2l_pww3l = pww2l * pww3l
    return (
        pww3r, pww3l, pww2r, pww2l,
        pww2r_pww3r, pww2r_pww3l, pww2l_pww3r, pww2l_pww3l,
        pwΔmp2_pwΔmp3,
        Δ3_ext_inv_r_pwΔmp2, Δ3_ext_inv_l_pwΔmp2,
        Δ2_ext_inv_r_pwΔmp3, Δ2_ext_inv_l_pwΔmp3,
    )
end

@fastmath function middleloop!(
    i3::Int, rb_tar_idx::Int, rb_src_idx::Int,
    pcs::PotentialCalculationSetup{T, S},
    update_even_points::Val{even_points},
    depletion_handling::Val{depletion_handling_enabled},
    is_weighting_potential::Val{_is_weighting_potential},
    only2d::Val{only_2d},
    idx3iseven::Val{idx3_is_even},
    is_r0_t::Val{is_r0}
)::Nothing where {T, S, even_points, depletion_handling_enabled, _is_weighting_potential, only_2d, idx3_is_even, is_r0}
    @inbounds begin
        in3 = i3 - 1

        geom_weights_3 = get_geom_weights_outerloop(pcs.geom_weights, in3, S)

        line_weights::Array{T, 2} = Array{T, 2}(undef, size(pcs.potential, 1) - 2, 6)
        # Even though this causes some allocations it
        # is faster than using a predefined array, e.g. stored in pcs
        # Especially when using multiple threads

        #=
            Splitting this into two loops over even and uneven i2
            seems to increase the performance quite a bit by hard-coding oscillating type
            instead of one loop using:
                rφi_is_even_t::Union{Val{true}, Val{false}} = Val(iseven(i3 + i2) ? true : false)
        =#
        for i2 in (only_2d ? (2,) : 2:2:(size(pcs.potential, 2) - 1))
            in2 = i2 - 1
            i23_is_even_t = Val(idx3_is_even ? true : false)

            geom_weights_2 = prepare_weights_in_middleloop(
                pcs.geom_weights, S, i2, in2,
                geom_weights_3...,
                is_r0_t
            )

            calculate_weights_for_innerloop!(line_weights, pcs, i2, in2, i3, in3,
                update_even_points, i23_is_even_t,
                geom_weights_2...
            )

            innerloop!(line_weights, pcs, i2, in2, i3, in3, rb_tar_idx, rb_src_idx,
                update_even_points, i23_is_even_t,
                depletion_handling, is_weighting_potential, only2d)
        end
        for i2 in 3:2:(size(pcs.potential, 2) - 1)
            in2 = i2 - 1
            i23_is_even_t = Val(idx3_is_even ? false : true)

            geom_weights_2 = prepare_weights_in_middleloop(
                pcs.geom_weights, S, i2, in2,
                geom_weights_3...,
                is_r0_t
            )

            calculate_weights_for_innerloop!(line_weights, pcs, i2, in2, i3, in3,
                update_even_points, i23_is_even_t,
                geom_weights_2...
            )

            innerloop!(line_weights, pcs, i2, in2, i3, in3, rb_tar_idx, rb_src_idx,
                update_even_points, i23_is_even_t,
                depletion_handling, is_weighting_potential, only2d)
        end
    end
end


text:
# Initial replacements to perform
replacements = Dict("á"=>"A","ch"=>"X","ç"=>"S","é"=>"E","í"=>"I","ó"=>"O","ú"=>"U","ñ"=>"NY","gü"=>"W","ü"=>"U","b"=>"V","ll"=>"Y")

# Vowels
vowels = Set(['A', 'E', 'I', 'O', 'U'])

# Consonants that make a single sound
single_sound_consonants = Set(['D','F','J','K','M','N','P','T','V','L','Y'])

function is_vowel(c::Char)
    in(c, vowels)
end

# Replace some easily confused spanish characters
function fixstr(st::AbstractString)
    for replacement in replacements
        st = replace(st, replacement)
    end
    st
end

"""
    metaphone(s::AbstractString, key_length=6)

Calculate the metaphone of a Spanish string `s` upto a given length `key_length`. Default `key_length` is 6.

# Examples
```jldoctest
julia> metaphone("consideración")
"KNSDRZ"

julia> metaphone("consideración", 7)
"KNSDRZN"

julia> metaphone("consideración", 8)
"KNSDRZN"

julia> metaphone("consideración", 2)
"KN"
```
"""
function metaphone(s::AbstractString, key_length=6)
    # initialize metaphone key string
    meta_key = ""
    # String length
    n = length(s)
    # Pointer to current character
    current_pos = 1
    #Set to the end of the string
    original_string = s* "    "
    # Replace some easily confused spanish characters
    original_string = uppercase(fixstr(lowercase(original_string)))

    # main loop
    while length(meta_key)<key_length
        #break out of the loop if greater or equal than the length
        if current_pos > length(original_string)
            break
        end
        #get character from the string
        current_char = original_string[current_pos]

        if is_vowel(current_char) && current_pos==1
            meta_key *= current_char
            current_pos = nextind(original_string, current_pos)
        else
            #Let's check for consonants  that have a single sound
            #or already have been replaced  because they share the same
            #sound like 'B' for 'V' and 'S' for 'Z'
            if in(current_char, single_sound_consonants)
                meta_key *= current_char
                current_pos = nextind(original_string, current_pos)

                #increment by two if a repeated letter is found
                if original_string[current_pos] == current_char
                    current_pos = nextind(original_string, current_pos)
                end
            else #check consonants with similar confusing sounds
                if current_char == 'C'
                    #special case 'macho', chato,etc.
                    #if original_string[nextind(original_string, current_pos)]=='H'
                    #    current_pos = nextind(original_string, current_pos, 2)

                    #special case 'acción', 'reacción',etc.
                    if original_string[nextind(original_string, current_pos)] == 'C'
                        meta_key*= 'X'
                        current_pos = nextind(original_string, current_pos, 2)
                    # special case 'cesar', 'cien', 'cid', 'conciencia'
                    elseif in(original_string[nextind(original_string, current_pos)], ['E','I'])
                        meta_key *= 'Z'
                        current_pos = nextind(original_string, current_pos, 2)
                    else
                        meta_key *= 'K'
                        current_pos = nextind(original_string, current_pos)
                    end
                elseif current_char == 'G'
                    # special case 'gente', 'ecologia',etc
                    if in(original_string[nextind(original_string, current_pos)], ['E','I'])
                        meta_key *= 'J'
                        current_pos = nextind(original_string, current_pos, 2)
                    else
                        meta_key *= 'G'
                        current_pos = nextind(original_string, current_pos)
                    end
                #since the letter 'h' is silent in spanish,
                #let's set the meta key to the vowel after the letter 'h'
                elseif current_char == 'H'
                    if is_vowel(original_string[nextind(original_string, current_pos)])
                        meta_key *= original_string[nextind(original_string, current_pos)]
                        current_pos = nextind(original_string, current_pos, 2)
                    else
                        meta_key *= 'H'
                        current_pos = nextind(original_string, current_pos)
                    end
                elseif current_char == 'Q'
                    if original_string[nextind(original_string, current_pos)]=='U'
                        current_pos = nextind(original_string, current_pos, 2)
                    else
                        meta_key *= 'K'
                        current_pos = nextind(original_string, current_pos)
                    end
                elseif current_char == 'W'
                    meta_key *= 'U'
                    current_pos = nextind(original_string, current_pos)

                # perro, arrebato, cara
                elseif current_char == 'R'
                    meta_key *= 'R'
                    current_pos = nextind(original_string, current_pos)
                # spain
                elseif current_char == 'S'
                    if !is_vowel(original_string[nextind(original_string, current_pos)]) && current_pos == 1
                        meta_key *= "ES"
                        current_pos = nextind(original_string, current_pos)
                    else
                        meta_key *= 'S'
                        current_pos = nextind(original_string, current_pos)
                    end
                # zapato
                elseif current_char == 'Z'
                    meta_key *= 'Z'
                    current_pos = nextind(original_string, current_pos)
                elseif current_char == 'X'
                    if !is_vowel(original_string[nextind(original_string, current_pos)]) && n>1 && current_pos == 1
                        meta_key *= "EX"
                        current_pos = nextind(original_string, current_pos)
                    else
                        meta_key *= 'X'
                        current_pos = nextind(original_string, current_pos)
                    end
                else
                    current_pos = nextind(original_string, current_pos)
                end
            end
        end
    end
    strip(meta_key)
end

ext: jl


text:
using Distributions

type PredictionResult
    proteinid::UTF8String
    decisionValues::Array{Float64}
    predictionLabels::Array{Int64}
end

type Graph
    proteinid::UTF8String
    graph::Dict{Int64, Array{Int64}}
end

function parsePredictionFile(predictionFile)
    predictions = Dict{UTF8String, PredictionResult}()
    decisionValues = Float64[]
    predictionLabels = Int64[]
    prevProteinid = ""
    c = 0
    for line in eachline(open(predictionFile))
        if length(rstrip(line)) == 0
            continue
        end
        proteinid, index, decisionValue, predictionLabel = split(line, '\t')
        if c != 0 && proteinid != prevProteinid
            predictionLabels = map(x-> if x == 0 return -1 else return 1 end, predictionLabels)
            predictionResult = PredictionResult(prevProteinid, decisionValues, predictionLabels)
            predictions[prevProteinid] = predictionResult
            dicisionValues = Float64[]
            predictionLabels = Int64[]
            push!(decisionValues, float(decisionValue))
            push!(predictionLabels, parse(Int, predictionLabel))
        else
            push!(decisionValues, float(decisionValue))
            push!(predictionLabels, parse(Int, predictionLabel))
        end
        prevProteinid = proteinid
        c += 1
    end
    if length(decisionValues) != 0 && length(predictionLabels) != 0
        predictionLabels = map(x-> if x == 0 return -1 else return 1 end, predictionLabels)
        predictionResult = PredictionResult(prevProteinid, decisionValues, predictionLabels)
        predictions[prevProteinid] = predictionResult
    end
    return predictions
end

function parseCorrectFile(correctFile)
    corrects = Dict{UTF8String, Array{Int64, 1}}()
    correctLabels = Int64[]
    prevProteinid = ""
    c = 0
    for line in eachline(open(correctFile))
        if length(rstrip(line)) == 0
            continue
        end
        proteinid, index, correctLabel = split(line, '\t')
        if c != 0 && proteinid != prevProteinid
            corrects[prevProteinid] = correctLabels
            correctLabels = Int64[]
        end
        push!(correctLabels, parse(Int, correctLabel))
        prevProteinid = proteinid
        c += 1
    end
    if length(correctLabels) != 0
        corrects[prevProteinid] = correctLabels
    end
    return corrects
end

function parseGraphFile(proteinid, graphFile)
    graph = Dict{Int64, Array{Int64}}()
    for line in eachline(open(graphFile))
        if length(rstrip(line)) == 0
            continue
        end
        left, right = map(x->parse(Int, x), split(line, ' '))

        #convert to Julia index (start from 1)
        left += 1
        right += 1
        if haskey(graph, left)
            push!(graph[left], right)
        else
            graph[left] = [right]
        end
    end
    return Graph(proteinid, graph)
end

function parseGraphFiles(predictions, graphDir)
    graphs = Dict{UTF8String, Graph}()
    for proteinid in keys(predictions)
        graph = parseGraphFile(proteinid, graphDir*proteinid*".graph")
        graphs[proteinid] = graph
    end
    return graphs
end

function getAllState(numNeighbor)
    allState = ASCIIString[]
    for i in 0:2^numNeighbor-1
        push!(allState, bits(i)[64-numNeighbor+1:end])
    end
    return allState
end

function expEnergy(y, x, σ, J, d)
    return exp(-(y-x)^2/(2*(σ^2)))*exp(-2*J*(d))
end

function disagreeingEdgeNumber(selfState, xs, neighbors)
    d = 0
    for neighbor in neighbors
        if xs[neighbor] != selfState
            d += 1
        end
    end
    return d
end

function encodeState(state)
    encodedStates =  map(x->if x == "0" return -1 elseif x == "1" return 1 end, split(state, ""))
    return encodedStates[1], encodedStates[2:end]
end

function isDisagreeing(selfState, neighborState)
    if selfState != neighborState
        return 1
    end
    return 0
end

function pseudolikelihood(xs, ys, σ, J, graph::Dict{Int64, Array{Int64}})
    logPseLike = 0
    for i in 1:length(xs)

        # For x
        if haskey(graph, i)
            neighbors = graph[i]
            #Z = 0
            s = xs[i] # current state
            s′ = -1*s # flipped state
            d  = disagreeingEdgeNumber(s, xs, neighbors)
            d′ = disagreeingEdgeNumber(s′, xs, neighbors)
            r  = exp((2*ys[i]*s)/(σ^2))*exp(-2*J*(d-d′))  # currentProb / flippedProb
            p = r / (1 + r)
            #d = disagreeingEdgeNumber(-xs[i], xs, neighbors)
            #Z += expEnergy(ys[i], -xs[i], σ, J, d)
            #d = disagreeingEdgeNumber(xs[i], xs, neighbors)
            #currExpEnergy = expEnergy(ys[i], xs[i], σ, J, d)
            #Z += currExpEnergy
            #p = currExpEnergy / Z
            logPseLike += log(p)
        else # has No Edge.
            #Z = 0
            s = xs[i] # current state
            s′ = -1*s     # flipped state
            r = exp((2*ys[i]*s)/(σ^2))  # currentProb / flippedProb
            p = r / (1 + r)
            #Z += expEnergy(ys[i], xs[i], σ, J, 0)
            #Z += expEnergy(ys[i], -xs[i], σ, J, 0)
            #currExpEnergy = expEnergy(ys[i], xs[i], σ, J, 0)
            #p = currExpEnergy / Z
            logPseLike += log(p)
        end
    end
    return logPseLike
end

function gibbsSampling(ys, iternum, σ, J, graph::Dict{Int64, Array{Int64}}, burnIn)
    maxIndex = length(ys)
    #dist = Binomial(1, 0.5)
    #currentStates = map(x->if x == 0 return -1 else return 1 end, rand(dist, maxIndex))
    currentStates = copy(ys)
    sampled = Array{Int8, 1}[]
    for i in 1:maxIndex
        push!(sampled, [])
    end
    for i in 1:iternum
        index = rand(1:maxIndex)
        if haskey(graph, index)
            neighbors = graph[index]
            s = currentStates[index] # current state
            s′ = -1*s     # flipped state
            d  = disagreeingEdgeNumber(s, currentStates, neighbors)
            d′ = disagreeingEdgeNumber(s′, currentStates, neighbors)
            r′ = exp((2*ys[index]*s′)/(σ^2))*exp(-2*J*(d′-d))  # flippedProb / currentProb
            p′ = r′ / (1 + r′)
        else # has No Edge.
            s = currentStates[index] # current state
            s′ = -1*s     # flipped state
            r′ = exp((2*ys[index]*s′)/(σ^2))  # flippedProb / currentProb
            p′ = r′ / (1 + r′)
        end
        dist = Binomial(1, p′)
        flip_or_remain = rand(dist) # 1 means flip, 0 means remain.
        if flip_or_remain == 1
            newState = s′
        else
            newState = s
        end
        if i > burnIn
            push!(sampled[index], newState)
        end
        currentStates[index] = newState
    end
    return map(mean, sampled)
end

function main()
    try
        predictionFile = ARGS[1]
        correctFile    = ARGS[2]
        graphDir       = ARGS[3]
    catch exception
        println(Usage)
        exit()
    end
    predictionFile = ARGS[1]
    correctFile    = ARGS[2]
    graphDir       = ARGS[3]
    iternum        = 100000
    burnIn         = 10000
    σRange         = (0.1, 2.0)
    JRange         = (0.1, 2.0)
    σGrid          = 100
    JGrid          = 100
    predictions =  parsePredictionFile(predictionFile)
    corrects = parseCorrectFile(correctFile)
    graphs = parseGraphFiles(predictions, graphDir)
    logPseudolikelihoods = Tuple[]
    maxPseudolikelihood = (0.0, 0.0, -Inf) # σ, J, logPseudolikelihood
    for σG in 0:σGrid
        σ = (σG/float(σGrid)) * (σRange[2]-σRange[1]) + σRange[1]
        for JG in 0:JGrid
            J = (JG/float(JGrid)) * (JRange[2]-JRange[1]) + JRange[1]
            logPseLike = 0
            for proteinid in keys(predictions)
                logPseLike += pseudolikelihood(corrects[proteinid], predictions[proteinid].predictionLabels,
                                                                        σ, J, graphs[proteinid].graph)
            end
            push!(logPseudolikelihoods, (σ, J, logPseLike))
            if logPseLike > maxPseudolikelihood[3]
                maxPseudolikelihood = (σ, J, logPseLike)
            end
        end
    end
    sort!(logPseudolikelihoods, by=x->x[3], rev=true)
    for logPseudolikelihood in logPseudolikelihoods[1:10]
        @show logPseudolikelihood
    end
    @show maxPseudolikelihood
    proteinApproxMeans = Dict{UTF8String, Array{Array{Float64, 1}}}()
    for proteinid in keys(predictions)
        proteinApproxMeans[proteinid] = Array{Array{Float64, 1}}[]
    end
    for i in 1:10
        for proteinid in keys(predictions)
            approxMean =  gibbsSampling(predictions[proteinid].predictionLabels, iternum, maxPseudolikelihood[1],
                                                    maxPseudolikelihood[2], graphs[proteinid].graph, burnIn)
            push!(proteinApproxMeans[proteinid], approxMean)
        end
    end
    for proteinid in keys(predictions)
        @show proteinid
        for eachMean in proteinApproxMeans[proteinid]
            @show eachMean
        end
        @show mean(proteinApproxMeans[proteinid])
    end
end

Usage = """

Usage:
    julia gibbs_sampler.jl <prediction_output> <correct_label_file> <graph_dir>

FORMAT:

<prediction_output>
    proteinid\tindex\t0
    proteinid\tindex\t1
    proteinid\tindex\t0
    proteinid\tindex\t0
    ...

<correct_label_file>
    proteinid\tindex\t0
    proteinid\tindex\t1
    proteinid\tindex\t0
    proteinid\tindex\t0
    ...

<graph_dir> is a directory that includes protein_graph_file
protein_graph_file must be "proteinid.graph"

protein_graph_file FORMAT:

    0 1
    1 4
    2 3
    2 4
    3 1
    ...
"""

main()

max_stars_repo_name: clclcocoro/GibbsSampling

max_issues_repo_name: clclcocoro/GibbsSampling

max_forks_repo_name: clclcocoro/GibbsSampling


text:
const START_TIME = time()

@show LoopVectorization.register_count()

@show RUN_SLOW_TESTS

@time @testset "LoopVectorization.jl" begin

  @time if LOOPVECTORIZATION_TEST == "all" || LOOPVECTORIZATION_TEST == "part1"
    @time include("broadcast.jl")
  end

  @time if LOOPVECTORIZATION_TEST == "all" || LOOPVECTORIZATION_TEST == "part2"
    using Aqua
    @time Aqua.test_all(LoopVectorization, ambiguities = VERSION ≥ v"1.6", stale_deps = false)
    @time Aqua.test_stale_deps(LoopVectorization, ignore = [:SIMDDualNumbers])
    # @test isempty(detect_unbound_args(LoopVectorization))

    @time include("printmethods.jl")

    @time include("can_avx.jl")

    @time include("fallback.jl")

    @time include("utils.jl")

    @time include("arraywrappers.jl")

    @time include("check_empty.jl")

    @time include("loopinductvars.jl")

    @time include("shuffleloadstores.jl")

    if VERSION < v"1.7-DEV"
      @time include("zygote.jl")
    else
      println("Skipping Zygote tests.")
    end

    @time include("map.jl")

    @time include("filter.jl")

    @time include("multiassignments.jl")

    @time include("reduction_untangling.jl")

    @time include("manyloopreductions.jl")

    @time include("simplemisc.jl")

    @time include("convolutions.jl")

    @time include("ifelsemasks.jl")

    @time include("gemv.jl")

    @time include("dot.jl")

    @time include("special.jl")

    @time include("mapreduce.jl")

    @time include("index_processing.jl")
  end

  @time if LOOPVECTORIZATION_TEST == "all" || LOOPVECTORIZATION_TEST == "part3"
    @time include("rejectunroll.jl")

    @time include("miscellaneous.jl")

    @time include("copy.jl")

    @time include("tensors.jl")

    @time include("staticsize.jl")
end

  @time if LOOPVECTORIZATION_TEST == "all" || LOOPVECTORIZATION_TEST == "part4"
    @time include("threading.jl")

    @time include("tullio.jl")

    @time include("iteration_bound_tests.jl")

    @time include("outer_reductions.jl")

    @time include("upperboundedintegers.jl")

    if VERSION ≥ v"1.6"
      @time include("quantum.jl")
    end

    @time include("offsetarrays.jl")
  end

  @time if LOOPVECTORIZATION_TEST == "all" || LOOPVECTORIZATION_TEST == "part5"
    @time include("steprange.jl")

    @time include("gemm.jl")

    @time include("inner_reductions.jl")
  end

end

const ELAPSED_MINUTES = (time() - START_TIME)/60
# @test ELAPSED_MINUTES < 180
@test ELAPSED_MINUTES < 300


text:
using Documenter, Simulation73
DocMeta.setdocmeta!(Simulation73, :DocTestSetup, :(using Simulation73); recursive=true)
makedocs(
    sitename="Simulation73",
    modules=[Simulation73],
    pages = [
       "Home" => "index.md",
       "Manual" => Any[
           "man/simulation.md"
       ],
       "Library" => Any[
           "Public" => "lib/public.md",
           hide("Internals" => "lib/internals.md",
            Any[joinpath("lib/internals", f) for f in readdir("docs/src/lib/internals")]
           )
       ]
   ]
)

deploydocs(
    repo = "github.com/grahamas/Simulation73.git"
)


text:
using FractionalCalculus
using SymbolicUtils
using Test

macro eqtest(expr)
    @assert expr.head == :call && expr.args[1] in [:(==), :(!=)]
    if expr.args[1] == :(==)
        :(@test isequal($(expr.args[2]), $(expr.args[3])))
    else
        :(@test !isequal($(expr.args[2]), $(expr.args[3])))
    end |> esc
end

@testset "Test DIFF CONSTANT AND POWERS" begin
    @syms x
    @eqtest semidiff(x^2) == gamma(3)/gamma(2.5)*(x^1.5)
    @eqtest semidiff(log(x)) == log(4*x) / sqrt(π*x)
    #@eqtest semidiff(sqrt(1+x)) == 1/sqrt(pi*x)+atan(sqrt(x))/sqrt(pi)
    #@eqtest semidiff(sqrt(1-x)) == 1/sqrt(pi*x)-atanh(sqrt(x))/sqrt(pi)
    @eqtest semidiff(1/sqrt(1+x)) == 1/(sqrt(pi*x)*(1+x))
    #@eqtest semidiff(1/sqrt(1-x)) == 1/(sqrt(pi*x)*(1-x))
    @eqtest semidiff(1/(1+x)) == (sqrt(1+x)-sqrt(x)*asinh(sqrt(x)))/(sqrt(pi*x)*(1+x)^(3/2))
    #@eqtest semidiff(1/(1-x)) == (sqrt(1-x)+sqrt(x)*asinh(sqrt(x)))/(sqrt(pi*x)*(1-x)^(3/2))
    #@eqtest semidiff(1/(1+x)^(3/2)) == (1-x)/(sqrt(pi*x)*(1+x)^2)
    #@eqtest semidiff(1/(1-x)^(3/2)) == (1+x)/(sqrt(pi*x)*(1-x)^2)
end


text:
import Base: show, product, @_inline_meta
import Base: eltype, getindex
import Base: +,-,*,/,sum

struct BernsteinBase{D, N, T}
    control_points::SVector{N, SVector{D, T}}
end

function BernsteinBase(cp::Vararg{SVector,N}) where {N}
    cpts = SVector{N}(cp...)
    BernsteinBase(cpts)
end

function BernsteinBase(data::Array{Array{T, 1}}) where {T}
    cptsz = length(data)
    dimsz = length(data[1])
    BernsteinBase(data, Val(dimsz), Val(cptsz))
end

function BernsteinBase(data::Array{Array{T, 1}}, ::Val{D}, ::Val{N}) where {D, N, T}
    cpts = SVector{N}([SVector{D}(d) for d in data])
    BernsteinBase(cpts)
end

Base.eltype(::Type{BernsteinBase{D, N, T}}) where {D, N, T} = T
Base.getindex(b::BernsteinBase, j::Int) = b.control_points[j]

function _eval(b::BernsteinBase{D, N}, t::Real) where {D, N}
    h = 1.0
    u = 1.0 - t
    q = b[1]
    if t < 0.5
        u = t/u
        for k = 1:N-1
            h = h*u*(N-k)
            h = h/(k + h)
            q = (1.0-h)*q + h*b[k+1]
        end
    else
        u = u/t
        for k = 1:N-1
            h = h*(N-k)
            h = h/(k*u + h)
            q = (1.0-h)*q + h*b[k+1]
        end
    end
    return q
end

function (b::BernsteinBase)(t::Real)
    _eval(b, t)
end

function (b::BernsteinBase{1})(t::Real)
    _eval(b, t)[1]
end

@generated function differentiate(b::BernsteinBase{D, N, T}) where {D, N, T}
    exprs = Array{Expr}(undef, N-1)
    for k = 1:N-1
        exprs[k] = :((b[$(k+1)] - b[$k])*(N-1))
    end

    return quote
        @_inline_meta
        @inbounds elements = tuple($(exprs...))
        @inbounds return BernsteinBase(SVector{N-1}(elements))
    end
end

@generated function integrate(b::BernsteinBase{D, N, T}) where {D, N, T}
    exprs = Array{Expr}(undef, N+1)
    for k = 1:N+1
        exprs[k] = :(zero(b[1]))
        for j=1:k-1
            exprs[k] = :($(exprs[k]) + b[$j])
        end
        exprs[k] = :($(exprs[k])/N)
    end

    return quote
        @_inline_meta
        @inbounds elements = tuple($(exprs...))
        @inbounds return BernsteinBase(SVector{N+1}(elements))
    end
end

function *(k::Real, b::BernsteinBase)
    BernsteinBase(b.control_points*k)
end
*(b::BernsteinBase, k::Real) = k*b
/(b::BernsteinBase, k::Real) = (1/k)*b
-(b::BernsteinBase) = (-1)*b
+(b::BernsteinBase) = (+1)*b

function +(a::BernsteinBase{D, Na}, b::BernsteinBase{D, Nb}) where {D, Na, Nb}
    if Na < Nb
        a = _degelevate(a)
    elseif Na > Nb
        b = _degelevate(b)
    end
    a + b
end

function +(a::BernsteinBase{D, N}, b::BernsteinBase{D, N}) where {D, N}
    BernsteinBase(a.control_points + b.control_points)
end

function -(a::BernsteinBase{D}, b::BernsteinBase{D}) where {D}
    a + -b
end

@generated function _degelevate(b::BernsteinBase{D, N, T}) where {D, N, T}
    exprs = Array{Expr}(undef, N+1)
    exprs[1] = :(b[1])
    exprs[N+1] = :(b[N])
    for k = 2:N
        exprs[k] = :($(k/N)*b[$(k-1)] + (1 - $(k/N))*b[$k])
    end

    return quote
        @_inline_meta
        @inbounds elements = tuple($(exprs...))
        @inbounds return BernsteinBase(SVector{N+1}(elements))
    end
end

@generated function _binomial(::Type{BernsteinBase{D, N, T}}) where {D, N, T}
    nck = Array{Expr}(undef, N)
    for i = 1:N
        nck[i] = (i==1) ? :(one($T)) : :($(nck[i-1])*($N-$i+1)/($i-1))
    end
    return nck
end

@generated function _magsquared(b::BernsteinBase{D, N, T}) where {D, N, T}
    exprs = Array{Expr}(undef, 2N-1)
    mck = Array{Expr}(undef, N)
    den = :nothing

    nck = _binomial(b)

    for k = 1:2N-1
        exprs[k] = :(zero($T))
        den = (k==1) ? :(one($T)) : :($den*$(2N-k)/$(k-1))
        for j = max(1, k-N+1):min(N, k)
            exprs[k] = :($(exprs[k]) + $(nck[j])*$(nck[k-j+1])*b[$j]'*b[$(k-j+1)])
        end
        exprs[k] = :(SVector{1}($(exprs[k])/$den))
    end

    return quote
        @_inline_meta
        @inbounds elements = tuple($(exprs...))
        @inbounds return BernsteinBase(SVector{2N-1}(elements))
    end
end

struct Bernstein{D, N, M, T} <: Curve{D, T}
    f::BernsteinBase{D, N, T}
    s::BernsteinBase{1, M, T}
    limits::Interval{T}
end

function Bernstein(cp::Vararg{SVector,N}; limits=Interval(0., 1.)) where {N}
    cpts = SVector{N}(cp...)
    f = BernsteinBase(cpts)
    df = differentiate(f)
    s = integrate(_magsquared(df))
    Bernstein(f, s, limits)
end

function Bernstein(f::BernsteinBase; limits=Interval(0., 1.)) where {N}
    df = differentiate(f)
    s = integrate(_magsquared(df))
    Bernstein(f, s, limits)
end

function Bernstein(data::Array{Array{T, 1}}; limits=Interval(0., 1.)) where {T}
    Bernstein(BernsteinBase(data), limits=limits)
end

scalet(b::Bernstein, t) = (t - b.limits.lo)/diam(b.limits)

function (b::Bernstein)(t::Real)
    t = scalet(b, t)
    b.f(t)
end

function arclength(b::Bernstein, t::Real)
    t = scalet(b, t)
    s = b.s(t)
    s > 0 ? sqrt(t*s) : 0.0
end

function arclength(b::Bernstein, θ::Interval)
    θ = scalet(b, θ)
    s = b.s(θ.hi) - b.s(θ.lo)
    s > 0 ? sqrt(diam(θ)*s) : 0.0
end

function arclength(b::Bernstein)
    s = b.s(1.0)
    s > 0 ? sqrt(s) : 0.0
end

Base.eltype(::Type{Bernstein{D, N, T}}) where {D, N, T} = T
Base.getindex(b::Bernstein, j::Int) = b.f.control_points[j]

function Base.show(io::IO, b::Bernstein{D, N, T}) where {D, N, T}
    ordind = (N-1)%10 == 1 ? "st" : "th"
    ordind = (N-1)%10 == 2 ? "nd" : ordind
    ordind = (N-1)%10 == 3 ? "rd" : ordind
    ordind = 11 ≤ (N-1)%100 ≤ 13 ? "th" : ordind

    print(io, "a ", (N-1), ordind, " order Bernstein polynomial with control points at:\n",
          b.f.control_points.data, "\nwith an arclength of ", arclength(b, b.limits.hi))
end

function *(k::Real, b::Bernstein)
    Bernstein(b.f*k, limits=b.limits)
end
*(b::Bernstein, k::Real) = k*b
/(b::Bernstein, k::Real) = (1/k)*b
-(b::Bernstein) = (-1)*b
+(b::Bernstein) = (+1)*b

function +(a::Bernstein{D}, b::Bernstein{D}) where {D}
    @assert a.limits == b.limits
    c = a.f + b.f
    Bernstein(c, limits=a.limits)
end

function -(a::Bernstein{D}, b::Bernstein{D}) where {D}
    a + -b
end

function differentiate(b::Bernstein)
    c = differentiate(b.f)/diam(b.limits)
    Bernstein(c, limits=b.limits)
end

import Random: rand
function rand(::Type{BernsteinBase{D, N}}) where {D, N}
    pts = [@SVector(rand(D)) for i=1:N]
    BernsteinBase(pts...)
end

function rand(::Type{Bernstein{D, N}}; limits=Interval(0., 1.)) where {D, N}
    bern = rand(BernsteinBase{D, N})
    Bernstein(bern, limits=limits)
end


text:
using Jecco
using Jecco.KG_3_1

par_base = ParamBase(
    which_potential = "square",
)

par_grid = ParamGrid(
    xmin        = -5.0,
    xmax        =  5.0,
    xnodes      =  128,
    ymin        = -5.0,
    ymax        =  5.0,
    ynodes      =  128,
    umin        =  0.0,
    umax        =  1.0,
    udomains    =  4,
    unodes      =  16,
)

par_id = ParamID(
    ID_type      = "sine2D",
    A0x          = 1.0,
    A0y          = 1.0,
    Lx           = par_grid.xmax - par_grid.xmin,
    Ly           = par_grid.ymax - par_grid.ymin,
)

par_evol = ParamEvol(
    ODE_method = "AB3",
    # ODE_method = "RK4",
    dt      = 0.008, # for AB3. runs faster than with RK4
    # tmax    = 4.0,
    tmax    = 1.0,
)

par_io = ParamIO(
    out_every   = 4,
    folder      = "./data",
    overwrite   = true,
)

Jecco.KG_3_1.run_model(par_base, par_grid, par_id, par_evol, par_io)


text:
# This file is a part of DensityInterface.jl, licensed under the MIT License (MIT).

"""
    DensityInterface.test_density_interface(object, x, ref_logd_at_x; kwargs...)

Test that `object` is compatible with `DensityInterface`.

Tests that either `DensityKind(object) isa IsOrHasDensity`.

Also tests that [`logdensityof(object, x)`](@ref) equals `ref_logd_at_x` and
that the behavior of [`logdensityof(object)`](@ref),
[`densityof(object, x)`](@ref) and [`densityof(object)`](@ref) is consistent.

The results of `logdensityof(object, x)` and `densityof(object, x)` are compared to
`ref_logd_at_x` and `exp(ref_logd_at_x)` using `isapprox`. `kwargs...` are
forwarded to `isapprox`.

Also tests that `d = logfuncdensity(logdensityof(object))` returns a density
(`DensityKind(d) == IsDensity()`) that is equivalent to `object` in respect to
`logdensityof` and `densityof`, and that `funcdensity(densityof(object))`
behaves the same way.
"""
function test_density_interface(object, x, ref_logd_at_x; kwargs...)
    @testset "test_density_interface: $object with input $x" begin
        ref_d_at_x = exp(ref_logd_at_x)

        @test DensityKind(object) isa IsOrHasDensity

        @test isapprox(logdensityof(object, x), ref_logd_at_x; kwargs...)
        log_f = logdensityof(object)
        @test isapprox(log_f(x), ref_logd_at_x; kwargs...)

        @test isapprox(densityof(object,x), ref_d_at_x; kwargs...)
        f = densityof(object)
        @test isapprox(f(x), ref_d_at_x; kwargs...)

        for d in (logfuncdensity(log_f), funcdensity(f))
            @test DensityKind(d) == IsDensity()
            @test isapprox(logdensityof(d, x), ref_logd_at_x; kwargs...)
            @test isapprox(logdensityof(d)(x), ref_logd_at_x; kwargs...)
            @test isapprox(densityof(d,x), ref_d_at_x; kwargs...)
            @test isapprox(densityof(d)(x), ref_d_at_x; kwargs...)
        end
    end
end


text:
function seed(s = 1234)
    """
    function: generates a seed
    #TODO: can it be nothing?

    # TODO:generate a seed

    """

    seed = s #or

    return seed

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
@testset "1556.thousand-separator.jl" begin
    @test thousand_separator(123456789) == "123.456.789"
    @test thousand_separator(123) == "123"
    @test thousand_separator(1234) == "1.234"
    @test thousand_separator(0) == "0"
end


text:
"""
    Assembly{TF, TP<:AbstractVector{<:AbstractVector{TF}}, TC<:AbstractVector{<:Integer}, TE<:AbstractVector{Element{TF}}}

Composite type that defines an assembly of connected nonlinear beam elements for
analysis.

# Fields
 - `points`: Array of all beam element endpoints
 - `start`: Array containing point index where each beam element starts
 - `stop`: Array containing point index where each beam element stops
 - `element`: Array of `Element`s
"""
struct Assembly{TF, TP<:AbstractVector{<:AbstractVector{TF}}, TC<:AbstractVector{<:Integer}, TE<:AbstractVector{Element{TF}}}
    points::TP
    start::TC
    stop::TC
    elements::TE
end
Base.eltype(::Assembly{TF, TP, TC, TE}) where {TF, TP, TC, TE} = TF

"""
    Assembly(points, start, stop; kwargs...)

Construct an assembly of connected nonlinear beam elements for analysis.  Beam lengths
and midpoints may be manually specified in case beam elements are curved rather than
straight.

# Arguments
 - `points`: Array of all beam element endpoints
 - `start`: Array containing point indices where each beam element starts
 - `stop`: Array containing point indices where each beam element stops

# Keyword Arguments
 - `stiffness`: Array of (6 x 6) stiffness matrices for each beam element, alternative to providing `compliance`
 - `compliance`: Array of (6 x 6) compliance matrices for each beam element, defaults to `zeros(6,6)` for each beam element
 - `mass`: Array of (6 x 6) mass matrices for each beam element, alternative to providing `minv`
 - `minv`: Array of (6 x 6) mass matrices for each beam element, defaults to the identity matrix for each beam element
 - `frames`: Array of (3 x 3) rotation matrices for each beam element (to transform
        to the global frame), defaults to the identity matrix for each beam element
 - `lengths`: Array containing the length of each beam, defaults to the distance between beam endpoints
 - `midpoints`: Array containing the midpoint of each beam element, defaults to the average of the beam element endpoints
"""
function Assembly(points, start, stop;
    stiffness = nothing,
    compliance = nothing,
    mass = nothing,
    minv = nothing,
    frames = nothing,
    lengths = norm.(points[stop] - points[start]),
    midpoints = (points[stop] + points[start])/2)

    nbeam = length(start)

    if isnothing(compliance)
        if isnothing(stiffness)
            compliance = fill((@SMatrix zeros(6,6)), nbeam)
        else
            compliance = fill((@MMatrix zeros(6,6)), nbeam)
            for i = 1:nbeam
                filled_cols = findall(vec(mapslices(col -> any(row -> !isapprox(row, 0), col), stiffness[i], dims = 1)))
                compliance[i][filled_cols,filled_cols] .= inv(Matrix(stiffness[i][filled_cols, filled_cols]))
            end
            compliance = SMatrix.(compliance)
        end
    end

    if isnothing(minv)
        if isnothing(mass)
            minv = fill(Diagonal(@SVector ones(6)), nbeam)
        else
            minv = fill(MMatrix{6,6}(Diagonal(@SVector ones(eltype(eltype(mass)), 6))), nbeam)
            for i = 1:nbeam
                filled_cols = findall(vec(mapslices(col -> any(row -> !isapprox(row, 0), col), mass[i], dims = 1)))
                minv[i][filled_cols,filled_cols] .= inv(Matrix(mass[i][filled_cols, filled_cols]))
            end
            minv = SMatrix{6,6}.(minv)
        end
    end

    if isnothing(frames)
        frames = fill(I3, nbeam)
    end

    elements = Element.(lengths, midpoints, compliance, minv, frames)

    return Assembly(points, promote(start, stop)..., elements)
end

"""
    discretize_beam(L, start, discretization; frame = Matrix(I,3,3)), curvature = zeros(3))

Discretize a beam according to the discretization provided in `discretization`
given the beam length (`L`), and starting point (`start`).

Return the lengths, endpoints, midpoints, and rotation matrices of the beam elements.

# Arguments
 - `L`: Beam length
 - `start`: Beam starting point
 - `discretization`: May be either an integer, representing the number of
        elements that the beam should be discretized into, or a vector containing
        the normalized endpoints of each beam element, where 0 is the beginning
        of the beam and 1 is the end of the beam.
 - `frame`: 3x3 beam rotation matrix which transforms from the local beam
        coordinate frame at the start of the beam to the global coordinate frame.
 - `curvature`: curvature vector
"""
discretize_beam(L, start, discretization::Integer; frame = I3, curvature = (@SVector zeros(3))) =
    discretize_beam(L, start, range(0, 1, length=discretization+1); frame=frame, curvature=curvature)

function discretize_beam(L, start, discretization; frame = I3, curvature = (@SVector zeros(3)))

    r1 = SVector{3}(start)
    Cab = SMatrix{3,3}(frame)
    k = SVector{3}(curvature)

    # discretize beam
    sp = L*discretization
    sm = (sp[1:end-1] .+ sp[2:end])/2
    ΔL = (sp[2:end] .- sp[1:end-1])

    # precompute some curvature quantities
    kkt = k*k'
    ktilde = tilde(k)
    kn = sqrt(k'*k)

    if curvature == zero(curvature)
        triads = fill(Cab, length(sm))
        xp = [r1 + s*Cab*e1 for s in sp]
        xm = [r1 + s*Cab*e1 for s in sm]
    else
        triads = curve_triad.(Ref(Cab), Ref(kkt), Ref(ktilde), Ref(kn), sm)
        xp = curve_coordinates.(Ref(r1), Ref(Cab), Ref(kkt), Ref(ktilde), Ref(kn), sp)
        xm = curve_coordinates.(Ref(r1), Ref(Cab), Ref(kkt), Ref(ktilde), Ref(kn), sm)
    end

    return ΔL, xp, xm, triads
end

"""
    curve_length(start, stop, curvature)

Calculate the length of a curve given its endpoints and its curvature vector
"""
function curve_length(start, stop, curvature)
    r1, r2, k = SVector{3}(start), SVector{3}(stop), SVector{3}(curvature)
    kn = sqrt(k'*k)
    r = r2 - r1
    rn = sqrt(r'*r)
    if kn == 0 || (k[2] == 0 && k[3] == 0)
        ΔL = rn
    else
        if k[1] == 0
            ΔL = 2*asin(kn*rn/2)/kn
        else
            lower_bound = rn
            upper_bound = rn*kn/k
            residual = (L) -> (2*(kn^2-k[1]^2)*(1-cos(kn*L)) + k[1]^2*(kn*L)^2)/kn^4 - rnorm
            ΔL = fzero(residual, lower_bound, upper_bound)
        end
    end
    return ΔL
end

"""
    curve_triad(Cab, k, s)
    curve_triad(Cab, kkt, ktilde, kn, s)

Return the rotation matrix at `s` along the length of the beam given
the curvature vector `k` and the initial rotation matrix `Cab`.
"""
@inline curve_triad(Cab, k, s) = curve_triad(Cab, k*k', tilde(k), sqrt(k'*k), s)
@inline curve_triad(Cab, kkt, ktilde, kn, s) = SMatrix{3,3}(Cab*((I - kkt/kn^2)*cos(kn*s) +
    ktilde/kn*sin(kn*s) + kkt/kn^2))

"""
    curve_coordiantes(r, Cab, k, s)
    curve_coordinates(r, Cab, kkt, ktilde, kn, s)

Return the coordinates at `s` along the length of the beam given the starting
point `r`, initial rotation matrix `Cab`, and curvature vector `k`.
"""
@inline curve_coordinates(r, Cab, k, s) = curve_coordinates(r, Cab, k*k', tilde(k), sqrt(k'*k), s)
@inline curve_coordinates(r, Cab, kkt, ktilde, kn, s) = r + SVector{3}(Cab*((I/kn - kkt/kn^2)*sin(kn*s) + ktilde/kn^2*(1-cos(kn*s)) + kkt/kn^2*s)*e1)


text:
function checkboxes(; options, value, multiple = true)
    if multiple
        update = js"""
        var vec = data.value;
        var val = option.value;
        var idx = vec.indexOf(val);
        idx >= 0 ? vec.splice(idx, 1) : vec.push(val);
        data.value = vec;
        """
        check = js"data.value.indexOf(option.value) >= 0"
    else
        update = js"data.value = option.value;"
        check = js"option.value == data.value"
    end
    template = js"""
    {
        view: function () {
            var list = data.options;
            var children = list.map(function (option) {
                return m("div.field", [
                    m("input", {
                        type: "checkbox",
                        value: option.value,
                        checked: $check,
                        onchange: function (e) {
                            $update
                        }}),
                    m("label", option.name)
                ])
            })
            return m("div", children);
        }
    }
    """
    MithrilWidget{:checkboxes}(MithrilComponent(template, (options = options, value = value)), value)
end


text:
#!/usr/bin/env julia

# Dipole chain DMRG example.

using LatticeSweeper

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "-R"
        metavar = "R"
        help = "separation distance"
        arg_type = Float64
        required = true
    "-N"
        metavar = "N"
        help = "number of rotors"
        arg_type = Int
        required = true
    "--l-max"
        metavar = "L"
        help = "local basis truncation"
        arg_type = Int
        required = true
    "--max-sweeps"
        metavar = "S"
        help = "maximum number of sweeps"
        arg_type = Int
        required = true
    "--output-file"
        metavar = "PATH"
        help = "path to output file"
        arg_type = String
end
c = parse_args(ARGS, s, as_symbols=true)

R = c[:R]
N = c[:N]
l_max = c[:l_max]
max_sweeps = c[:max_sweeps]
output_file = c[:output_file]

# l m = 0 0, 1 -1, 1 0, 1 1, ...
basis = [(l, m) for l in 0:l_max for m in -l:l]
basis_size = length(basis)

# One-site operators.
id = zeros(basis_size, basis_size)
L2 = zeros(basis_size, basis_size)
Lp = zeros(basis_size, basis_size)
Lm = zeros(basis_size, basis_size)
LpMp = zeros(basis_size, basis_size)
LpMm = zeros(basis_size, basis_size)
LmMp = zeros(basis_size, basis_size)
LmMm = zeros(basis_size, basis_size)
for (col, (l, m)) in enumerate(basis)
    for (row, (lp, mp)) in enumerate(basis)
        if lp == l && mp == m
            id[row, col] = 1.0
            L2[row, col] = l*(l+1)
        end
        if lp == l+1 && mp == m
            Lp[row, col] = sqrt((l-m+1)*(l+m+1)/((2l+1)*(2l+3)))
        end
        if lp == l-1 && mp == m
            Lm[row, col] = sqrt((l-m)*(l+m)/((2l-1)*(2l+1)))
        end
        if lp == l+1 && mp == m+1
            LpMp[row, col] = -sqrt((l+m+1)*(l+m+2)/((2l+1)*(2l+3)))
        end
        if lp == l+1 && mp == m-1
            LpMm[row, col] = -sqrt((l-m+1)*(l-m+2)/((2l+1)*(2l+3)))
        end
        if lp == l-1 && mp == m+1
            LmMp[row, col] = sqrt((l-m-1)*(l-m)/((2l-1)*(2l+1)))
        end
        if lp == l-1 && mp == m-1
            LmMm[row, col] = sqrt((l+m-1)*(l+m)/((2l-1)*(2l+1)))
        end
    end
end

# Dipole-dipole potential operator.
ks = Float64[]
op1s = Matrix{Float64}[]
op2s = Matrix{Float64}[]
for (k, OP1, OP2) in [(-2.0, [Lp, Lm], [Lp, Lm]),
                      (-0.5, [LpMm, LmMm], [LpMp, LmMp]),
                      (-0.5, [LpMp, LmMp], [LpMm, LmMm])]
    for op1 in OP1
        for op2 in OP2
            push!(ks, k)
            push!(op1s, op1)
            push!(op2s, op2)
        end
    end
end
K = length(ks)

println("[ ] Constructing Hamiltonian.")
H_tnsr = zeros(2+(N-1)*K, 2+(N-1)*K, basis_size, basis_size)
H_tnsr[1, 1, :, :] = id
H_tnsr[2, 1, :, :] = L2
H_tnsr[2, 2, :, :] = id
for i in 1:(N-1)
    for k in 1:K
        H_tnsr[2, 2+(i-1)*K+k, :, :] = ks[k]*op1s[k]/(R*i)^3
    end
end
for k in 1:K
    H_tnsr[2+k, 1, :, :] = op2s[k]
end
for i in 1:(N-2)
    for k in 1:K
        H_tnsr[2+i*K+k, 2+(i-1)*K+k, :, :] = id
    end
end
@time H = MPO(H_tnsr, N)
@time mpo_mean_bond_dim = compress!(H, 1e-15)
println("[+] Constructed Hamiltonian: $(mpo_mean_bond_dim).")

println("[ ] Constructing wavefunction.")
wf = zeros(basis_size)
wf[1] = 1.0
@time psi = MPS(wf, N)
println("[+] Constructed wavefunction.")

println("[ ] Sweeping.")
outputs = SweepOutput[]
if !isnothing(output_file)
    push!(outputs, SweepOutputFile(output_file))
end
push!(outputs, SweepOutputDynamic())
@time hist = dmrg!(psi, H, SweepSchedule(max_sweeps; tolerance=1e-7); outputs=outputs)
println("[+] Swept.")

hist.converged || @warn("Ground state not converged.")
if hist[end].mps_max_bond_dim >= hist[end].set.bond_max
    @warn("Bond dimension saturated at final sweep.")
end

# Ground state energy.
println_result("E0 = $(hist[end].energy)")
# Von Neumann entanglement entropy.
println_result("SvN = $(S_vn(hist[end].middle_eigvals))")


text:
using Printf

# Source: https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/src/trials.jl
function prettytime(t)
    if t < 1e3
        value, units = t, "ns"
    elseif t < 1e6
        value, units = t / 1e3, "μs"
    elseif t < 1e9
        value, units = t / 1e6, "ms"
    else
        s = t / 1e9
        if s < 60
            value, units = s, "s"
        else
            value, units = (s / 60), "min"
        end
    end
    return string(@sprintf("%.3f", value), " ", units)
end

function Base.zeros(T, ::GPU, g)
    a = CuArray{T}(undef, g.Nx, g.Ny, g.Nz)
    a .= 0
    return a
end

Base.zeros(T, ::CPU, g) = zeros(T, size(g))

# Default to type of Grid
Base.zeros(arch, g::Grid{T}) where T = zeros(T, arch, g)


text:
using CMPFit, GFit, Gnuplot
using QSFit, DataStructures, Statistics, Dierckx

import QSFit: known_spectral_lines, add_qso_continuum!, LineComponent

export T2Recipe

abstract type T2Recipe <: DefaultRecipe end

function known_spectral_lines(source::QSO{T}) where T <: T2Recipe
    list = [
        NarrowLine(                      :Lyb                       ),
        NarrowLine(                      :Lya                       ),
        NarrowLine(                      :NV_1241                   ),
        NarrowLine(                      :CIV_1549                  ),
        NarrowLine(                      :CIII_1909                 ),
        NarrowLine(                      :MgII_2798                 ),
        NarrowLine(                      :NeV_3426                  ),
        NarrowLine(                      :OII_3727a                 ),
        NarrowLine(                      :OII_3727b                 ),
        NarrowLine(                      :NeIII_3870                ),
        NarrowLine(                      :Hg                        ),
        NarrowLine(                      :Hb                        ),
        NarrowLine(                      :OIII_4959                 ),
        NarrowLine(                      :OIII_5007                 ),
        AsymmTailLine(                   :OIII_5007      ,     :blue),
        NarrowLine(                      :OI_6300                   ),
        NarrowLine(                      :OI_6364                   ),
        NarrowLine(                      :NII_6549                  ),
        NarrowLine(                      :Ha                        ),
        NarrowLine(                      :NII_6583                  ),
        NarrowLine(                      :SII_6716                  ),
        NarrowLine(                      :SII_6731                  )]
    return list
end

function add_qso_continuum(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: T2Recipe
    λ = domain(model)[:]

    comp = QSFit.powerlaw(median(λ))
    comp.alpha.val = -1.8

    model[:qso_cont] = comp
    push!(model[:Continuum].list, :qso_cont)
    evaluate!(model)
end

function LineComponent(source::QSO{T}, line::NarrowLine, multicomp::Bool) where T <: T2Recipe
    lc = LineComponent(parent_recipe(source), line, multicomp) # invoke parent recipe
    lc.comp.fwhm.low  = 10
    return lc
end

function default_unk_line(source::QSO{T}) where T <: T2Recipe
    comp = SpecLineGauss(5e3)
    comp.norm.val = 0.
    comp.center.fixed = false
    comp.center.low = 0
    comp.center.high = Inf
    comp.fwhm.val  = 500
    comp.fwhm.low  = 10
    comp.fwhm.high = 2e3
    comp.voff.fixed = true
    return comp
end


text:
# import Pkg; Pkg.add(["FiniteDifferences", "Plots"])
using Faust, FiniteDifferences, Plots

process = compile("""
import("stdfaust.lib");

process = pm.ks(pm.f2l(ba.midikey2hz(60)), 0.1);
""")

function f(x)
    init!(process, block_size=size(x, 1))
    process.inputs = x
    compute!(process)
end

x = randn(1024, 1)
∇f = jacobian(central_fdm(5, 1), f, x)[1]
heatmap(∇f)


text:
module runexamples

import Modia3D

# Test Modia3D - Examples

# examples/dynamics
dynamicsPath = joinpath(Modia3D.path, "examples", "dynamics")
include(joinpath(dynamicsPath, "Simulate_DoublePendulum.jl"))
include(joinpath(dynamicsPath, "Simulate_DoublePendulumWithDampers.jl"))
include(joinpath(dynamicsPath, "Simulate_FallingBall1.jl"))
#include(joinpath(dynamicsPath, "Simulate_FallingBall2.jl"))
#include(joinpath(dynamicsPath, "Simulate_FallingBall3.jl"))
#include(joinpath(dynamicsPath, "Simulate_FallingBall4.jl"))
#include(joinpath(dynamicsPath, "Simulate_FallingBall4.jl"))
#include(joinpath(dynamicsPath, "Simulate_FallingBodies.jl"))
include(joinpath(dynamicsPath, "Simulate_Pendulum.jl"))
include(joinpath(dynamicsPath, "Simulate_Pendulum_explained.jl"))
include(joinpath(dynamicsPath, "Simulate_Pendulum_withoutInertia.jl"))
include(joinpath(dynamicsPath, "Simulate_PendulumWithController.jl"))
include(joinpath(dynamicsPath, "Simulate_PendulumWithDamper.jl"))

# examples/kinematics
kinematicsPath = joinpath(Modia3D.path, "examples", "kinematics")
include(joinpath(kinematicsPath, "Move_DoublePendulum.jl"))
include(joinpath(kinematicsPath, "Move_FourBar.jl"))
include(joinpath(kinematicsPath, "Move_FourBar2.jl"))
#include(joinpath(kinematicsPath, "Move_FourBar_Collision1.jl"))
#include(joinpath(kinematicsPath, "Move_FourBar_Collision2.jl"))

# examples/signalToFlange
signalToFlangePath = joinpath(Modia3D.path, "examples", "SignalToFlange")
include(joinpath(signalToFlangePath, "Test_SignalAngle.jl"))
include(joinpath(signalToFlangePath, "Test_SignalTorque.jl"))

# examples/visual
visualPath = joinpath(Modia3D.path, "examples", "visual")
include(joinpath(visualPath, "Move_AllVisualObjects.jl"))
include(joinpath(visualPath, "Move_SolidFileMesh.jl"))
include(joinpath(visualPath, "Visualize_AllVisualObjects.jl"))
include(joinpath(visualPath, "Visualize_Assembly.jl"))
include(joinpath(visualPath, "Visualize_GeometriesWithMaterial.jl"))
include(joinpath(visualPath, "Visualize_GeometriesWithoutMaterial.jl"))
include(joinpath(visualPath, "Visualize_SolidFileMesh.jl"))
include(joinpath(visualPath, "Visualize_Solids.jl"))
include(joinpath(visualPath, "Visualize_Text.jl"))
include(joinpath(visualPath, "Visualize_TextFonts.jl"))

println("\n... success of runexamples.jl")
end


text:
function export_matpowerm_instance!(db::SQLite.DB, export_dir, name, scenario, source_type, source_path, pfn_path)
    println("Exporting instance to MATPOWERM file: ($name, $scenario)")
    export_path = ""
    if source_type == "MATPOWERM"
        export_path = joinpath(export_dir, basename(source_path))
        cp(source_path, export_path; force=true)
    else
        network = ismissing(pfn_path) ? PowerFlowNetwork(source_path, source_type) : deserialize(pfn_path)
        export_path = joinpath(export_dir, "$(name)_$(scenario).m")
        write_pfn(export_path, network)
    end
    query = "UPDATE instances SET matpowerm_path = '$export_path' WHERE name = '$name' AND scenario = $scenario"
    execute_query(db, query)
end

function export_matpowerm_instance_dfrow!(db::SQLite.DB, export_dir, row)
    export_matpowerm_instance!(db, export_dir, row[:name], row[:scenario], row[:source_type], row[:source_path], row[:pfn_path])
end

function export_matpowerm_instances!(db::SQLite.DB; min_nv, max_nv, export_dir, recompute=false, subset=nothing)
    !isdir(export_dir) && mkpath(export_dir)
    query = "SELECT name, scenario, source_type, source_path, pfn_path FROM instances WHERE nbus >= $min_nv AND nbus <= $max_nv"
    if !recompute
        query *= " AND matpowerm_path IS NULL"
    end
    if !isnothing(subset)
        query *= " AND id IN ($(join(subset, ',')))"
    end
    results = DBInterface.execute(db, query) |> DataFrame
    println("Exporting instance to MATPOWERM file: recompute=$recompute, min_nv=$min_nv, max_nv=$max_nv, export_dir=$export_dir")
    println("subset\n$subset\nend subset")
    export_function!(row) = export_matpowerm_instance_dfrow!(db, export_dir, row)
    export_function!.(eachrow(results[!, [:name, :scenario, :source_type, :source_path, :pfn_path]]))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
@testset "periodic boundary condition" begin
    @test pbc(0, 3) == 3
    @test pbc(1, 3) == 1
    @test pbc(2, 3) == 2
    @test pbc(3, 3) == 3
    @test pbc(4, 3) == 1
    @test pbc([0, 1, 2, 3, 4] ,3) == [3, 1, 2, 3, 1]
    @test pbc(-7:7, 3) == [2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1]
end


text:
#Demonstrate Saramaki half-band filter design
using RSDeltaSigmaPort
using RSDeltaSigmaPort.EasyPlot #set, cons
import RSDeltaSigmaPort: csdsize2, linlin
import RSDeltaSigmaPort: BoundingBox
import RSDeltaSigmaPort: fft
import Printf: @sprintf
j=im

use_canned_example = false
println("\n*** Half-band filter design (ALPHA VERSION)")
f1=f2=0 #Define variables
title=nothing
if use_canned_example
	(f1, f2) = exampleHBF(2)
else
	fp = 0.9*0.25
	@warn("Need firls() equivalent to converge for small dB values.")
	#delta = undbv( -100 )
	delta = undbv( -25 )
	title = "designHBF Iterations"
	(f1, f2, info) = designHBF(fp, delta=delta, debug=true)
end
n1 = length(f1); n2 = length(f2)
complexity = sum(csdsize2(f1)) + (2*n1-1)*(n2+sum(csdsize2(f2))-1) #VERIFYME
@show complexity

println()
@info("Running simulations..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
#interleave the even and odd decimated impulse responses
Nimp = 2^11
imp = simulateHBF(vcat(1, zeros(Nimp-1)),f1,f2)
mag = abs.(fft(imp))
mag = mag[1:array_round(end/2+1)]
t = range(0, stop=0.5, length=length(mag))
println("\tdone.")

println()
@info("Plotting resulting filter"); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
plot = cons(:plot, linlin, title="designHBF() result", legend=false,
	xyaxes=set(xmin=0, xmax=0.5, ymin=-150, ymax=10),
	labels=set(xaxis="Normalized Frequency", yaxis="|H| [dB]"),
)
magdBw = waveform(t, dbv.(mag))
push!(plot,
	cons(:wfrm, magdBw, line=set(style=:solid, color=:blue, width=2), label=""),
)
displaygui(plot)

:END_OF_DEMO


text:
using UniqueVectors: UniqueVector

export Blackbox
"""
    Blackbox{T, U}(ω::Vector{Float64}, Y::Vector{U}, P::U, Q::U,
        ports::UniqueVector{T}) where {T, U}
    Blackbox(ω::Vector{<:Real}, Y::Vector{U}, P::U, Q::U,
        ports::AbstractVector{T}) where {T, U}

An admittance model with constant `P` and `Q` but with varying `Y`, indexed by a real
variable `ω`.
"""
struct Blackbox{T, U<:AbstractMatrix{<:Number}} <: AdmittanceModel{T}
    ω::Vector{Float64}
    Y::Vector{U}
    P::U
    Q::U
    ports::UniqueVector{T}
    function Blackbox{T, U}(ω::Vector{Float64}, Y::Vector{U}, P::U, Q::U,
        ports::UniqueVector{T}) where {T, U}
        if length(Y) > 0
            l = size(Y[1], 1)
            @assert size(P) == (l, length(ports))
            @assert size(Q) == (l, length(ports))
            @assert all([size(y) == (l, l) for y in Y])
        end
        @assert length(ω) == length(Y)
        return new{T, U}(ω, Y, P, Q, ports)
    end
end

Blackbox(ω::Vector{<:Real}, Y::Vector{U}, P::U, Q::U,
    ports::AbstractVector{T}) where {T, U} = Blackbox{T, U}(ω, Y, P, Q, UniqueVector(ports))

get_Y(bbox::Blackbox) = bbox.Y

get_P(bbox::Blackbox) = bbox.P

get_Q(bbox::Blackbox) = bbox.Q

get_ports(bbox::Blackbox) = bbox.ports

function partial_copy(bbox::Blackbox{T, U};
    Y::Union{Vector{V}, Nothing}=nothing,
    P::Union{V, Nothing}=nothing,
    Q::Union{V, Nothing}=nothing,
    ports::Union{AbstractVector{W}, Nothing}=nothing) where {T, U, V, W}
    Y = Y == nothing ? get_Y(bbox) : Y
    P = P == nothing ? get_P(bbox) : P
    Q = Q == nothing ? get_Q(bbox) : Q
    ports = ports == nothing ? get_ports(bbox) : ports
    return Blackbox(bbox.ω, Y, P, Q, ports)
end

function compatible(bboxes::AbstractVector{Blackbox{T, U}}) where {T, U}
    if length(bboxes) == 0
        return true
    end
    ω = bboxes[1].ω
    return all([ω == bbox.ω for bbox in bboxes[2:end]])
end

# TODO for sparse matrices use sparse linear solver instead of \ and /

function impedance_to_scattering(Z::AbstractMatrix{<:Number}, Z0::AbstractVector{<:Real})
    U = diagm(0 => Z0)
    return (Z + U) \ (Z - U)
end

function scattering_to_impedance(S::AbstractMatrix{<:Number}, Z0::AbstractVector{<:Real})
    U = diagm(0 => Z0)
    return U * (I + S) / (I - S)
end

function admittance_to_scattering(Y::AbstractMatrix{<:Number}, Z0::AbstractVector{<:Real})
    U = diagm(0 => Z0)
    return (I + Y * U) \ (I - Y * U)
end

function scattering_to_admittance(S::AbstractMatrix{<:Number}, Z0::AbstractVector{<:Real})
    V = diagm(0 => 1 ./ Z0)
    return (I - S) * ((I + S) \ V)
end

export impedance_matrices
"""
    impedance_matrices(bbox::Blackbox)

Find the impedance matrices `Z` with `y = Zx` for each `ω`.
"""
impedance_matrices(bbox::Blackbox) = [transpose(bbox.Q) * (Y \ collect(bbox.P))
    for Y in bbox.Y]

export admittance_matrices
"""
    admittance_matrices(bbox::Blackbox)

Find the admittance matrices `Y` with `Yy = x` for each `ω`.
"""
function admittance_matrices(bbox::Blackbox)
    if bbox.P == I # avoid unnecessary matrix inversion
        return map(collect, bbox.Y)
    else
        return [inv(Z) for Z in impedance_matrices(bbox)]
    end
end

export scattering_matrices
"""
    scattering_matrices(bbox::Blackbox, Z0::AbstractVector{<:Real})

Find the scattering matrices that map incoming waves to outgoing waves on transmission
lines with characteristic impedance `Z0` for each `ω`.
"""
function scattering_matrices(bbox::Blackbox, Z0::AbstractVector{<:Real})
    if bbox.P == I # avoid unnecessary matrix inversion
        return [admittance_to_scattering(Y, Z0) for Y in admittance_matrices(bbox)]
    else
        return [impedance_to_scattering(Z, Z0) for Z in impedance_matrices(bbox)]
    end
end

function canonical_gauge(bbox::Blackbox)
    n = length(bbox.ports)
    Y = admittance_matrices(bbox)
    P = Matrix{eltype(eltype(Y))}(I, n, n)
    Q = Matrix{eltype(eltype(Y))}(I, n, n)
    return Blackbox(bbox.ω, Y, P, Q, bbox.ports)
end

"""
    Blackbox(ω::Vector{<:Real}, pso::PSOModel)

Create a Blackbox for a PSOModel where `ω` represents angular frequency.
"""
function Blackbox(ω::Vector{<:Real}, pso::PSOModel)
    K, G, C = get_Y(pso)
    Y = [K ./ s + G + C .* s for s in 1im * ω]
    P = get_P(pso) * (1.0 + 0im)
    Q = get_Q(pso) * (1.0 + 0im)
    return Blackbox(ω, Y, P, Q, get_ports(pso))
end


text:
using EmpiricalCDFs: EmpiricalCDF
using MaximumLikelihoodPower: mle

@testset "ordinary RW, first return" begin
    nsamples = 10^4
    cutoff = 200
    first_return_actor = FirstReturnActor()
    ecdf_actor = ECDFValueActor(first_return_actor)
    for w in (WalkB(), WalkF())
        @time trial!(WalkPlan(w, ActorSet(StepLimitActor(10^10), first_return_actor)), SampleLoopActor(nsamples, ecdf_actor))
        alpha, err = mle(ecdf_actor[end - cutoff:end])
        @test isapprox(alpha, 3/2, atol=0.2)
    end
end

@testset "ordinary RW, mean std" begin
    nsteps = 10^4
    ntrials = 10^5
    for w in (WalkB(), WalkF())
        wp = WalkPlan(w, StepLimitActor(nsteps))
        ecdf = EmpiricalCDF{Int}()
        cdf_actor = ECDFActor(get_x, ecdf)
        trial_loop = SampleLoopActor(ntrials, cdf_actor)
        @time trial!(wp, trial_loop)
        ecdf =  ecdf ./ sqrt(nsteps)
        tol = 1e-2
        @test isapprox(std(ecdf), 1; atol=tol)
        @test isapprox(mean(ecdf), 0; atol=tol)
    end
end


text:
function include_parameters!(
    optimization_container::OptimizationContainer,
    constraint_infos::Vector{DeviceTimeSeriesConstraintInfo},
    param_reference::UpdateRef,
    expression_name::Symbol,
    multiplier::Float64 = 1.0,
)
    @assert model_has_parameters(optimization_container)
    time_steps = model_time_steps(optimization_container)
    names = [get_component_name(r) for r in constraint_infos]
    @debug "adding" param_reference "parameter"
    container =
        add_param_container!(optimization_container, param_reference, names, time_steps)
    param = get_parameter_array(container)
    mult = get_multiplier_array(container)
    expr = get_expression(optimization_container, expression_name)
    for t in time_steps, r in constraint_infos
        param[get_component_name(r), t] =
            add_parameter(optimization_container.JuMPmodel, r.timeseries[t])
        mult[get_component_name(r), t] = r.multiplier * multiplier
        add_to_expression!(
            expr,
            r.bus_number,
            t,
            param[get_component_name(r), t],
            r.multiplier * multiplier,
        )
    end
    return container
end

function include_parameters!(
    optimization_container::OptimizationContainer,
    constraint_infos::Vector{DeviceTimeSeriesConstraintInfo},
    param_reference::UpdateRef,
    multiplier::Float64 = 1.0,
)
    @assert model_has_parameters(optimization_container)
    time_steps = model_time_steps(optimization_container)
    names = [get_component_name(r) for r in constraint_infos]
    container =
        add_param_container!(optimization_container, param_reference, names, time_steps)
    param = get_parameter_array(container)
    mult = get_multiplier_array(container)
    for t in time_steps, r in constraint_infos
        param[get_component_name(r), t] =
            add_parameter(optimization_container.JuMPmodel, r.timeseries[t])
        mult[get_component_name(r), t] = r.multiplier * multiplier
    end
    return container
end


text:
#!/usr/bin/julia

# Trizen
# 17 January 2017
# https://github.com/trizen

# Generates a Mandelbrot-like set.

# See also:
#   https://en.wikipedia.org/wiki/Mandelbrot_set
#   https://trizenx.blogspot.ro/2017/01/mandelbrot-set.html

using Images

@inline function hsv2rgb(h, s, v)
    c = v * s
    x = c * (1 - abs(((h/60) % 2) - 1))
    m = v - c

    if h < 60
        r,g,b = (c, x, 0)
    elseif h < 120
        r,g,b = (x, c, 0)
    elseif h < 180
        r,g,b = (0, c, x)
    elseif h < 240
        r,g,b = (0, x, c)
    elseif h < 300
        r,g,b = (x, 0, c)
    else
        r,g,b = (c, 0, x)
    end

    (r + m), (b + m), (g + m)
end

#
## A specific function which maps the value of `c` to some other value.
#
@inline function f(c)
    1/c
end

function mandelbrot()

    w, h = 1000, 1000

    zoom  = 1         # the zoom factor
    moveX = 0         # the amount of shift on the x axis
    moveY = 0         # the amount of shift on the y axis

    L = 100           # the maximum value of |z|
    I = 30            # the maximum number of iterations

    img = zeros(RGB{Float64}, h, w)

    for x in 1:w, y in 1:h
        c = Complex(
            (2*x - w) / (w * zoom) + moveX,
            (2*y - h) / (h * zoom) + moveY
        )
        z = c
        n = 0
        q = f(c)
        while (abs(z) < L && (n += 1) < I)
            z = z^q
        end
        v = (I - n) / I
        r,g,b = hsv2rgb(v*360, 1, v)
        img[y,x] = RGB{Float64}(r, g, b)
    end

    println("Generating image...")
    save("mandelbrot_like_set.png", img)
end

mandelbrot()

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
export NwtdirectionCG, NewtonCG
using Krylov

function NwtdirectionCG(H, ∇f; kwargs...)
    e=1e-6
    n = length(∇f)
    τ = 0.5 # need parametrization
    cgtol = max(e, min(0.7, 0.01 * norm(∇f)^(1.0 + τ)))

    (d, cg_stats) = cg(H, -∇f,
                       atol=cgtol, rtol=0.0,
                       itmax=max(2 * n, 50),
                       linesearch = true,
                       kwargs...)

    return d
end

function NewtonCG(nlp :: AbstractNLPModel;
                  x :: Vector{T}=copy(nlp.meta.x0),
                  stp :: NLPStopping = NLPStopping(nlp,
                                                   NLPAtX(nlp.meta.x0)),
                  kwargs...) where T
    return  Newton_G_Stop(nlp;
                          x = x,
                          stp = stp,
                          NwtDirection = NwtdirectionCG,
                          hessian_rep = hessian_operator,
                          kwargs...)
end

function NewtonCGLS(nlp :: AbstractNLPModel;
                  x :: Vector{T}=copy(nlp.meta.x0),
                  stp :: NLPStopping = NLPStopping(nlp,
                                                   NLPAtX(nlp.meta.x0)),
                  LS_algo   :: Function = bracket,
                  LS_logger :: AbstractLogger = Logging.NullLogger(),
                  kwargs...) where T
    return  Newton_G_StopLS(nlp;
                          x = x,
                          stp = stp,
                          LS_algo = LS_algo,
                          LS_logger = LS_logger,
                          NwtDirection = NwtdirectionCG,
                          hessian_rep = hessian_operator,
                          kwargs...)
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
Base.@kwdef struct PackedManifoldKernelDensity <: PackedSamplableBelief
  _type::String            = "IncrementalInference.PackedManifoldKernelDensity"
  varType::String
  pts::Vector{Vector{Float64}}
  bw::Vector{Float64}      = Float64[]
  partial::Vector{Int}     = Int[]
  infoPerCoord::Vector{Float64} = zeros(length(pts[1]))
end

Base.@kwdef struct PackedAliasingScalarSampler <: PackedSamplableBelief
  _type::String           = "IncrementalInference.PackedAliasingScalarSampler"
  domain::Vector{Float64} = [0;1.0;]
  weights::Vector{Float64}= [0.5;0.5;]
end

Base.@kwdef mutable struct PackedHeatmapGridDensity <: PackedSamplableBelief
  _type::String = "IncrementalInference.PackedHeatmapGridDensity"
  data::Vector{Vector{Float64}}
  domain::Tuple{Vector{Float64}, Vector{Float64}}
  hint_callback::String
  bw_factor::Float64
  N::Int
  # _densityFnc::String = "" # only use if storing parched belief data entry label/id
end

Base.@kwdef mutable struct PackedLevelSetGridNormal <: PackedSamplableBelief
  _type::String = "IncrementalInference.PackedLevelSetGridNormal"
  level::Float64
  sigma::Float64
  sigma_scale::Float64
  # make sure the JSON nested packing works with the serialization overlords
  heatmap::PackedHeatmapGridDensity
end


text:
module ThreePointBound

using ClusteredLowRankSolver, AbstractAlgebra, BasesAndSamples

export three_point_spherical_cap

function Q(n::Int, k::Int, u, v, t)
    R, x = PolynomialRing(RealField, "x")
    p = basis_gegenbauer(k, n, x)[end]
    sum(coeff(p, i) * ((1-u^2)*(1-v^2))^div(k-i, 2)*(t-u*v)^i for i=0:length(p)-1)
end

function m(w, d)
    [w^k for k=0:d]
end

p(u, costheta) = (u+1)*(costheta-u)

function three_point_spherical_cap(n, d, costheta, prec=256, all_free=false; basis=1, grid_option=2, N=2d, verbose=false, kwargs...)
    #basis: 1 = monomial
    #       2 = chebyshev
    #       3 = gegenbauer
    # These are orthogonalized with respect to the sample points, so it shouldn't really matter
    setprecision(BigFloat,prec)
    FF = RealField

    R, (u,v,t) = PolynomialRing(FF, ["u","v","t"])
    if verbose
        println("Creating the polynomials for F_k...")
    end
    F = Dict()
    for k=0:d
        #Instead of having both e.g. m(v,d-k) * m(u,d-k)^T and m(t,d-k) * m(u,d-k)^T, we combine them
        F[Block((:F,k))] = LowRankMatPol([R(1) for i=1:3],
                                    [m(u,d-k), m(v,d-k), m(t,d-k)], #combine the vectors which are used multiple times
                                    [1//6*(Q(n-1, k, u, v, t) .* m(v,d-k) + Q(n-1, k, u, t, v) .* m(t,d-k)),
                                    1//6*(Q(n-1, k, v, u, t) .* m(u,d-k) + Q(n-1, k, v, t, u) .* m(t,d-k)),
                                    1//6*(Q(n-1, k, t, v, u) .* m(v,d-k) + Q(n-1, k, t, u, v) .* m(u,d-k))])
    end
    #The vectors which construct the Pi_pi matrices. We have to tensor them with a symmetric basis for the E_pi matrices
    symmetry_weights = [[[R(1)]],
                        [[R((u-v)*(v-t)*(t-u))]],
                        [[1/sqrt(FF(2))*(2u-v-t),1/sqrt(FF(2))*(2v*t-u*t-u*v)],
                            [sqrt(FF(3)/FF(2))*(v-t),sqrt(FF(3)/FF(2))*(u*t-u*v)]]]
    #The polynomial weights which describe the domain
    weights = [R(1),
        p(u,costheta)+p(v,costheta)+p(t,costheta),
        p(u,costheta)*p(v,costheta)+p(v,costheta)*p(t,costheta)+p(t,costheta)*p(u,costheta),
        p(u,costheta)*p(v,costheta)*p(t,costheta),
        2u*v*t + 1 - u^2 - v^2 - t^2]
    #The basis in 1 variable used for the symmetric basis
    R1, (x,) = PolynomialRing(FF,["x"])
    if basis == 1
        basis1d_sym = m(x, 2d)
    elseif basis == 2
        basis1d_sym = basis_chebyshev(2d, x)
    elseif basis == 3
        basis1d_sym = basis_gegenbauer(2d, n, x)
    elseif basis == 4
        basis1d_sym = basis_laguerre(2d, BigFloat(n)/2-1,x)
        max_coef = [max(coefficients(basis1d_sym[i])...) for i = 1:length(basis1d_sym)]
        basis1d_sym = [max_coef[i]^(-1) * basis1d_sym[i] for i = 1:length(basis1d_sym)]
    else
        # also monomial basis for wrong input
        println("The option 'basis = $basis' was not recognized, we use the monomial basis as initial basis.")
        basis1d_sym = m(x, 2d)
    end
    if verbose
        println("Creating the polynomials for the SOS basis and the potential samples...")
    end
    # with the smart_eval basis, we don't compute the polynomials explicitely, but just the evaluations.
    sym_basis_smart_eval = [(u,v,t) -> basis1d_sym[deg-3k-2j+1](u+v+t) * basis1d_sym[j+1](u*v+v*t+u*t) * basis1d_sym[k+1](u*v*t) for deg=0:2d for k=0:div(deg,3) for j=0:div(deg-3k,2)]

    # we don't have a lot of double samples point due to symmetry, because we don't use exactly the same point sets for the x,y and z
    #We create a grid of sample points. either equally spaced or chebyshev-like
    if grid_option == 1
        samples = [-1 .+ (1+costheta) .* [i/N,j/N,k/N] for i=0:N for j=i:N for k=j:N] #this is kinda useless
    elseif grid_option == 2
        cheb_points = [vcat(sample_points_chebyshev(N+k,-1,costheta)...) for k=0:2]
        samples = [[cheb_points[1][i+1], cheb_points[2][j+1], cheb_points[3][k+1]] for i=0:N for j=0:N+1 for k=0:N+2]
    end
    samples = [
        s for s in samples if (1 + 2 * u * v * t - u^2 - v^2 - t^2)(s...) >= 0
    ] # We only keep the samples in the domain. This is not really needed
    #The degrees of the polynomials. They are ordered on degree, so we can look for the last one to determine the length of the bases we need later.
    degrees = [deg for deg=0:2d for k=0:div(deg,3) for j=0:div(deg-3k,2)]
    last_degree = [findlast(x -> x == k, degrees) for k=0:2d]

    if verbose
        println("Computing the approximate Fekete points and a corresponding basis...")
    end
    sym_basis, samples = approximatefekete(sym_basis_smart_eval,samples)

    if verbose
        println("Creating the SOS part of the trivariate constraint...")
    end
    for wi=1:length(weights)
        for swi=1:length(symmetry_weights)
            rank = length(symmetry_weights[swi])
            #this has in general too many entries, so we remove the ones with too high degree
            #kron(a,b) does [a[1]*b, a[2]*b,..., a[end] * b]. So the degree of vecs[r][i] is total_degree(symmetry_weights[swi][r][div(i,length(b))]+degree[(i-1)%length(b)+1]
            vecs = [kron(symmetry_weights[swi][r], sym_basis[1:last_degree[1+2d-total_degree(weights[wi])]]) for r=1:rank]
            for r=1:rank
                len = last_degree[1+2d-total_degree(weights[wi])] # the length of the symmetric basis we tensored the symmetric weight with
                #We keep the entries if the resulting diagonal element has degree at most 2d-deg(weight), since the weight will be set as the 'eigenvalue'. Note that this allows us to use sample points outside the domain
                keep_idx = [i for i=1:length(vecs[r])
                    if 2*(total_degree(symmetry_weights[swi][r][div(i-1,len)+1]) + degrees[(i-1)%len+1]) <= 2d-total_degree(weights[wi])]
                vecs[r] = vecs[r][keep_idx]
            end
            F[Block((:trivariatesos, wi,swi))] = LowRankMatPol([weights[wi] for r=1:rank], vecs)
        end
    end
    trivariatecon = Constraint(R(0), F, Dict(), samples)

    if verbose
        println("Creating the univariate constraint")
    end
    W, (w,) = PolynomialRing(FF, ["w"])
    f = Dict()
    for k=0:d # F(u,u,1)
        # from the (normally) three vu^T items, we only need 2 if we sum the ones with common u^T
        f[Block((:F,k))] = LowRankMatPol([W(1), W(1)],
                            [m(w,d-k), m(W(1),d-k)],
                            [Q(n-1,k,w,w,1) .* m(w,d-k)+Q(n-1,k,w,1,w) .* m(W(1),d-k),Q(n-1,k,1,w,w) .* m(w,d-k)])
    end
    basis = basis_gegenbauer(2d, n, w)
    for k=0:2d # a_k
        f[Block((:a,k))] = LowRankMatPol([basis[k+1]], [[W(1)]])
    end
    uni_basis = basis_chebyshev(2d, x)
    samples1d = sample_points_chebyshev(2d, -1, costheta) #we can take more samples, now we just get a better basis based on these samples. But since chebyshev points are relatively good, this is not a problem
    uni_basis,samples1d = approximatefekete(uni_basis, samples1d)

    f[Block((:uni_SOS,1))] = LowRankMatPol([W(1)], [uni_basis[1:d+1]])
    f[Block((:uni_SOS,2))] = LowRankMatPol([p(w,costheta)], [uni_basis[1:d]])

    univariatecon = Constraint(W(-1), f, Dict(), samples1d)

    obj = Dict()
    obj[Block((:F, 0))] = ones(d+1, d+1)
    for k=0:2d
        obj[Block((:a,k))] = ones(1, 1)
    end
    objective = Objective(1, obj, Dict())

    sos = LowRankPolProblem(false, objective, [univariatecon, trivariatecon])
    if verbose
        println("Converting to a clustered low-rank SDP...")
    end
    if all_free
        sdp = ClusteredLowRankSDP(sos,[(:F,k) for k=0:d])
    else
        sdp = ClusteredLowRankSDP(sos)
    end
    status, sol, time, errorcode = solvesdp(sdp;kwargs...)
end

end # of module


text:
#!/usr/bin/env julia

using MyFunctions: fact

for n in 0:5
    println("fact($n) = $(fact(n))")
end


text:
#!/usr/bin/julia
doc = """Usage:
    dpd_sim.py [--L <L> --N <N> --steps <ns> --dt <dt> --thermo <th> --gamma <g>]
               [--read <xyz> --config <config>]

Options:
    --N <N>            Number of atoms [default: 500]
    --L <L>            Box size [default: 5]
    --steps <ns>       Number of steps [default: 300]
    --dt <dt>          Timestep [default: 0.01]
    --thermo <th>      Frequency of printing to screen and file [default: 10]
    --gamma <g>        Friction [default: 0.0]
    --read <xyz>       Read the xyz file
    --config <config>  Read the DL_MESO CONFIG file

08/07/17
"""
using DocOpt
include("dpd_functions.jl")
include("dpd_io.jl")
using dpd_functions.init_pos, dpd_functions.init_vel,
    dpd_functions.integrate_verlet, dpd_functions.integrate_euler

args = docopt(doc)
#for (k, v) in args; println(k, "  ", v); end
L = parse(Float64, args["--L"])
N = parse(Int, args["--N"])
dt = parse(Float64, args["--dt"])
Ns = parse(Int, args["--steps"])
thermo = parse(Int, args["--thermo"])
gamma = parse(Float64, args["--gamma"])
kT = 1.0
rc = 1.0

rho = N / L^3
seed = 1234
srand(seed)

@printf "===== DPD simulations =====\n"
@printf "Beads: %i | Box: %.2f | rho: %.2f | gamma: %.1f\n" N L rho gamma
@printf "Steps: %i | dt: %.3f | thermo: %i\n" Ns dt thermo

A = 25.0 * ones(2, 2)

if args["--read"] != nothing
    fname = args["--read"]
    println("Reading initial xyz from $fname...")
    blist, X = read_xyz(fname)
    N = size(X, 1)
    V = init_vel(N, kT)
elseif args["--config"] != nothing
    fname = args["--config"]
    println("Reading initial positions and velocities from $fname...")
    X, V = read_config(fname)
    N = size(X, 1)
    blist = ones(Int, N)

else
    println("Random initialisation of positions and velocities.")
    blist = ones(Int, N)
    X = init_pos(N, L)
    V = init_vel(N, kT)
end

ti = time()
#T, KE, PE = integrate_euler(X, V, A, blist, L, gamma, kT, dt, rc, Ns, thermo)
T, KE, PE = integrate_verlet(X, V, A, blist, L, gamma, kT, dt, rc, Ns, thermo)

println("Simulation finished.Time: $(time() - ti) s.")
writedlm("data.out", [1:Ns+1 T KE PE])

max_stars_repo_name: petervanya/DPDsim

max_issues_repo_name: petervanya/DPDsim

max_forks_repo_name: petervanya/DPDsim


text:
# p8.m - eigenvalues of harmonic oscillator -u"+x^2 u on R

L = 8;                             # domain is [-L L], periodic
for N = 6:6:36
    h = 2*pi/N; x = h*(1:N); x = L*(x-pi)/pi;
    column = [-pi^2/(3*h^2)-1/6; @. -.5*(-1)^(1:N-1)/sin(h*(1:N-1)/2)^2 ];
    D2 = (pi/L)^2*toeplitz(column);  # 2nd-order differentiation
    eigenvalues = sort(eigvals(-D2 + diagm(x.^2)));
    @show N
    [ println(eigenvalues[i]) for i=1:4 ];
    println("");
end


text:
abstract type TA_Gaugefields_4D{NC} <: TA_Gaugefields{NC,4}
end

include("./TA_gaugefields_4D_serial.jl")

#=
function __init__()
    @require MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" begin
        include("./TA_gaugefields_4D_mpi.jl")
    end
end

=#

function TA_Gaugefields(NC,NX,NY,NZ,NT;mpi=false)
    if mpi
        return TA_Gaugefields_4D_mpi(NC,NX,NY,NZ,NT)
        #error("mpi = $mpi is not supoorted")
    else
        return TA_Gaugefields_4D_serial(NC,NX,NY,NZ,NT)
    end
end

function clear_U!(U::Array{T,1}) where T <: TA_Gaugefields_4D
    for μ=1:4
        clear_U!(U[μ])
    end
end


text:
using FileIO, JLD2

function addToFile(filename::String, varName::String, var::Any)
	if isfile(filename)
		data = load(filename)
		entry = joinpath("simulation", varName)
		if haskey(data, entry)
			delete!(data, entry)
		end
		data[entry] = var
		save(filename, data)
		#=jldopen(filename, "a+") do file
			entry = joinpath("simulation", varName)
			if haskey(file, entry)
				delete!(file, entry)
			end
			file[entry] = var
			#file["simulation"][varName] = var
		end=#
	else
		jldopen(filename, "w") do file
    		simulation = JLD2.Group(file, "simulation")
    		simulation[varName] = var
		end
	end
end

addToDict(dict::Dict, varName::String, var::Any) = merge(dict, Dict(varName => var))

loadFile(filename::String) = jldopen(filename, "r")["simulation"]

loadFromFile(filename::String, varName::String) = loadFile(filename)[varName]

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
using Pkg
Pkg.activate("./")
# load pkgs
using Ipopt, HarmonicPowerModels, PowerModels
using JuMP #avoids problems with Revise

# pkg const
const _PMs = PowerModels
const _HPM = HarmonicPowerModels

# path to the data
# path = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_basf_simplified_no_filter.m")
path = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_basf_simplified_with_filter.m")
# path = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_basf.m")

# transformer excitation data
exc_1 = Dict("voltage_harmonics" => [1,5],
            "current_harmonics" => [1,3,5,7,9,13],
            "N" => 50,
            "current_type" => :rectangular,
            "excitation_type" => :sigmoid,
            "inom" => 0.13,
            "ψmax" => 1,
            "voltage_type" => :rectangular,
            "dv" => [0.1,0.1],
            "vmin" => [-1.1,-1.1],
            "vmax" => [1.1,1.1],
            "dθ" => [π/5,π/5],
            "θmin" => [0.0,0.0],
            "θmax" => [2π,2π])

# exc_2 = Dict("voltage_harmonics" => [1,3],
#             "current_harmonics" => [1,3],
#             "N" => 50,
#             "current_type" => :rectangular,
#             "excitation_type" => :sigmoid,
#             "inom" => 0.4,
#             "ψmax" => 0.5,
#             "voltage_type" => :rectangular,
#             "dv" => [0.1,0.1],
#             "vmin" => [-1.1,-1.1],
#             "vmax" => [1.1,1.1],
#             "dθ" => [π/5,π/5],
#             "θmin" => [0.0,0.0],
#             "θmax" => [2π,2π])

# xfmr_exc = Dict("1" => exc_1, "2" => exc_2)

# load data
data  = _PMs.parse_file(path)

vmmin = 0.8
for (b,branch) in data["branch"]
    branch["c_rating"] = branch["rate_a"]/vmmin
end
for (d,load) in data["load"]
    load["c_rating"] = abs(load["pd"] + im* load["qd"])/vmmin
end
for (g,gen) in data["gen"]
    gen["c_rating"] = abs(gen["pmax"] + im* gen["qmax"])/vmmin
end

# hdata = _HPM.replicate(data, xfmr_exc=exc_1)
hdata = _HPM.replicate(data)

for (n, nw) in hdata["nw"]
    for (t, xfmr) in nw["xfmr"]
        if t != "1"
            # @show xfmr
            xfmr["voltage_harmonics"] = String[]
            xfmr["current_harmonics"] = String[]
            xfmr["voltage_harmonics_ntws"] = []
            xfmr["current_harmonics_ntws"] = []
        end
    end
end

# set the solver
solver = Ipopt.Optimizer

#solve power flow
# resultpf = run_hpf_iv(hdata, _PMs.IVRPowerModel, solver)
# @assert resultpf["termination_status"] == LOCALLY_SOLVED
# _HPM.append_indicators!(resultpf, hdata)

# println("Harmonic 13")
# _PMs.print_summary(resultpf["solution"]["nw"]["6"])
# println("Harmonic 9")
# _PMs.print_summary(resultpf["solution"]["nw"]["5"])
# println("Harmonic 7")
# _PMs.print_summary(resultpf["solution"]["nw"]["4"])
# println("Harmonic 5")
# _PMs.print_summary(resultpf["solution"]["nw"]["3"])
# vm = resultpf["solution"]["nw"]["3"]["bus"]["1"]

# println("Harmonic 3")
# _PMs.print_summary(resultpf["solution"]["nw"]["2"])
# vm = resultpf["solution"]["nw"]["2"]["bus"]["6"]
# println("Harmonic 1")
# _PMs.print_summary(resultpf["solution"]["nw"]["1"])

# pg = resultpf["solution"]["nw"]["1"]["gen"]["1"]

##
# solve the hopf
result = run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
pm = _PMs.instantiate_model(hdata, _PMs.IVRPowerModel, _HPM.build_hopf_iv; ref_extensions=[_HPM.ref_add_xfmr!]);
result = optimize_model!(pm, optimizer=solver, solution_processors=[ _HPM.sol_data_model!])
@assert result["termination_status"] == LOCALLY_SOLVED
_HPM.append_indicators!(result, hdata)

pg = result["solution"]["nw"]["1"]["gen"]["1"]
println("Harmonic 13")
_PMs.print_summary(result["solution"]["nw"]["6"])
println("Harmonic 9")
_PMs.print_summary(result["solution"]["nw"]["5"])
println("Harmonic 7")
_PMs.print_summary(result["solution"]["nw"]["4"])
println("Harmonic 5")
_PMs.print_summary(result["solution"]["nw"]["3"])
println("Harmonic 3")
_PMs.print_summary(result["solution"]["nw"]["2"])
println("Harmonic 1")
_PMs.print_summary(result["solution"]["nw"]["1"])
result["objective"]
result["termination_status"]

Dict(n=>(nw["gen"]["2"]["crg"]+ im*nw["gen"]["2"]["cig"] ) for (n,nw) in result["solution"]["nw"])
result["solution"]["nw"]["1"]["gen"]["1"]["pg"] + im* result["solution"]["nw"]["1"]["gen"]["1"]["qg"]
# @show result["solution"]["nw"]["1"]["xfmr"]["1"]["pexc"]


text:
# testtools.jl

function simulatedBEDMatrix(n::Integer, p::Integer, datatype::DataType=UInt8, navalue=BEDMatrices.NA_byte,
                            missingvalues=true, flips=false)
    srand(0)  # for reproducibility
    if missingvalues
        M = rand(0b00:0b11, n, p)
    else
        M = rand([0b00, 0b10, 0b11], n, p)
    end

    byten = cld(n, 4)
    numlastsnps = n - 4*(byten - 1)

    X = Matrix{UInt8}(byten, p)
    for col in 1:p
        for r in 1:(byten - 1)
            X[r, col] = quarterstobyte(M[(4r - 3):(4r), col])
        end
        X[byten, col] = quarterstobyte([q <= numlastsnps ? M[end-numlastsnps + q, col] : 0b0 for q in 1:4])
    end

    path = "[memory]"
    rownames = [string("row_", j) for j in 1:n]
    colnames = [string("col_", j) for j in 1:p]
    if flips
        flipvec = BitVector(rand(Bool, p))
    else
        flipvec = falses(p)
    end

    bed = BEDMatrix{datatype, typeof(X)}(n, p, X, convert(datatype, navalue),
                                         path, colnames, rownames, BEDMatrices.getbytemap(convert(datatype, navalue)), flipvec)
    navalue_typed = convert(datatype, navalue)

    if flips
        data = Matrix{datatype}(n, p)
        for col in 1:p
            if flipvec[col]
                for row in 1:n
                    x = M[row, col] === 0b01 ? navalue_typed : convert(datatype, BEDMatrices.rawformat(M[row, col]))
                    data[row, col] = x == navalue_typed ? x : 2 - x
                end
            else
                for row in 1:n
                    x = M[row, col] === 0b01 ? navalue_typed : convert(datatype, BEDMatrices.rawformat(M[row, col]))
                    data[row, col] = x
                end
            end
        end
    else
        data = map(x -> x === 0b01 ? navalue_typed : convert(datatype, BEDMatrices.rawformat(x)), M)
    end
    return (data, bed)
end

function quarterstobyte(v::Vector{UInt8})
    byte = 0x0
    for (qidx, q) in enumerate(v)
        byte += q << (2*(qidx - 1))
    end
    byte
end

function readRAW(filename::String)
    n = countlines(filename) - 1

    lines = readlines(filename)

    header = split(lines[1])
    p = length(header) - 6
    colnames = header[7:end]
    data = Matrix{UInt8}(n, p)
    rownames = Vector{String}(n)

    for (row, line) in enumerate(lines[2:end])
        fields = split(line)
        rownames[row] = fields[1]*"_"*fields[2]
        data[row, :] = map(fld -> fld == "NA" ? BEDMatrices.NA_byte : parse(UInt8, fld), fields[7:end])
    end

    return rownames, colnames, data
end


text:
@testset "qraddcol with β = 0" begin
    m = 100
    A = randn(m,0)
    R = zeros(0,0)

    for i in 1:m
        a = randn(m)
        R = qraddcol(A, R, a)
        A = [A a]
        @test  norm( R'*R - A'*A ) <1e-5
    end
end

@testset "qraddcol with β > 0" begin
    m = 100
    A = randn(m,0)
    R = zeros(0,0)
    β = 0.1
    for i in 1:m
        a = randn(m)
        B = [A a]
        R = qraddcol(A, R, a, β)
        @test norm( R'*R - B'*B - β^2*I ) <1e-5
    end
end

@testset "qraddrow" begin
    m, n = 3, 3
    A = randn(m,m)
    Q, R = qr(A)
    for i in 1:100
        a = randn(m)'
        A = [A; a]
        R = qraddrow(R, a)
        @test norm( R'R - A'*A ) <1e-5
    end

end

@testset "qrdelcol" begin
    m = 100
    A = randn(m,m)
    Q, R = qr(A)
    for i in 100:-1:1
        k = rand(1:i)
        A = A[:,1:i .!= k]
        R = qrdelcol(R, k)
        @test norm( R'*R - A'*A ) <1e-5#--> less_than(1e-5)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SamplingMethod

A method for sampling from geometries.
"""
abstract type SamplingMethod end

"""
    sample(geometry, method)

Sample elements from `geometry` with `method`.
"""
function sample end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("sampling/regular.jl")
include("sampling/uniform.jl")
include("sampling/weighted.jl")
include("sampling/ball.jl")

# ----------
# UTILITIES
# ----------

"""
    sample(object, nsamples, [weights], replace=false)

Generate `nsamples` samples from spatial `object`
uniformly or using `weights`, with or without
replacement depending on `replace` option.
"""
function sample(object::Union{Domain,Data}, nsamples::Int,
                weights::AbstractVector=[]; replace=false)
  if isempty(weights)
    sample(object, UniformSampling(nsamples, replace))
  else
    sample(object, WeightedSampling(nsamples, weights, replace))
  end
end


text:
using Dtree
using Base.Test

if dt_nodeid == 1
    println("Dtree: $dt_nnodes nodes")
end


text:
function _ensure_not_running_on_functor(f, name)
    # if x itself is a Type, then it is a constructor, thus not a functor.
    # This also catchs UnionAll constructors which have a `:var` and `:body` fields
    f isa Type && return

    if fieldcount(typeof(f)) > 0
        throw(ArgumentError(
            "$name cannot be used on closures/functors (such as $f)"
        ))
    end
end

"""
    _make_fdm_call(fdm, f, ȳ, xs, ignores) -> Tuple

Call `FiniteDifferences.j′vp`, with the option to ignore certain `xs`.

# Arguments
- `fdm::FiniteDifferenceMethod`: How to numerically differentiate `f`.
- `f`: The function to differentiate.
- `ȳ`: adjoint w.r.t. output of `f`.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ignores`: Collection of `Bool`s, the same length as `xs`.
  If `ignores[i] === true`, then `xs[i]` is ignored; `∂xs[i] === nothing`.

# Returns
- `∂xs::Tuple`: Derivatives estimated by finite differencing.
"""
function _make_fdm_call(fdm, f, ȳ, xs, ignores)
    sig = Expr(:tuple)
    call = Expr(:call, f)
    newxs = Any[]
    arginds = Int[]
    i = 1
    for (x, ignore) in zip(xs, ignores)
        if ignore
            push!(call.args, x)
        else
            push!(call.args, Symbol(:x, i))
            push!(sig.args, Symbol(:x, i))
            push!(newxs, x)
            push!(arginds, i)
        end
        i += 1
    end
    fdexpr = :(j′vp($fdm, $sig -> $call, $ȳ, $(newxs...)))
    fd = eval(fdexpr)
    args = Any[nothing for _ in 1:length(xs)]
    for (dx, ind) in zip(fd, arginds)
        args[ind] = dx
    end
    return (args...,)
end

"""
    test_scalar(f, x; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

Given a function `f` with scalar input and scalar output, perform finite differencing checks,
at input point `x` to confirm that there are correct `frule` and `rrule`s provided.

# Arguments
- `f`: Function for which the `frule` and `rrule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).

All keyword arguments except for `fdm` is passed to `isapprox`.
"""
function test_scalar(f, x; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    _ensure_not_running_on_functor(f, "test_scalar")

    r_res = rrule(f, x)
    f_res = frule((Zero(), 1), f, x)
    @test r_res !== nothing  # Check the rule was defined
    @test f_res !== nothing
    r_fx, prop_rule = r_res
    f_fx, f_∂x = f_res
    @testset "$f at $x, $(nameof(rule))" for (rule, fx, ∂x) in (
        (rrule, r_fx, prop_rule(1)),
        (frule, f_fx, f_∂x)
    )
        @test fx == f(x)  # Check we still get the normal value, right

        if rule == rrule
            ∂self, ∂x = ∂x
            @test ∂self === NO_FIELDS
        end
        @test isapprox(∂x, fdm(f, x); rtol=rtol, atol=atol, kwargs...)
    end
end

"""
    frule_test(f, (x, ẋ)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function for which the `frule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `ẋ`: differential w.r.t. `x` (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function frule_test(f, (x, ẋ); rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    return frule_test(f, ((x, ẋ),); rtol=rtol, atol=atol, fdm=fdm, kwargs...)
end

function frule_test(f, xẋs::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    _ensure_not_running_on_functor(f, "frule_test")
    xs, ẋs = first.(xẋs), last.(xẋs)
    Ω, dΩ_ad = frule((NO_FIELDS, ẋs...), f, xs...)
    @test f(xs...) == Ω

    # Correctness testing via finite differencing.
    dΩ_fd = jvp(fdm, xs->f(xs...), (xs, ẋs))
    @test isapprox(
        collect(extern.(dΩ_ad)),  # Use collect so can use vector equality
        collect(dΩ_fd);
        rtol=rtol,
        atol=atol,
        kwargs...
    )
end

"""
    rrule_test(f, ȳ, (x, x̄)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function to which rule should be applied.
- `ȳ`: adjoint w.r.t. output of `f` (should generally be set randomly).
  Should be same structure as `f(x)` (so if multiple returns should be a tuple)
- `x`: input at which to evaluate `f` (should generally be set to an arbitary point in the domain).
- `x̄`: currently accumulated adjoint (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function rrule_test(f, ȳ, (x, x̄)::Tuple{Any, Any}; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    _ensure_not_running_on_functor(f, "rrule_test")

    # Check correctness of evaluation.
    fx, pullback = rrule(f, x)
    @test collect(fx) ≈ collect(f(x))  # use collect so can do vector equality
    (∂self, x̄_ad) = if fx isa Tuple
        # If the function returned multiple values,
        # then it must have multiple seeds for propagating backwards
        pullback(ȳ...)
    else
        pullback(ȳ)
    end

    @test ∂self === NO_FIELDS  # No internal fields
    # Correctness testing via finite differencing.
    x̄_fd = only(j′vp(fdm, f, ȳ, x))  # j′vp returns a tuple, but `f` is a unary function.
    @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)
end

# case where `f` takes multiple arguments
function rrule_test(f, ȳ, xx̄s::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    _ensure_not_running_on_functor(f, "rrule_test")

    # Check correctness of evaluation.
    xs, x̄s = collect(zip(xx̄s...))
    y, pullback = rrule(f, xs...)
    @test f(xs...) == y

    @assert !(isa(ȳ, Thunk))
    ∂s = pullback(ȳ)
    ∂self = ∂s[1]
    x̄s_ad = ∂s[2:end]
    @test ∂self === NO_FIELDS

    # Correctness testing via finite differencing.
    x̄s_fd = _make_fdm_call(fdm, f, ȳ, xs, x̄s .== nothing)
    for (x̄_ad, x̄_fd) in zip(x̄s_ad, x̄s_fd)
        if x̄_fd === nothing
            # The way we've structured the above, this tests the propagator is returning a DoesNotExist
            @test x̄_ad isa DoesNotExist
        else
            @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)
        end
    end
end


text:
using Descartes
using GeometryBasics: Mesh

# params
beam_size = [50,10,10]
hole_ct = 5
hole_d = 3

# computations
hole_interval = beam_size[1]/(hole_ct + 1)

c = Cuboid(beam_size)

for i = 1:hole_ct
    h = translate([hole_interval*i, -1, beam_size[3]/2])*
            rotate(-pi/2, [1,0,0])*
                Cylinder(hole_d/2, beam_size[2]+2, center=false)
    global c = diff(c, h)
end

save("fea.ply",Mesh(c))


text:
# Tests nothing currently
# Maybe use FactCheck.jl

using IsingLite
using Base.Test

@test 1 == 1


text:
mutable struct scratch_data
    n_particles::Int
    polygon::Vector{Point2D{Float64}}
    extremePoint::Point2D{Float64}
    pointOnWall::Point2D{Float64}
    w::Vector{Float64}
    uEl::Vector{Float64}
    vEl::Vector{Float64}
    a::Vector{Float64}
    b::Vector{Float64}
    c::Vector{Float64}
    d::Vector{Float64}
    cfX::Vector{Float64}
    cfY::Vector{Float64}
    afX::Vector{Float64}
    afY::Vector{Float64}
    gfX::Vector{Float64}
    gfY::Vector{Float64}
    sfX::Vector{Float64}
    sfY::Vector{Float64}

    function scratch_data(n_particles::Int)
        polygon      = [Point2D(0.0,0.0),Point2D(0.0,0.0),Point2D(0.0,0.0),Point2D(0.0,0.0)]
        extremePoint = Point2D(100_000.0,0.0)
        pointOnWall  = Point2D(0.0,0.0)
        w   = zeros(4)
        uEl = zeros(4)
        vEl = zeros(4)
        a   = zeros(3)
        b   = zeros(3)
        c   = zeros(3)
        d   = zeros(3)
        cfX = zeros(n_particles)
        cfY = zeros(n_particles)
        afX = zeros(n_particles)
        afY = zeros(n_particles)
        gfX = zeros(n_particles)
        gfY = zeros(n_particles)
        sfX = zeros(n_particles)
        sfY = zeros(n_particles)
        return new(n_particles,polygon,extremePoint,pointOnWall,
                    w,uEl,vEl,a,b,c,d,
                    cfX,cfY,afX,afY,gfX,gfY,sfX,sfY)
    end
end

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart