text:
"""
    manifold_element_forwards(T, field::Symbol)
    manifold_element_forwards(T, Twhere, field::Symbol)

Introduce basic fallbacks for type `T` (which can be a subtype of `Twhere`) that represents
points or vectors for a manifold.
Fallbacks will work by forwarding to the field passed in `field``

List of forwarded functions:
* [`allocate`](@ref),
* [`copy`](@ref),
* [`copyto!`](@ref),
* [`number_eltype`](@ref) (only for values, not the type itself),
* `similar`,
* `size`,
* `==`.
"""
macro manifold_element_forwards(T, field::Symbol)
    return esc(quote
        ManifoldsBase.@manifold_element_forwards ($T) _ ($field)
    end)
end
macro manifold_element_forwards(T, Twhere, field::Symbol)
    return esc(
        quote
            ManifoldsBase.allocate(p::$T) where {$Twhere} = $T(allocate(p.$field))
            function ManifoldsBase.allocate(p::$T, ::Type{P}) where {P,$Twhere}
                return $T(allocate(p.$field, P))
            end
            function ManifoldsBase.allocate(
                p::$T,
                ::Type{P},
                dims::Tuple,
            ) where {P,$Twhere}
                return $T(allocate(p.$field, P, dims))
            end

            @inline Base.copy(p::$T) where {$Twhere} = $T(copy(p.$field))

            function Base.copyto!(q::$T, p::$T) where {$Twhere}
                copyto!(q.$field, p.$field)
                return q
            end

            function ManifoldsBase.number_eltype(p::$T) where {$Twhere}
                return typeof(one(eltype(p.$field)))
            end

            Base.similar(p::$T) where {$Twhere} = $T(similar(p.$field))
            Base.similar(p::$T, ::Type{P}) where {P,$Twhere} = $T(similar(p.$field, P))

            Base.size(p::$T) where {$Twhere} = size(p.$field)

            Base.:(==)(p::$T, q::$T) where {$Twhere} = (p.$field == q.$field)
        end,
    )
end

"""
    default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol)

Introduce default fallbacks for all basic functions on manifolds, for manifold of type `TM`,
points of type `TP`, tangent vectors of type `TV`, with forwarding to fields `pfield` and
`vfield` for point and tangent vector functions, respectively.
"""
macro default_manifold_fallbacks(TM, TP, TV, pfield::Symbol, vfield::Symbol)
    block = quote
        function ManifoldsBase.allocate_result(::$TM, ::typeof(log), p::$TP, ::$TP)
            a = allocate(p.$vfield)
            return $TV(a)
        end
        function ManifoldsBase.allocate_result(
            ::$TM,
            ::typeof(inverse_retract),
            p::$TP,
            ::$TP,
        )
            a = allocate(p.$vfield)
            return $TV(a)
        end
        function ManifoldsBase.allocate_coordinates(M::$TM, p::$TP, T, n::Int)
            return ManifoldsBase.allocate_coordinates(M, p.$pfield, T, n)
        end

        function ManifoldsBase.angle(M::$TM, p::$TP, X::$TV, Y::$TV)
            return angle(M, p.$pfield, X.$vfield, Y.$vfield)
        end

        function ManifoldsBase.check_point(M::$TM, p::$TP; kwargs...)
            return check_point(M, p.$pfield; kwargs...)
        end

        function ManifoldsBase.check_vector(M::$TM, p::$TP, X::$TV; kwargs...)
            return ManifoldsBase.check_vector(M, p.$pfield, X.$vfield; kwargs...)
        end

        function ManifoldsBase.distance(M::$TM, p::$TP, q::$TP)
            return distance(M, p.$pfield, q.$pfield)
        end

        function ManifoldsBase.embed!(M::$TM, q::$TP, p::$TP)
            embed!(M, q.$pfield, p.$pfield)
            return q
        end

        function ManifoldsBase.embed!(M::$TM, Y::$TV, p::$TP, X::$TV)
            embed!(M, Y.$vfield, p.$pfield, X.$vfield)
            return Y
        end

        function ManifoldsBase.exp!(M::$TM, q::$TP, p::$TP, X::$TV)
            exp!(M, q.$pfield, p.$pfield, X.$vfield)
            return q
        end

        function ManifoldsBase.inner(M::$TM, p::$TP, X::$TV, Y::$TV)
            return inner(M, p.$pfield, X.$vfield, Y.$vfield)
        end

        function ManifoldsBase.inverse_retract!(
            M::$TM,
            X::$TV,
            p::$TP,
            q::$TP,
            m::LogarithmicInverseRetraction,
        )
            inverse_retract!(M, X.$vfield, p.$pfield, q.$pfield, m)
            return X
        end

        function ManifoldsBase.isapprox(M::$TM, p::$TP, q::$TP; kwargs...)
            return isapprox(M, p.$pfield, q.$pfield; kwargs...)
        end

        function ManifoldsBase.isapprox(M::$TM, p::$TP, X::$TV, Y::$TV; kwargs...)
            return isapprox(M, p.$pfield, X.$vfield, Y.$vfield; kwargs...)
        end

        function ManifoldsBase.log!(M::$TM, X::$TV, p::$TP, q::$TP)
            log!(M, X.$vfield, p.$pfield, q.$pfield)
            return X
        end

        function ManifoldsBase.norm(M::$TM, p::$TP, X::$TV)
            return norm(M, p.$pfield, X.$vfield)
        end

        function ManifoldsBase.retract!(
            M::$TM,
            q::$TP,
            p::$TP,
            X::$TV,
            m::ExponentialRetraction,
        )
            retract!(M, q.$pfield, p.$pfield, X.$vfield, m)
            return X
        end

        function ManifoldsBase.vector_transport_along!(
            M::$TM,
            Y::$TV,
            p::$TP,
            X::$TV,
            c::AbstractVector,
        )
            vector_transport_along!(M, Y.$vfield, p.$pfield, X.$vfield, c)
            return Y
        end

        function ManifoldsBase.zero_vector(M::$TM, p::$TP)
            return $TV(zero_vector(M, p.$pfield))
        end

        function ManifoldsBase.zero_vector!(M::$TM, X::$TV, p::$TP)
            zero_vector!(M, X.$vfield, p.$pfield)
            return X
        end
    end
    for f_postfix in [:default, :orthogonal, :orthonormal, :vee, :cached, :diagonalizing]
        ca = Symbol("get_coordinates_$(f_postfix)")
        cm = Symbol("get_coordinates_$(f_postfix)!")
        va = Symbol("get_vector_$(f_postfix)")
        vm = Symbol("get_vector_$(f_postfix)!")
        B_types = if f_postfix in [:default, :orthogonal, :orthonormal, :vee]
            [:AbstractNumbers, :RealNumbers, :ComplexNumbers]
        elseif f_postfix === :cached
            [:CachedBasis]
        elseif f_postfix === :diagonalizing
            [:DiagonalizingOrthonormalBasis]
        else
            [:Any]
        end
        for B_type in B_types
            push!(
                block.args,
                quote
                    function ManifoldsBase.$ca(M::$TM, p::$TP, X::$TV, B::$B_type)
                        return ManifoldsBase.$ca(M, p.$pfield, X.$vfield, B)
                    end
                    function ManifoldsBase.$cm(M::$TM, Y, p::$TP, X::$TV, B::$B_type)
                        ManifoldsBase.$cm(M, Y, p.$pfield, X.$vfield, B)
                        return Y
                    end
                    function ManifoldsBase.$va(M::$TM, p::$TP, X, B::$B_type)
                        return $TV(ManifoldsBase.$va(M, p.$pfield, X, B))
                    end
                    function ManifoldsBase.$vm(M::$TM, Y::$TV, p::$TP, X, B::$B_type)
                        ManifoldsBase.$vm(M, Y.$vfield, p.$pfield, X, B)
                        return Y
                    end
                end,
            )
        end
    end
    # TODO  forward retraction / inverse_retraction
    for f_postfix in [:polar, :project, :qr, :softmax]
        ra = Symbol("retract_$(f_postfix)")
        rm = Symbol("retract_$(f_postfix)!")
        push!(block.args, quote
            function ManifoldsBase.$ra(M::$TM, p::$TP, X::$TV)
                return $TP(ManifoldsBase.$ra(M, p.$pfield, X.$vfield))
            end
            function ManifoldsBase.$rm(M::$TM, q, p::$TP, X::$TV)
                ManifoldsBase.$rm(M, q.$pfield, p.$pfield, X.$vfield)
                return q
            end
        end)
    end
    push!(
        block.args,
        quote
            function ManifoldsBase.retract_exp_ode(
                M::$TM,
                p::$TP,
                X::$TV,
                m::AbstractRetractionMethod,
                B::ManifoldsBase.AbstractBasis,
            )
                return $TP(ManifoldsBase.retract_exp_ode(M, p.$pfield, X.$vfield, m, B))
            end
            function ManifoldsBase.retract_exp_ode!(
                M::$TM,
                q::$TP,
                p::$TP,
                X::$TV,
                m::AbstractRetractionMethod,
                B::ManifoldsBase.AbstractBasis,
            )
                ManifoldsBase.retract_exp_ode!(M, q.$pfield, p.$pfield, X.$vfield, m, B)
                return q
            end
            function ManifoldsBase.retract_pade(M::$TM, p::$TP, X::$TV, n)
                return $TP(ManifoldsBase.retract_pade(M, p.$pfield, X.$vfield, n))
            end
            function ManifoldsBase.retract_pade!(M::$TM, q::$TP, p::$TP, X::$TV, n)
                ManifoldsBase.retract_pade!(M, q.$pfield, p.$pfield, X.$vfield, n)
                return q
            end
            function ManifoldsBase.retract_embedded(
                M::$TM,
                p::$TP,
                X::$TV,
                m::AbstractRetractionMethod,
            )
                return $TP(ManifoldsBase.retract_embedded(M, p.$pfield, X.$vfield, m))
            end
            function ManifoldsBase.retract_embedded!(
                M::$TM,
                q::$TP,
                p::$TP,
                X::$TV,
                m::AbstractRetractionMethod,
            )
                ManifoldsBase.retract_embedded!(M, q.$pfield, p.$pfield, X.$vfield, m)
                return q
            end
        end,
    )
    for f_postfix in [:polar, :project, :qr, :softmax]
        ra = Symbol("inverse_retract_$(f_postfix)")
        rm = Symbol("inverse_retract_$(f_postfix)!")
        push!(block.args, quote
            function ManifoldsBase.$ra(M::$TM, p::$TP, q::$TP)
                return $TV((ManifoldsBase.$ra)(M, p.$pfield, q.$pfield))
            end
            function ManifoldsBase.$rm(M::$TM, Y::$TV, p::$TP, q::$TP)
                ManifoldsBase.$rm(M, Y.$vfield, p.$pfield, q.$pfield)
                return Y
            end
        end)
    end
    push!(
        block.args,
        quote
            function ManifoldsBase.inverse_retract_embedded(
                M::$TM,
                p::$TP,
                q::$TP,
                m::AbstractInverseRetractionMethod,
            )
                return $TV(
                    ManifoldsBase.inverse_retract_embedded(M, p.$pfield, q.$pfield, m),
                )
            end
            function ManifoldsBase.inverse_retract_embedded!(
                M::$TM,
                X::$TV,
                p::$TP,
                q::$TP,
                m::AbstractInverseRetractionMethod,
            )
                ManifoldsBase.inverse_retract_embedded!(
                    M,
                    X.$vfield,
                    p.$pfield,
                    q.$pfield,
                    m,
                )
                return X
            end
            function ManifoldsBase.inverse_retract_nlsolve(
                M::$TM,
                p::$TP,
                q::$TP,
                m::NLSolveInverseRetraction,
            )
                return $TV(
                    ManifoldsBase.inverse_retract_nlsolve(M, p.$pfield, q.$pfield, m),
                )
            end
            function ManifoldsBase.inverse_retract_nlsolve!(
                M::$TM,
                X::$TV,
                p::$TP,
                q::$TP,
                m::NLSolveInverseRetraction,
            )
                ManifoldsBase.inverse_retract_nlsolve!(
                    M,
                    X.$vfield,
                    p.$pfield,
                    q.$pfield,
                    m,
                )
                return X
            end
        end,
    )
    # forward vector transports

    for sub in [:project, :diff]
        # project & diff
        vtaa = Symbol("vector_transport_along_$(sub)")
        vtam = Symbol("vector_transport_along_$(sub)!")
        vtta = Symbol("vector_transport_to_$(sub)")
        vttm = Symbol("vector_transport_to_$(sub)!")
        push!(
            block.args,
            quote
                function ManifoldsBase.$vtaa(M::$TM, p::$TP, X::$TV, c::AbstractVector)
                    return $TV(ManifoldsBase.$vtaa(M, p.$pfield, X.$vfield, c))
                end
                function ManifoldsBase.$vtam(
                    M::$TM,
                    Y::$TV,
                    p::$TP,
                    X::$TV,
                    c::AbstractVector,
                )
                    ManifoldsBase.$vtam(M, Y.$vfield, p.$pfield, X.$vfield, c)
                    return Y
                end
                function ManifoldsBase.$vtta(M::$TM, p::$TP, X::$TV, q::$TP)
                    return $TV(ManifoldsBase.$vtta(M, p.$pfield, X.$vfield, q.$pfield))
                end
                function ManifoldsBase.$vttm(M::$TM, Y::$TV, p::$TP, X::$TV, q::$TP)
                    ManifoldsBase.$vttm(M, Y.$vfield, p.$pfield, X.$vfield, q.$pfield)
                    return Y
                end
            end,
        )
    end
    # parallel transports
    push!(
        block.args,
        quote
            function ManifoldsBase.parallel_transport_along(
                M::$TM,
                p::$TP,
                X::$TV,
                c::AbstractVector,
            )
                return $TV(
                    ManifoldsBase.parallel_transport_along(M, p.$pfield, X.$vfield, c),
                )
            end
            function ManifoldsBase.parallel_transport_along!(
                M::$TM,
                Y::$TV,
                p::$TP,
                X::$TV,
                c::AbstractVector,
            )
                ManifoldsBase.parallel_transport_along!(
                    M,
                    Y.$vfield,
                    p.$pfield,
                    X.$vfield,
                    c,
                )
                return Y
            end
            function ManifoldsBase.parallel_transport_direction(
                M::$TM,
                p::$TP,
                X::$TV,
                d::$TV,
            )
                return $TV(
                    ManifoldsBase.parallel_transport_direction(
                        M,
                        p.$pfield,
                        X.$vfield,
                        d.$vfield,
                    ),
                )
            end
            function ManifoldsBase.parallel_transport_direction!(
                M::$TM,
                Y::$TV,
                p::$TP,
                X::$TV,
                d::$TV,
            )
                ManifoldsBase.parallel_transport_direction!(
                    M,
                    Y.$vfield,
                    p.$pfield,
                    X.$vfield,
                    d.$vfield,
                )
                return Y
            end
            function ManifoldsBase.parallel_transport_to(M::$TM, p::$TP, X::$TV, q::$TP)
                return $TV(
                    ManifoldsBase.parallel_transport_to(M, p.$pfield, X.$vfield, q.$pfield),
                )
            end
            function ManifoldsBase.parallel_transport_to!(
                M::$TM,
                Y::$TV,
                p::$TP,
                X::$TV,
                q::$TP,
            )
                ManifoldsBase.parallel_transport_to!(
                    M,
                    Y.$vfield,
                    p.$pfield,
                    X.$vfield,
                    q.$pfield,
                )
                return Y
            end
        end,
    )
    return esc(block)
end

@doc raw"""
    manifold_vector_forwards(T, field::Symbol)
    manifold_vector_forwards(T, Twhere, field::Symbol)

Introduce basic fallbacks for type `T` that represents vectors from a vector bundle for a
manifold. `Twhere` is put into `where` clause of each method. Fallbacks work by forwarding
to field passed as `field`.

List of forwarded functions:
* basic arithmetic (`*`, `/`, `\`, `+`, `-`),
* all things from [`@manifold_element_forwards`](@ref),
* broadcasting support.

# example

    @eval @manifold_vector_forwards ValidationFibreVector{TType} TType value
"""
macro manifold_vector_forwards(T, field::Symbol)
    return esc(quote
        ManifoldsBase.@manifold_vector_forwards ($T) _ ($field)
    end)
end
macro manifold_vector_forwards(T, Twhere, field::Symbol)
    return esc(
        quote
            Base.:*(X::$T, s::Number) where {$Twhere} = $T(X.$field * s)
            Base.:*(s::Number, X::$T) where {$Twhere} = $T(s * X.$field)
            Base.:/(X::$T, s::Number) where {$Twhere} = $T(X.$field / s)
            Base.:\(s::Number, X::$T) where {$Twhere} = $T(s \ X.$field)
            Base.:+(X::$T, Y::$T) where {$Twhere} = $T(X.$field + Y.$field)
            Base.:-(X::$T, Y::$T) where {$Twhere} = $T(X.$field - Y.$field)
            Base.:-(X::$T) where {$Twhere} = $T(-X.$field)
            Base.:+(X::$T) where {$Twhere} = $T(X.$field)
            Base.zero(X::$T) where {$Twhere} = $T(zero(X.$field))

            @eval ManifoldsBase.@manifold_element_forwards $T $Twhere $field

            Base.axes(p::$T) where {$Twhere} = axes(p.$field)

            function Broadcast.BroadcastStyle(::Type{<:$T}) where {$Twhere}
                return Broadcast.Style{$T}()
            end
            function Broadcast.BroadcastStyle(
                ::Broadcast.AbstractArrayStyle{0},
                b::Broadcast.Style{$T},
            ) where {$Twhere}
                return b
            end

            function Broadcast.instantiate(
                bc::Broadcast.Broadcasted{Broadcast.Style{$T},Nothing},
            ) where {$Twhere}
                return bc
            end
            function Broadcast.instantiate(
                bc::Broadcast.Broadcasted{Broadcast.Style{$T}},
            ) where {$Twhere}
                Broadcast.check_broadcast_axes(bc.axes, bc.args...)
                return bc
            end

            Broadcast.broadcastable(X::$T) where {$Twhere} = X

            @inline function Base.copy(
                bc::Broadcast.Broadcasted{Broadcast.Style{$T}},
            ) where {$Twhere}
                return $T(Broadcast._broadcast_getindex(bc, 1))
            end

            Base.@propagate_inbounds function Broadcast._broadcast_getindex(
                X::$T,
                I,
            ) where {$Twhere}
                return X.$field
            end

            @inline function Base.copyto!(
                dest::$T,
                bc::Broadcast.Broadcasted{Broadcast.Style{$T}},
            ) where {$Twhere}
                axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))
                # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
                if bc.f === identity && bc.args isa Tuple{$T} # only a single input argument to broadcast!
                    A = bc.args[1]
                    if axes(dest) == axes(A)
                        return copyto!(dest, A)
                    end
                end
                bc′ = Broadcast.preprocess(dest, bc)
                # Performance may vary depending on whether `@inbounds` is placed outside the
                # for loop or not. (cf. https://github.com/JuliaLang/julia/issues/38086)
                copyto!(dest.$field, bc′[1])
                return dest
            end
        end,
    )
end


text:
mutable struct Waypoint
    x::Int
    y::Int
end

struct Coordinate
    x::Int
    y::Int
end

struct Action
    instruction::Char
    amount::Int
end

function getRoute(file)
    actions::Vector{Action} = []
    for line in readlines(file)
        push!(actions, Action(line[1], parse(Int,line[2:length(line)])))
    end
    #println(actions)
    return actions
end

function rotate(wp, degrees)
    x::Int = wp.x
    y::Int = wp.y
    clockwise = degrees > 0
    quadrants = degrees ÷ 90
    newWp = Waypoint(x,y)
    if (abs(quadrants) == 1)
        newWp.x = y
        newWp.y = x
        if clockwise
            newWp.y *= -1
        else
            newWp.x *= -1
        end
    end
    if (abs(quadrants) == 2)
        newWp.x *= -1
        newWp.y *= -1
    end
    if(abs(quadrants) == 3)
        newWp.x = y
        newWp.y = x
        if clockwise
            newWp.x *= -1
        else
            newWp.y *= -1
        end
    end
    return newWp
end

# function testRotate()
#     a = Waypoint(2,1)
#     println("a90", rotate(a,90))
#     println("a180", rotate(a,180))
#     println("a270", rotate(a,270))
#     println("a-90", rotate(a,-90))
#     println("a-180", rotate(a,-180))
#     println("a-270", rotate(a,-270))
# end
# testRotate()

function moveWaypoint(currentWaypoint, action)
    x::Int = currentWaypoint.x
    y::Int = currentWaypoint.y
    newWp = Waypoint(x,y)
    if action.instruction == 'N'
        newWp.y += action.amount
    end
    if action.instruction == 'E'
        newWp.x += action.amount
    end
    if action.instruction == 'S'
        newWp.y -= action.amount
    end
    if action.instruction == 'W'
        newWp.x -= action.amount
    end
    if action.instruction == 'L'
        degrees = -1 * action.amount
        newWp = rotate(currentWaypoint, degrees)
    end
    if action.instruction == 'R'
        degrees = action.amount
        newWp = rotate(currentWaypoint, degrees)
    end
    return newWp
end

function moveShip(currPos, wp, amount)
    x::Int = currPos.x
    y::Int = currPos.y
    newX = x + wp.x * amount
    newY = y + wp.y * amount
    newPos = Coordinate(newX, newY)
    return newPos
end

function moveWaypointWithShip(currPos, currWp)
    x::Int = currPos.x
    y::Int = currPos.y
    newWaypoint = Waypoint(x + currWp.x, y + currWp.y)
    return newWaypoint
end

function navigate(route)
    currPosition = Coordinate(0,0)
    currWaypoint = Waypoint(10,1)
    for action in route
        println(currPosition, " ", currWaypoint, " ", action)
        if action.instruction == 'F'
            currPosition = moveShip(currPosition, currWaypoint, action.amount)
        else
            currWaypoint = moveWaypoint(currWaypoint, action)
        end
    end
    manhattanDistance = abs(currPosition.x) + abs(currPosition.y)
    println("Manhattan Distance = ", manhattanDistance)
end

route = getRoute("data.txt")
navigate(route)

# function testMoveShip()
#     cp = Coordinate(0,0)
#     cw = Waypoint(10,-1)
#     newPos = moveShip(cp, cw, 2)
#     println(newPos)
# end
# testMoveShip()

ext: jl


text:
import FFTW: fft!, ifft!
import LinearAlgebra: transpose

export Ampere

struct Ampere end

"""

    advection!( fᵀ, mesh1, mesh2, E, dt, type, axis )

- axis == 1 Advection in x and compute electric field

    ∂f/∂t − υ∂f/∂x  = 0
    ∂E/∂t = −J = ∫ fυ dυ

- axis == 2 Advection in υ

    ∂f/∂t − E(x) ∂f/∂υ  = 0

"""
function advection!( f     :: Array{Complex{Float64},2},
                     fᵀ    :: Array{Complex{Float64},2},
		             mesh1 :: UniformMesh,
		             mesh2 :: UniformMesh,
		             e     :: Vector{Complex{Float64}},
		             dt    :: Float64,
                     atype :: Ampere,
		             axis  :: Int64 )

    @assert ( axis == 1 || axis == 2 )

    if ( axis == 2 )

        n = mesh2.length ÷ 2
        L = mesh2.stop - mesh2.start
        k = 2π/L * [0:1:n-1;-n:1:-1]
        ek = exp.(-1im * dt * k * transpose(e))

        transpose!(fᵀ, f)
        fft!(fᵀ, 1)
        fᵀ .= fᵀ .* ek
        ifft!(fᵀ, 1)
        transpose!(f, fᵀ)

    else

        L = mesh1.stop - mesh1.start
        m = mesh1.length ÷ 2
        k = 2π/L * [0:1:m-1;-m:1:-1]
        v = mesh2.points
        ev = exp.(-1im * dt * k * transpose(v))

        fft!(f,1)
        f   .= f .* ev
        k[1] = 1.0
        e   .= -1im * mesh2.step * vec(sum(f,dims=2)) ./ k
        e[1] = 0im
        ifft!(f,1)
        ifft!(e)

     end

end


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

@inline initMethod(settings::Settings = defaultSettings) =
    initMethod(Val(settings.initMethod), settings)

"""
    initMethod(::Val{:uniform_min_max})

Returns an initializer function.

The values are uniformly distributed over ``[\\min(c), \\max(c)]`` of the data
column.
"""
function initMethod(::Val{:uniform_min_max}, settings::Settings = defaultSettings)
    col -> rand(Uniform(minimum(col), maximum(col)), prod(settings.latticeSize))
end

"""
    initMethod(::Val{:uniform_mean_std})

Returns an initializer function.

The values are uniformly distributed between ``[μ-2σ, μ+2σ]`` of the data
column.
"""
function initMethod(::Val{:uniform_mean_std}, settings::Settings = defaultSettings)
    m = mean(col)
    s = std(col)
    col -> rand(Uniform(m - 2s, m + 2s), prod(settings.latticeSize))
end

"""
    initMethod(::Val{:normal_mean_std})

Returns an initializer function.

The values are Normal distributed with the same mean and std as the data column.
"""
function initMethod(::Val{:normal_mean_std}, settings::Settings = defaultSettings)
    col -> rand(Normal(mean(col), std(col)), prod(settings.latticeSize))
end

"""
    initMethod(::Val{:zeros})

Returns an initializer function.

The values are set to 0.
"""
function initMethod(::Val{:zeros}, settings::Settings = defaultSettings)
    col -> zeros(prod(settings.latticeSize))
end


text:
# Autogenerated wrapper script for CBFlib_small_jll for aarch64-apple-darwin
export libcbf

JLLWrappers.@generate_wrapper_header("CBFlib_small")
JLLWrappers.@declare_library_product(libcbf, "@rpath/libcbf_small.3.0.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcbf,
        "lib/libcbf_small.3.0.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libspatialite_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libspatialite")
JLLWrappers.@generate_main_file("libspatialite", UUID("a4c9d1fc-8f32-52c2-b7e7-251980156b63"))
end  # module libspatialite_jll


text:
function check(uf :: UnionFinder, groups :: Vector{Vector{Int}})
    # Check node count.

    if sum(length, groups) != length(uf)
        println("Found $(length(uf)), but expected $(sum(length, groups)).")
        return false
    end

    # Check that all nodes which should be in groups together are, in fact,
    # in groups together.

    for group in groups
        if !all(id -> find!(uf, id) == find!(uf, group[1]), group)
            ids = map(id -> find!(uf, id), group)
            println("Group $group has ids $ids.")
            return false
        elseif !all(id -> length(group) == size!(uf, id), group)
            sizes = map(id -> size!(uf, id), group)
            println("Group $group has sizes $sizes.")
            return false
        end
    end

    # Check uniqueness of groups.

    for i in 1:length(groups)
        for j in 1:length(groups)
            if i == j
                continue
            elseif find!(uf, groups[i][1]) == find!(uf, groups[j][1])
                println("Group $(groups[i]) and $(groups[j]) have same id.")
                return false
            end
        end
    end

    return true
end

function check(cf :: CompressedFinder, group_set :: Vector{Vector{Int}})
    # Check node count and group count.

    if sum(length, group_set) != length(cf)
        println("Found $(length(cf)), but expected $(sum(length, group_set)).")
        return false
    elseif length(group_set) != groups(cf)
        println("Expected $(length(groups)) groups, but got $(groups(cf)).")
        return false
    end

    # Check that all nodes which should be in groups together are, in fact,
    # in groups together.

    for group in group_set
        if !all(id -> find(cf, id) == find(cf, group[1]), group)
            ids = map(id -> find(cf, id), group)
            println("Group $group has ids $ids.")
            return false
        elseif find(cf, group[1]) <= 0 || find(cf, group[1]) > groups(cf)
            println("Group $group has invlaid id, $(find(cf, group[1]))")
            return false
        end
    end

    # Check uniqueness of groups.

    for i in 1:length(group_set)
        for j in 1:length(group_set)
            if i == j
                continue
            elseif find(cf, group_set[i][1]) == find(cf, group_set[j][1])
                println("Group $(group_set[i]) and $(group_set[j]) " +
                        "have same id.")
                return false
            end
        end
    end

    return true
end

function sc_println(test_num :: Integer, uf :: UnionFinder)
    println("Test $test_num:")
    print("parents: [")
    for (i, id) in enumerate(uf.parents)
        print("$i: $id")
        if i != length(uf.parents)
            print(", ")
        end
    end
    println("]")

    print("sizes: [")
    for (i, size) in enumerate(uf.sizes)
        print("$i: $size")
        if i != length(uf.sizes)
            print(", ")
        end
    end
    println("]")

    return true
end

function sc_println(test_num :: Integer, cf :: CompressedFinder)
    println("Test $test_num:")
    print("ids: [")
    for (i, id) in enumerate(cf.ids)
        print("$i: $id")
        if i != length(cf.ids)
            print(", ")
        end
    end
    println("]")
    return true
end

tests = [# singleton
         (Int[], Int[], [[1]]),
         # self loop
         ([1], [1], [[1]]),
         # unconnected pair
         (Int[], Int[], [[1],[2]]),
         # connected pair
         ([1], [2], [[1,2]]),
         # floater
         ([1], [2], [[1,2],[3]]),
         # "V"
         ([1,2], [2,3], [[1,2,3]]),
         # triangle
         ([1,2,3], [3,1,2], [[1,2,3]]),
         # "Z": order 1
         ([1,2,3], [2,3,4], [[1,2,3,4]]),
         # "Z": order 2
         ([1,4,2], [2,3,3], [[1,2,3,4]]),
         # diamond
         ([1,2,4,3,1], [2,4,3,4,3], [[1,2,3,4]]),
         # K-5
         ([1,1,1,1], [2,3,4,5], [[1,2,3,4,5]]),
         ([2,2,2,2], [1,3,4,5], [[1,2,3,4,5]]),
         ([3,3,3,3], [1,2,4,5], [[1,2,3,4,5]]),
         ([4,4,4,4], [1,2,3,5], [[1,2,3,4,5]]),
         ([5,5,5,5], [1,2,3,4], [[1,2,3,4,5]]),
         # jewish star
         ([1,3,5,6,4,2], [3,5,1,2,6,4], [[1,3,5],[2,4,6]]),
         # barbell
         ([1,2,3,4,5,6,6], [2,3,1,5,6,4,3], [[1,2,3,4,5,6]])
         ]

function table_main()
    for (i, (us, vs, groups)) in enumerate(tests)
        uf = UnionFinder(sum(length, groups))
        union!(uf, us, vs)
        cf = CompressedFinder(uf)

        @test check(uf, groups) || !sc_println(i, uf)
        @test check(cf, groups) || !sc_println(i, cf)

        reset!(uf)
        for (u, v) in zip(us, vs)
            union!(uf, u, v)
        end

        @test check(uf, groups) || !sc_println(i, uf)

        reset!(uf)
        union!(uf, zip(us, vs))

        @test check(uf, groups) || !sc_println(i, uf)
    end
end

table_main()


text:
function getImportantWordsInTopic( topic, dictionary, numToGet)
#gets the words in the topic with the heighest weights, and outputs the
#actual words in them and their associated weights
    topWords = Array{String}(numToGet, 1);
    weights = zeros(numToGet, 1);
#    [topic inds] = sort(topic, 'descend');

    inds = reverse(sortperm(topic));
    topic = topic[inds];

    for i = 1:numToGet
        if topic[i] == 0
            topWords = topWords[1:i-1];
            weights = weights[1:i-1];
            break;
        end
        topWords[i] = dictionary[inds[i]];
        weights[i] = topic[i];
    end

topWords, weights
end

function getImportantWordsInAllTopics(sample, dictionary, numToGet)
    numTopics = size(sample.phi,2);
    topWords = Array{String}(numToGet, numTopics);
    topicWordProbs = zeros(numToGet, numTopics);
    for i = 1:numTopics
        topWords[:,i], topicWordProbs[:,i] = getImportantWordsInTopic( sample.phi[:,i], dictionary, numToGet);
    end
    topWords, topicWordProbs
end


text:
##
using MRIgeneralizedBloch
using DifferentialEquations
using StaticArrays
using Test

## choose random parameters
# this test should only be passed for small α and long TRF
α = rand() * π/2
TRF = 500e-6 + 1000e-6 * rand()
ω1 = α / TRF
B1 = 0.7 + 0.6 * rand()
ω0 = 0 # Graham's spectral model is only implemented for on-resonance pulses
m0s = 0.2 * rand()
m0f = 1 - m0s
R1f = 0.5 + 0.2 * rand()
R1s = 2 + rand()
R2f = 1 / (40e-3 + 100e-3 * rand())
T2s = 5e-6 + 10e-6 * rand()
R2s = 1 / T2s
Rx = 100 * rand()

ms = m0s * rand()
mf = (1 - m0s) * rand()
ϑ = rand() * π / 2
φ = rand() * 2π

## Solve Sled's solution
g = interpolate_greens_function(greens_superlorentzian, 0, TRF/T2s)
p = (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g)
alg = Tsit5()
m0 = [mf * sin(ϑ) * cos(φ), mf * sin(ϑ) * sin(φ), mf * cos(ϑ), ms, 1]

sol = solve(
    ODEProblem(MRIgeneralizedBloch.apply_hamiltonian_sled!, m0, (0.0, TRF), p),
    alg,
)

u_Sled = sol[end]

## Solve generalized Bloch-McConnell with super-Lorentzian lineshape
mfun(p, t; idxs = nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5)

p = (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g)
alg = MethodOfSteps(DP8())
sol = solve(
    DDEProblem(MRIgeneralizedBloch.apply_hamiltonian_gbloch!, m0, mfun, (0.0, TRF), p),
    alg,
)

u_gBloch = sol[end]

## Test!
max_error = 1e-2
@test u_gBloch[1] ≈ u_Sled[1] atol = max_error
@test u_gBloch[2] ≈ u_Sled[2] atol = max_error
@test u_gBloch[3] ≈ u_Sled[3] atol = max_error
@test u_gBloch[4] ≈ u_Sled[4] atol = max_error
@test u_gBloch[5] ≈ u_Sled[5] atol = max_error


text:
## Least squares solutions
a = [ones(20) 1:20 1:20]
b = reshape(eye(8, 5), 20, 2)
for elty in (Float32, Float64, Complex64, Complex128)
    a = convert(Matrix{elty}, a)
    b = convert(Matrix{elty}, b)

    # Vector rhs
    x = a[:,1:2]\b[:,1]
    @test_approx_eq ((a[:,1:2]*x-b[:,1])'*(a[:,1:2]*x-b[:,1]))[1] convert(elty, 2.546616541353384)

    # Matrix rhs
    x = a[:,1:2]\b
    @test_approx_eq det((a[:,1:2]*x-b)'*(a[:,1:2]*x-b)) convert(elty, 4.437969924812031)

    # Rank deficient
    x = a\b
    @test_approx_eq det((a*x-b)'*(a*x-b)) convert(elty, 4.437969924812031)

    # Underdetermined minimum norm
    x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1])
    @test_approx_eq x convert(Vector{elty}, [1, 0.5, -0.5])

    # symmetric, positive definite
    @test_approx_eq inv(convert(Matrix{elty}, [6. 2; 2 1])) convert(Matrix{elty}, [0.5 -1; -1 3])

    # symmetric, indefinite
    @test_approx_eq inv(convert(Matrix{elty}, [1. 2; 2 1])) convert(Matrix{elty}, [-1. 2; 2 -1]/3)
end

## Test Julia fallbacks to BLAS routines

# matrices with zero dimensions
@test ones(0,5)*ones(5,3) == zeros(0,3)
@test ones(3,5)*ones(5,0) == zeros(3,0)
@test ones(3,0)*ones(0,4) == zeros(3,4)
@test ones(0,5)*ones(5,0) == zeros(0,0)
@test ones(0,0)*ones(0,4) == zeros(0,4)
@test ones(3,0)*ones(0,0) == zeros(3,0)
@test ones(0,0)*ones(0,0) == zeros(0,0)

# 2x2
A = [1 2; 3 4]
B = [5 6; 7 8]
@test A*B == [19 22; 43 50]
@test At_mul_B(A, B) == [26 30; 38 44]
@test A_mul_Bt(A, B) == [17 23; 39 53]
@test At_mul_Bt(A, B) == [23 31; 34 46]
Ai = A+(0.5*im).*B
Bi = B+(2.5*im).*A[[2,1],[2,1]]
@test Ai*Bi == [-21+53.5im -4.25+51.5im; -12+95.5im 13.75+85.5im]
@test Ac_mul_B(Ai, Bi) == [68.5-12im 57.5-28im; 88-3im 76.5-25im]
@test A_mul_Bc(Ai, Bi) == [64.5+5.5im 43+31.5im; 104-18.5im 80.5+31.5im]
@test Ac_mul_Bc(Ai, Bi) == [-28.25-66im 9.75-58im; -26-89im 21-73im]

# 3x3
A = [1 2 3; 4 5 6; 7 8 9].-5
B = [1 0 5; 6 -10 3; 2 -4 -1]
@test A*B == [-26 38 -27; 1 -4 -6; 28 -46 15]
@test Ac_mul_B(A, B) == [-6 2 -25; 3 -12 -18; 12 -26 -11]
@test A_mul_Bc(A, B) == [-14 0 6; 4 -3 -3; 22 -6 -12]
@test Ac_mul_Bc(A, B) == [6 -8 -6; 12 -9 -9; 18 -10 -12]
Ai = A+(0.5*im).*B
Bi = B+(2.5*im).*A[[2,1,3],[2,3,1]]
@test Ai*Bi == [-44.75+13im 11.75-25im -38.25+30im; -47.75-16.5im -51.5+51.5im -56+6im; 16.75-4.5im -53.5+52im -15.5im]
@test Ac_mul_B(Ai, Bi) == [-21+2im -1.75+49im -51.25+19.5im; 25.5+56.5im -7-35.5im 22+35.5im; -3+12im -32.25+43im -34.75-2.5im]
@test A_mul_Bc(Ai, Bi) == [-20.25+15.5im -28.75-54.5im 22.25+68.5im; -12.25+13im -15.5+75im -23+27im; 18.25+im 1.5+94.5im -27-54.5im]
@test Ac_mul_Bc(Ai, Bi) == [1+2im 20.75+9im -44.75+42im; 19.5+17.5im -54-36.5im 51-14.5im; 13+7.5im 11.25+31.5im -43.25-14.5im]

# Generic integer matrix multiplication
A = [1 2 3; 4 5 6] .- 3
B = [2 -2; 3 -5; -4 7]
@test A*B == [-7 9; -4 9]
@test At_mul_Bt(A, B) == [-6 -11 15; -6 -13 18; -6 -15 21]
A = ones(Int, 2, 100)
B = ones(Int, 100, 3)
@test A*B == [100 100 100; 100 100 100]
A = rand(1:20, 5, 5) .- 10
B = rand(1:20, 5, 5) .- 10
@test At_mul_B(A, B) == A'*B
@test A_mul_Bt(A, B) == A*B'

# Preallocated
C = Array(Int, size(A, 1), size(B, 2))
@test A_mul_B!(C, A, B) == A*B
@test At_mul_B!(C, A, B) == A'*B
@test A_mul_Bt!(C, A, B) == A*B'
@test At_mul_Bt!(C, A, B) == A'*B'

# matrix algebra with subarrays of floats (stride != 1)
A = reshape(float64(1:20),5,4)
Aref = A[1:2:end,1:2:end]
Asub = sub(A, 1:2:5, 1:2:4)
b = [1.2,-2.5]
@test (Aref*b) == (Asub*b)
@test At_mul_B(Asub, Asub) == At_mul_B(Aref, Aref)
@test A_mul_Bt(Asub, Asub) == A_mul_Bt(Aref, Aref)
Ai = A .+ im
Aref = Ai[1:2:end,1:2:end]
Asub = sub(Ai, 1:2:5, 1:2:4)
@test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref)
@test A_mul_Bc(Asub, Asub) == A_mul_Bc(Aref, Aref)

# syrk & herk
A = reshape(1:1503, 501, 3).-750.0
res = float64([135228751 9979252 -115270247; 9979252 10481254 10983256; -115270247 10983256 137236759])
@test At_mul_B(A, A) == res
@test A_mul_Bt(A',A') == res
cutoff = 501
A = reshape(1:6*cutoff,2*cutoff,3).-(6*cutoff)/2
Asub = sub(A, 1:2:2*cutoff, 1:3)
Aref = A[1:2:2*cutoff, 1:3]
@test At_mul_B(Asub, Asub) == At_mul_B(Aref, Aref)
Ai = A .- im
Asub = sub(Ai, 1:2:2*cutoff, 1:3)
Aref = Ai[1:2:2*cutoff, 1:3]
@test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref)

# Matrix exponential
for elty in (Float32, Float64, Complex64, Complex128)
        A1  = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4])
        eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181  127.781085523182;
        183.765138646367 183.765138646366  163.679601723179;
        71.797032399996  91.8825693231832 111.968106246371]')
        @test_approx_eq expm(A1) eA1

        A2  = convert(Matrix{elty},
            [29.87942128909879    0.7815750847907159 -2.289519314033932;
            0.7815750847907159 25.72656945571064    8.680737820540137;
            -2.289519314033932   8.680737820540137  34.39400925519054])
        eA2 = convert(Matrix{elty},
            [  5496313853692458.0 -18231880972009236.0 -30475770808580460.0;
             -18231880972009252.0  60605228702221920.0 101291842930249760.0;
             -30475770808580480.0 101291842930249728.0 169294411240851968.0])
        @test_approx_eq expm(A2) eA2

        A3  = convert(Matrix{elty}, [-131 19 18;-390 56 54;-387 57 52])
        eA3 = convert(Matrix{elty}, [-1.50964415879218 -5.6325707998812  -4.934938326092;
        0.367879439109187 1.47151775849686  1.10363831732856;
        0.135335281175235 0.406005843524598 0.541341126763207]')
        @test_approx_eq expm(A3) eA3

        # issue 5116
        A4  = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
        eA4 = [-0.999786072879326  -0.065407069689389   0.0   0.0
                0.006540706968939  -0.999786072879326   0.0   0.0
                0.0                 0.0                 1.0   0.0
                0.013081413937878  -3.999572145758650   0.0   1.0]
        @test_approx_eq expm(A4) eA4

        # issue 5116
        A5  = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.]
        eA5 = [ 1.0+0.0im   0.0+0.0im                 0.0+0.0im                0.0+0.0im
                0.0+0.0im   1.543080634815244+0.0im   0.0-1.175201193643801im  0.0+0.0im
                0.0+0.0im   0.0+1.175201193643801im   1.543080634815243+0.0im  0.0+0.0im
                0.0+0.0im   0.0+0.0im                 0.0+0.0im                1.0+0.0im]
        @test_approx_eq expm(A5) eA5

        # Hessenberg
        @test_approx_eq hessfact(A1)[:H] convert(Matrix{elty},
                        [4.000000000000000  -1.414213562373094  -1.414213562373095
                        -1.414213562373095   4.999999999999996  -0.000000000000000
                                         0  -0.000000000000002   3.000000000000000])
end

# Hermitian matrix exponential
A1 = randn(4,4) + im*randn(4,4)
A2 = A1 + A1'
@test_approx_eq expm(A2) expm(Hermitian(A2))

# matmul for types w/o sizeof (issue #1282)
A = Array(Complex{Int},10,10)
A[:] = complex(1,1)
A2 = A^2
@test A2[1,1] == 20im

# test sparse matrix norms
Ac = sprandn(10,10,.1) + im* sprandn(10,10,.1)
Ar = sprandn(10,10,.1)
Ai = int(ceil(Ar*100))
@test_approx_eq norm(Ac,1)     norm(full(Ac),1)
@test_approx_eq norm(Ac,Inf)   norm(full(Ac),Inf)
@test_approx_eq vecnorm(Ac)    vecnorm(full(Ac))
@test_approx_eq norm(Ar,1)     norm(full(Ar),1)
@test_approx_eq norm(Ar,Inf)   norm(full(Ar),Inf)
@test_approx_eq vecnorm(Ar)    vecnorm(full(Ar))
@test_approx_eq norm(Ai,1)     norm(full(Ai),1)
@test_approx_eq norm(Ai,Inf)   norm(full(Ai),Inf)
@test_approx_eq vecnorm(Ai)    vecnorm(full(Ai))

# scale real matrix by complex type
@test_throws ErrorException scale!([1.0], 2.0im)
@test isequal(scale([1.0], 2.0im),             Complex{Float64}[2.0im])
@test isequal(scale(Float32[1.0], 2.0f0im),    Complex{Float32}[2.0im])
@test isequal(scale(Float32[1.0], 2.0im),      Complex{Float64}[2.0im])
@test isequal(scale(Float64[1.0], 2.0f0im),    Complex{Float64}[2.0im])
@test isequal(scale(Float32[1.0], big(2.0)im), Complex{BigFloat}[2.0im])
@test isequal(scale(Float64[1.0], big(2.0)im), Complex{BigFloat}[2.0im])
@test isequal(scale(BigFloat[1.0], 2.0im),     Complex{BigFloat}[2.0im])
@test isequal(scale(BigFloat[1.0], 2.0f0im),   Complex{BigFloat}[2.0im])

# issue #6450
@test dot({1.0,2.0},{3.5,4.5}) === 12.5

# Issue 7181
A = reshape([1:9], (3, 3))
@test_throws BoundsError diag(A, 3)
@test [7] == diag(A, 2)
@test [4, 8] == diag(A, 1)
@test [1,5,9] == diag(A, 0)
@test [2,6] == diag(A, -1)
@test [3] == diag(A, -2)
@test_throws BoundsError diag(A, -3)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# Compute smFISH histogram - classic telegraph model for arbtirary numbers of gene states
function telegraph(n::Int,rin::Vector,total::Int,tol::Float64,nhist::Int,nalleles::Int,count=false)

	# n-state telegraph model for nalleles alleles
	# Compute steady state mRNA histogram
	# n+1 total gene states
	# gene states are labeled from 1 to n+1
	# There are 2n gene reactions, 1 mRNA creation reaction, and 1 mRNA decay reaction
	# rin = reaction rates
	# forward gene reactions are labeled by odd numbers <2n, backward reactions are even numbers <= 2n
	# reactions 2n+1 mRNA creation reaction, reaction 2n+2 is mRNA decay reaction

	# Gillespie direct method with speed up from Gibson and Bruck

	#Initialize m histogram
	mhist = zeros(nhist+1)

	# active start site is n+1
	mu = n
	#print(n," ",mu," ",zeta,'\n')

	# Initialize rates
	gammaf = zeros(n+1)
	gammab = zeros(n+1)

	# assign gene rates
	for i = 1:n
		gammaf[i] = rin[2*(i-1)+1]
		gammab[i+1] = rin[2*i]
	end

	# mRNA creation rate
	nu = rin[2*n + 1]

	# assign decay rate
	# decay reaction index
	reactiondecay = 2*n + 2
	delta = rin[reactiondecay]

	# allele coupling strength
	if length(rin) > reactiondecay
		rcorr = rin[reactiondecay+1]
	else
		rcorr = 0.
	end

	#initialize mRNA number to 0
	m = 0

	# Initialize gene state to most inactive state
	state = ones(Int,nalleles)

	# initialize propensities
	af = Array{Float64,1}(undef,nalleles)  # forward
	ab = Array{Float64,1}(undef,nalleles)  # backward
	afree = nu*ones(nalleles).*(state .== 1)    # mRNA creation
	adecay = 0.       # mRNA decay

	for i = 1:nalleles
		af[i] = gammaf[state[i]]
		ab[i] = gammab[state[i]]
	end
	astate = af + ab

	steps = 0

	# Begin simulation
	t=0.
	tsample = 20/minimum(rin)
	sample = 0
	err = 10*tol
	mhist0 = ones(nhist+1)

	while err > tol && steps < total

		steps += 1

		# sum of all propensities
		asum = sum(afree) + sum(astate) + adecay

		if asum > 0.
			# choose which state to update and time of update
			r = asum*rand()
			tau = -log(rand())/asum

			# update sample time
			t += tau
			#print(tau,' ',asum,' ',state,'\n')

			# update m histogram
			if m + 1 <= nhist
				mhist[m+1] += tau
			else
				mhist[nhist+1] += tau
			end

			# Compute change in histogram from last sample

			if t > tsample
				err = norm(mhist/sum(mhist)-mhist0/sum(mhist0),Inf)
				mhist0 = copy(mhist)
				t = 0.
			end

			if r > sum(afree + astate)
				# mRNA decay reaction
				m -= 1
				adecay = m*delta
			else
				# gene state reaction
				# find gene that updated
				agenecs = afree[1] + astate[1]
				agene0 = 0.
				l = 1
				while r > agenecs
					agene0 = agenecs
					agenecs += afree[1+l] + astate[l+1]
					l += 1
				end
				if l > 1
					r -= agene0
				end

				if r > afree[l]
					# state update
					if r > ab[l] + afree[l]
						# update forward reaction
						state[l] += 1
					else
						# update backward reaction
						state[l] -= 1
					end

					# update propensities
					af[l] = gammaf[state[l]]
					ab[l] = gammab[state[l]]
					astate[l] = af[l] + ab[l]

					# gene state is active if > mu
					if state[l] > mu
						afree[l] = nu
					else
						afree[l] = 0.
					end

					# Activate coupling if any alleles are in state mu + 1
					if rcorr > 0.
						ba = false
						for k = 1:nalleles
							ba |= state[k] > mu
						end
						if ba
							for k = 1:nalleles
								if state[k] == mu
									af[k] = gammaf[state[k]] + rcorr
									astate[k] = af[k] + ab[k]
								end
							end
						else
							for k = 1:nalleles
								if state[k] == mu
									af[k] = gammaf[state[k]]
									astate[k] = af[k] + ab[k]
								end
							end
						end
					end
				else
					# increase free mRNA
					m += 1
					adecay = m*delta
				end # if r > free[l]
			end #if r > sum(afree + astate)

		end  #if asum > 0

	end  # while

	counts = max(sum(mhist),1)
	mhist /= counts

	if count
		return mhist[1:nhist],counts,steps,err
	else
		return mhist[1:nhist]
	end

end

function prep_params(n,r)
	gammap = zeros(n+2)
	gamman = zeros(n+2)
	for i = 1:n
		gammap[i+1] = r[2*(i-1)+1]
		gamman[i+1] = r[2*i]
	end
	return gammap, gamman
end

function transitiontelegraph(n,rin,nhist)
	nT = n+1
	total = nhist + 2
	r = rin/rin[2*n + 2]
	gammap,gamman = prep_params(n,r)  # State Transition rates
	nu = r[2*n+1] # mRNA creation rate
	T = zeros(nT,nT)
	B = zeros(nT,nT)
	S = zeros(total,total)
	Sminus = similar(S)
	Splus = similar(S)
	# Generate T = A + B matrix and B matrix
	for ip=1:nT, i=1:nT
		T[i,ip] =  gammap[i]*float(i-1==ip) - (gamman[i]+gammap[i+1])*float(i==ip) + gamman[i+1]*float(i+1==ip)
		B[i,ip] =  nu*float(i==ip)*float(i==n+1)
	end
	# Generate matrices for m transitions
	for m=1:total,mp=1:total
		S[m,mp] = -(mp)*float(m==mp)
		Sminus[m,mp] = float(m==mp+1)
		Splus[m,mp] = (mp-1)*float(m==mp-1)
	end
	# Construct m truncated Master equation rate matrix
	M = kron(S,Matrix{Float64}(I,nT,nT)) + kron(Matrix{Float64}(I,total,total),T-B) + kron(Sminus,B) + kron(Splus,Matrix{Float64}(I,nT,nT))
	M[end,end]+=nu  # boundary condition to ensure probability is conserved
	return M
end

function transitiontelegraph(r,nhist)
	nu = r[1]/r[2]
	total = nhist + 2
	S = zeros(total,total)
	for m=1:total,mp=1:total
		S[m,mp] = nu *(m==mp+1) - nu*(m==mp) + (mp-1)*(m==mp-1) - (mp-1)*(m==mp)
	end
	S[end,end] += nu
	return S
end

# function marginalize(p,n,nhist)
# 	mhist = zeros(nhist)
# 	nT = n+1
# 	# Marginalize over G states
# 	for m in 1:nhist
# 		i = (m-1)*nT
# 		mhist[m] = sum(p[i+1:i+nT])
# 	end
# 	return mhist
# end

function transientMaster(t::Float64,n::Int,P0,Mvals::Vector,Mvects::Matrix)
	Lambda = Diagonal(exp.(Mvals*t))
	Q = V/P0
	P = V*Lambda*Q
end

# function sstel(n,r,nhist,nalleles)
# 	M = transitiontelegraph(n,r,nhist)
# 	p = normalizednullspace(M)
# 	mhist = marginalize(p,n,nhist)
# 	alleleConvolve(mhist,nalleles)
# end
#
# function sstel(r,nhist,nalleles)
# 	S = transitiontelegraph(r,nhist)
# 	mhist = normalizednullspace(S)
# 	alleleConvolve(mhist,nalleles)
# end

function sstel2(n::Int,rin::Vector{Float64},nhist::Int,nalleles::Int)
	# Closed form steady state mRNA distribution of generalized n-state classic telegraph model
	nT = n+1
	r = rin/rin[2*n + 2]
	total = nhist + 2
	# Transition rates
	gammap,gamman = prep_params(n,r)
	# mRNA creation rate
	nu = r[2*n+1]
	# Declare Transition matrices
	T = zeros(nT,nT)
	B = zeros(nT,nT)
	S = zeros(total,total)
	Sminus = similar(S)
	Splus = similar(S)
	# Generate T = A + B matrix and B matrix
	for i=1:n+1,  ip=1:n+1
		T[i,ip] =  gammap[i]*float(i-1==ip) - (gamman[i]+gammap[i+1])*float(i==ip) + gamman[i+1]*float(i+1==ip)
		B[i,ip] =  nu*float(i==ip)*float(i==n+1)
	end
	# Generate matrices for m transitions
	for m=1:total,mp=1:total
		S[m,mp] = -(mp-1)*float(m==mp)
		Sminus[m,mp] = float(m==mp+1)
		Splus[m,mp] = (mp-1)*float(m==mp-1)
	end
	# Construct m truncated Master equation rate matrix
	M = kron(S,Matrix{Float64}(I,nT,nT)) + kron(Matrix{Float64}(I,total,total),T-B) + kron(Sminus,B) + kron(Splus,Matrix{Float64}(I,nT,nT))
	# boundary condition to ensure probability is conserved
	M[end,end]+=nu

	# Compute steady state solution of Master equation
	# MS = try svd(M)
	# p = nullspace(M)[:,1]
	# p /= sum(p)
	p = normalizednullspace(M)
	# Marginalize over G states
	mhist = marginalize(p,n,nhist)
	alleleConvolve(mhist,nalleles)
end

function sstelFISH(n::Int,rin::Vector{Float64},nhist::Int,nalleles::Int)
	# Closed form steady state mRNA distribution of generalized n-state classic telegraph model
	r = copy(rin)
	nT = n+1
	reactiondecay = 2*n + 2
	r /= rin[reactiondecay]
	total = nhist + 10
	mhist = zeros(nhist)
	# Transition rates
	gammap = zeros(n+2)
	gamman = zeros(n+2)
	for i = 1:n
		gammap[i+1] = r[2*(i-1)+1]
		gamman[i+1] = r[2*i]
	end
	# mRNA creation rate
	nu = r[2*n+1]
	if n > 0
		# Declare Transition matrices
		T = zeros(nT,nT)
		B = zeros(nT,nT)
		S = zeros(total,total)
		Sminus = similar(S)
		Splus = similar(S)
		# Generate T = A + B matrix and B matrix
		for i=1:n+1,  ip=1:n+1
			T[i,ip] =  gammap[i]*float(i-1==ip) - (gamman[i]+gammap[i+1])*float(i==ip) + gamman[i+1]*float(i+1==ip)
			B[i,ip] =  nu*float(i==ip)*float(i==n+1)
		end
		# Generate matrices for m transitions
		for m=1:total,mp=1:total
			S[m,mp] = -(mp-1)*float(m==mp)
			Sminus[m,mp] = float(m==mp+1)
			Splus[m,mp] = (mp-1)*float(m==mp-1)
		end
		# Construct m truncated Master equation rate matrix
		M = kron(S,Matrix{Float64}(I,nT,nT)) + kron(Matrix{Float64}(I,total,total),T-B) + kron(Sminus,B) + kron(Splus,Matrix{Float64}(I,nT,nT))
		# boundary condition to ensure probability is conserved
		M[end,end]+=nu

		# Compute steady state solution of Master equation
		# MS = try svd(M)
		p = nullspace(M)[:,1]
		# catch
		# 	tol = 1e-8
		# 	return telegraph(n,r,100000000,tol,nhist,nalleles)
		# end
		# p = MS.Vt[end,:]
		p /= sum(p)
		# Marginalize over G states
		for n in 1:nhist
			i = (n-1)*nT
			mhist[n] = sum(p[i+1:i+nT])
		end
	else
		# 0-state telegraph model
		# total = nhist + 10
		S = zeros(total,total)
		for m=1:total,mp=1:total
			S[m,mp] = nu *(m==mp+1) - nu*(m==mp) + (mp-1)*(m==mp-1) - (mp-1)*(m==mp)
		end
		S[end,end] += nu
		mhist = nullspace(S)[:,1]
		# catch
		# 	tol = 1e-5
		# 	mhist = telegraph(n,r,10000000,tol,nhist,2)
		# end
		mhist /= sum(mhist)
	end
	# Convolve to compute distribution for contributions from multiple alleles
	mhists = Array{Array{Float64,1}}(undef,nalleles)
	mhists[1] = convert.(Float64,mhist)
	for i = 2:nalleles
		mhists[i] = zeros(nhist)
		for m = 0:nhist-1
			for m2 = 0:min(nhist-1,m)
				mhists[i][m+1] += mhists[i-1][m-m2+1]*mhist[m2+1]
			end
		end
	end
	return mhists[nalleles]
end

function transientFISH(t::Float64,n::Int,P0,Mvals,Mvects,rin::Vector{Float64},nhist::Int,nalleles::Int)
	# Closed form steady state mRNA distribution of generalized n-state classic telegraph model
	P = transientMaster(t,Mvals,Mvects)
	mhist= marginalize(P,nhist,n)
	alleleConvolve(mhist,nalleles)
end

max_stars_repo_name: nih-niddk-mbs/StochasticGene

max_issues_repo_name: nih-niddk-mbs/StochasticGene

max_forks_repo_name: nih-niddk-mbs/StochasticGene


text:
__precompile__(true)
module Baz
import Foo, Qux
this = "Baz"
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
function selector!(z,x)
  @inbounds for i ∈ eachindex(x)
    z[i] = x[i] > 0.5 ? i : x[i]
  end
  z
end
function selector_turbo!(z,x)
  @turbo for i ∈ eachindex(x)
    z[i] = x[i] > 0.5 ? i : x[i]
  end
  z
end

# testset for using them in loops
@testset "Loop Induction Variables" begin
  @show @__LINE__
  f(x) = cos(x) * log(x)
  function avxmax(v)
    max_x = -Inf
    @turbo for i ∈ eachindex(v)
      x = f(i)
      max_x = max(max_x, x)
    end
    max_x
  end
  function avxextrema(v)
    max_x = -Inf
    min_x = Inf
    @turbo for i ∈ eachindex(v)
      x = f(i)
      max_x = max(max_x, x)
      min_x = min(min_x, x)
    end
    min_x, max_x
  end

  v = 1:19
  minref, maxref = extrema(f, v)
  @test maxref ≈ avxmax(v)
  minavx, maxavx = avxextrema(v);
  @test minref ≈ minavx
  @test maxref ≈ maxavx

  x = rand(119); y0 = similar(x); y1 = similar(y0);
  @test selector!(y0, x) == selector_turbo!(y1, x)
end


text:
using Revise

includet("src/Spring.jl")
IArray = InternalArray

# 3. Validating the solution --------------------------------------------------------------

const step = 1e-3
const nit = Int(40/step)  # round to avoid small rounding errors

# We define a problem with no friction
prob = SpringProblem(p0 = 0.8, v0 = 0, m = 1, k = 1.5)

x = eulermethod(prob, step, nit)

timeline = [prob.t0 + i*step for i in 0:length(x)-1]
plot(timeline, IArray(x, 1), label = "euler method")

# And also define the solution we would have obtained by analytic methods
# Reference: https://en.wikipedia.org/wiki/Simple_harmonic_motion
ω0 = sqrt(prob.k/prob.m)
g(t) = prob.p0*cos(ω0*t) + prob.v0/ω0*sin(ω0*t)

y = g.(timeline)
plot!(timeline, y, label = "analytic solution")

# (The comparison in my computer is very close)
# (Which means that for the second step we could have taken a smaller L very likely)
# Let's check the norm of the difference:

err = maximum(norm.(y-IArray(x, 1)))

# 4. Choosing the step to obtain a final error in the order of 1e-2  ----------------------

# Reference: https://en.wikipedia.org/wiki/Euler_method#Local_truncation_error
# The global truncation error is bounded by $\frac{hM}{2L}(ℯ^{L(t-t_0)}-1)$ where
# M is an upper bound for the second derivative of the solution (in the interval) and
# L is the lipschitz constant of f.

# This will be explained later, but L < k+b and M can be taken to be maximum of the
# second derivative for the approximated solution.

# <-- Copied from the definitions
n = 1
function f(t, x)
    p = x[1:n]
    v = x[n+1:2*n]
    return [v; (-prob.k*p -prob.b*v + [prob.A*sin(prob.ω*t)])/prob.m]
end

M = maximum((norm ∘ f).(timeline, x))
L = prob.k+prob.b
err_ct = M/L*(ℯ^(L*40)-1) # error constant
target_step = 10e-2/err_ct

# After checking that with 1e-3 the error is approximately $0.02$
# for the solution without friction.
# It seems natural to admit that an appropiate step to obtain
# an error less than 1e-2 is approximately 5e-4.
# This can be verified using this same script but
# with a smaller error for the method without friction.

# After doing it, I saw that it was not enough and
# the error was still slightly larger than 1e-2.
# With the step equal to 1e-4, the error for the problem without friction is $0.0023$.
# Hence, we can assume that it won't be much larger for the problem with friction.

savefig("media/spring_validation.png")

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
struct Policy
    letter::Char
    min_count::Int
    max_count::Int
end

satisfies_policy(pw::AbstractString, p::Policy) = p.min_count <= count(l -> l == p.letter, pw) <= p.max_count

passwords = split.(readlines("day2/input.txt"))
out = 0
for (range, char_text, password) in passwords
    min_count, max_count = parse.(Int, split(range, "-"))
    pw_policy = Policy(char_text[1], min_count, max_count)
    global out += satisfies_policy(password, pw_policy)
end

print(out)


text:
"wavenumber from Challis, R. E., et al. Ultrasound techniques for characterizing colloidal dispersions. Reports on progress in physics 68.7 (2005): 1541."
wavenumber_challis(ωs::AbstractArray,medium::Medium, species::Array{Specie{T}}; kws...) where T<:AbstractFloat = [wavenumber_challis(ω, medium, species; kws...) for ω in ωs]

# will sum Hankel orders only up to n =2 as in the paper "Ultrasound techniques for characterizing colloidal dispersions"
function wavenumber_challis(ω::Number, medium::Medium{T}, species::Array{Specie{T}};
    radius_multiplier = 1.005, verbose = false, hankel_order=2) where T<:AbstractFloat
  # tol=0.004; radius_multiplier = 1.005
  # a12 = radius_multiplier(a1+a2)

  volume_fraction = sum(pi*sp.r^2.0*sp.num_density for sp in species)
  if volume_fraction >= 0.4
    @warn("the volume fraction $(round(100*volume_fraction))% is a bit too high, expect a relative error of approximately $(round(100*volume_fraction^3.0))%")
  end

  kT2 = (ω/medium.c)^2.0
  kT2 += 4.0im*sum(sp.num_density*Zn(ω,sp,medium,m) for sp in species
  , m in -hankel_order:hankel_order)

  # second order number fraction, sum up too same hankel order
  kT2 += -(8.0/(ω/medium.c)^2)*sum(
    sp1.num_density^2*abs(m2-m1)*Zn(ω,sp1,medium,m1)*Zn(ω,sp1,medium,m2)
  for sp1 in species, m1 = -hankel_order:hankel_order, m2 = -hankel_order:hankel_order)

  return sqrt(kT2)
end

"low volfrac and low wavenumber. Fails badly for strong scatterers"
function one_species_low_wavenumber(ωs, medium, sp)
  φ = pi*sp.r^2*sp.num_density
  kS = 1.0/sp.c
  k = 1.0/medium.c
  P = 1 - kS^2*medium.ρ/(k^2*medium.ρ)
  Q = (medium.ρ - sp.ρ)/(medium.ρ + sp.ρ)
  (ωs./medium.c).*(1 - (φ/2)*(P+2*Q)-(φ^2/8)*(2*P^2-(P+2*Q)^2))
end

function wavenumber_far_field_low_volfrac(ω::Complex{T}, medium::Medium{T}, species::Array{Specie{T}};
        tol=0.0002, radius_multiplier = 1.005,
        verbose = false) where T<:AbstractFloat
  # tol=0.004; radius_multiplier = 1.005
  # a12 = radius_multiplier(a1+a2)

  volume_fraction = sum(pi*sp.r^2.0*sp.num_density for sp in species)
  if volume_fraction >= 0.4
    @warn("the volume fraction $(round(100*volume_fraction))% is a bit too high, expect a relative error of approximately $(round(100*volume_fraction^3.0))%")
  end
  kT2 = (ω/medium.c)^2.0
  next_order = 4.0im*sum(sp.num_density*Zn(ω,sp,medium,0) for sp in species)
  hankel_order=1

  # sum more hankel orders until the relative error < tol
  while abs(next_order/kT2) > tol
    kT2 += next_order
    hankel_order +=1
    next_order = 4.0im*sum(sp.num_density*Zn(ω,sp,medium,m) for sp in species, m in (-hankel_order,hankel_order))
  end
  kT2 += next_order
  if verbose println("max Hankel order = $hankel_order") end

  # second order number fraction, sum up too same hankel order
  kT2 += -8.0*sum(
    sp1.num_density*sp2.num_density*abs(m2-m1)/(ω/medium.c)^2*Zn(ω,sp1,medium,m1)*Zn(ω,sp2,medium,m2)
  for sp1 in species, sp2 in species, m1 = -hankel_order:hankel_order, m2 = -hankel_order:hankel_order)

  return sqrt(kT2)
end


text:
mutable struct PVSystem
    capacity::Float64
    monthly_pu_energy::Array{Float64}
    life_loss::Float64
    cost::Array{Any}
    maintenance::Float64
    time_series::Array{Float64}
end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
function make_fft_plan(x)
    xc = zeros(2 * length(x))
    p = plan_rfft(xc)

    return (xc, p)
end

function _fft_with_plan(x, r, p)
    n = length(x)
    r[1:n] = x
    result = imag(conj!(p * r))

    return result[1:(end - 1)]
end

function fft_oz(f, r, rmax, nr, p)
    for i in 1:nr
        f[i] *= (i - 1)
    end

    result = _fft_with_plan(f, r, p)

    normalization = (4.0 * rmax^3) / nr^2
    for i in 2:nr
        result[i] *= normalization / (i - 1)
    end

    return result
end

function ifft_oz(f, r, rmax, nr, p)
    for i in 1:nr
        f[i] *= (i - 1)
    end

    result = _fft_with_plan(f, r, p)

    normalization = nr * (0.5 / rmax^3)
    for i in 2:nr
        result[i] *= normalization / (i - 1)
    end

    return result
end


text:
@testset "Renewable data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type HydroDispatch, consider changing the device models"
    model = DeviceModel(PSY.HydroDispatch, PSI.HydroDispatchRunOfRiver)
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5)
    @test_logs (:warn, warn_message) construct_device!(op_model, :Hydro, model)
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys14)
    @test_logs (:warn, warn_message) construct_device!(op_model, :Hydro, model)
end

@testset "Hydro DCPLossLess FixedOutput" begin
    model = DeviceModel(PSY.HydroFix, PSI.HydroFixed)

    # Parameters Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_hy ; parameters = true)
    construct_device!(op_model, :Hydro, model)
    moi_tests(op_model, true, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)

    # No Parameters Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_hy)
    construct_device!(op_model, :Hydro, model);
    moi_tests(op_model, false, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)

    # No Forecast Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_hy ; parameters = true, forecast = false)
    construct_device!(op_model, :Hydro, model);
    moi_tests(op_model, true, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)

    # No Forecast - No Parameters Testing
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_hy ; forecast = false)
    construct_device!(op_model, :Hydro, model);
    moi_tests(op_model, false, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(op_model, GAEVF)

end

#=
@testset " Hydro Tests" begin
    PSI.activepower_variables(ps_model, generators_hg, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, DCPPowerModel, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, ACPPowerModel, 1:24)
    PSI.reactivepower_variables(ps_model, generators_hg, 1:24)
    PSI.reactivepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, ACPPowerModel, 1:24)
end

@testset " Hydro Tests" begin
    ps_model = PSI._canonical_model_init(length(sys5b.buses), nothing, PM.AbstractPowerModel, sys5b.time_periods)
    PSI.activepower_variables(ps_model, generators_hg, 1:24)
    PSI.commitment_variables(ps_model, generators_hg, 1:24);
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, DCPPowerModel, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, ACPPowerModel, 1:24)
    PSI.reactivepower_variables(ps_model, generators_hg, 1:24)
    PSI.reactivepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, ACPPowerModel, 1:24)
end
=#


text:
"""
    get_proportion_2d(for_prop_arr::AbstractArray{Float64, 2})::AbstractArray{Float64, 2}
Compute the proportional values for each element in the 2d array.
returns a AbstractArray{Float64, 2}
# Example
```
julia> get_proportion_2d([1.0 2.0 ; 4.0 5.0])
=[0.08333333333333333 0.16666666666666666; 0.3333333333333333 0.4166666666666667]
```
"""
function get_proportion_2d(for_prop_arr::AbstractArray{Float64,2})::AbstractArray{Float64,2}
    return for_prop_arr ./ sum(for_prop_arr)
end

"""
    get_proportion_round_2d(for_prop_arr::AbstractArray{Float64, 2};round_digits::Int64)::AbstractArray{Float64, 2}
Compute the proportional values for each element in the 2d array and round to a specified decimal digits.
If the sum of all the elements is not 1 , then the adjustment is done in the 1st largest element
returns a AbstractArray{Float64, 1}
# Example
```
julia> get_proportion_round_2d([1.1101 2.243 ; 4.9898 5.87],round_digits=5)
=[0.07811 0.15781; 0.35108 0.413]
```
"""
function get_proportion_round_2d(for_prop_arr::AbstractArray{Float64,2};round_digits::Int64 = 2)::AbstractArray{Float64,2}
    temparr=  round.(for_prop_arr ./ sum(for_prop_arr), digits=round_digits, base = 10)
    diff = 1 - sum(temparr) #delta naming
    if diff != 0.0
        max = argmax(temparr)
        temparr[max] =  round(temparr[max] + diff, digits=round_digits, base = 10)
    end
    return temparr
end

"""
    get_proportion_round_add1251(for_prop_arr::AbstractArray{Float64, 2};round_digits::Int64)::AbstractArray{Float64, 2}
Compute the proportional values for each element in the 2d array and round to a specified decimal digits.
If the sum of all the elements is not 1 , then add 1251 to the array and try looking for proportional correctness
after 2 itterations , the adjustment is done in the 1st largest element
returns a AbstractArray{Float64, 2}
# Example
```
julia> get_proportion_round_add1251_2d([1.1101 2.243 ; 4.9898 5.87 ; 34.333 6.78],round_digits=8)
=[0.16559399 0.16574381; 0.16610708 0.16622349; 0.16998779 0.16634384]
```
"""
function get_proportion_round_add1251_2d(for_prop_arr::AbstractArray{Float64,2};round_digits::Int64 = 2)::AbstractArray{Float64,2}
    temparr=  round.(for_prop_arr ./ sum(for_prop_arr), digits=round_digits, base = 10)
    diff = 1 - sum(temparr)
    if diff != 0.0
        check = 1
        while check < 2
            sqarr= (for_prop_arr) .+1251
            temparr=  round.(sqarr ./ sum(sqarr), digits=round_digits, base = 10)
            diff = 1 - sum(temparr)
            if diff != 0.00
                check = check + 1
                if check == 2
                    max = argmax(temparr)
                    temparr[max] = round(temparr[max]+diff,digits=round_digits, base = 10)
                end
            else
                break
            end
        end
    end
    return temparr
end

max_stars_repo_name: hillelawaskar/Proportions

max_issues_repo_name: hillelawaskar/Proportions


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function StrArrayStruct{T<:AbstractString}(strs::T...)
    strcount = length(strs)
    for s in strs
        if Base.containsnul(s)
            throw("embedded NULs are not allowed in C strings: $(repr(s))")
        end
    end
    sa_strings = convert(Ptr{Cstring}, Libc.malloc(sizeof(Cstring) * strcount))
    for i=1:strcount
        len = length(strs[i])
        in_ptr  = pointer(String(strs[i]))
        out_ptr = convert(Ptr{UInt8}, Libc.malloc(sizeof(UInt8) * (len + 1)))
        unsafe_copy!(out_ptr, in_ptr, len)
        unsafe_store!(out_ptr, zero(UInt8), len + 1) # NULL byte
        unsafe_store!(sa_strings, convert(Cstring, out_ptr), i)
    end
    return StrArrayStruct(sa_strings, strcount)
end
StrArrayStruct{T<:AbstractString}(strs::Vector{T}) = StrArrayStruct(strs...)

function Base.convert(::Type{Vector{AbstractString}}, sa::StrArrayStruct)
    arr = Array{AbstractString}(sa.count)
    for i=1:sa.count
        arr[i] = unsafe_string(unsafe_load(sa.strings, i))
    end
    return arr
end

function Base.copy(src::StrArrayStruct)
    dst_ptr = Ref(StrArrayStruct())
    @check ccall((:git_strarray_copy, :libgit2), Cint,
                  (Ptr{StrArrayStruct}, Ptr{StrArrayStruct}),
                   dst_ptr, Ref(src))
    return dst_ptr[]
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using DashBootstrapComponents

email_input = html_div(
    [
        dbc_label("Email", html_for = "example-email"),
        dbc_input(type = "email", id = "example-email", placeholder = "Enter email"),
        dbc_formtext(
            "Are you on email? You simply have to be these days",
            color = "secondary",
        ),
    ],
    className = "mb-3",
);

password_input = html_div(
    [
        dbc_label("Password", html_for = "example-password"),
        dbc_input(
            type = "password",
            id = "example-password",
            placeholder = "Enter password",
        ),
        dbc_formtext("A password stops mean people taking your stuff", color = "secondary"),
    ],
    className = "mb-3",
);

form = dbc_form([email_input, password_input]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
type SingleVariate <: MCTrait end

type MonteCarloModel{P <: AbstractPathPricer} <: AbstractMonteCarloModel
  pathGenerator::PathGenerator
  pathPricer::P
  sampleAccumulator::RiskStatistics
  isAntitheticVariate::Bool
end

function add_samples!(mcmodel::MonteCarloModel, samples::Int, idx::Int=1)
  # re-init the risk data
  adding_data!(mcmodel.sampleAccumulator, samples)
  for j = 1:samples
    path = get_next!(mcmodel.pathGenerator)
    price = mcmodel.pathPricer(path.value)

    # TODO Control Variate
    if mcmodel.isAntitheticVariate
      path2 = get_antithetic!(mcmodel.pathGenerator)
      price2 = mcmodel.pathPricer(path2.value)
      add_sample!(mcmodel.sampleAccumulator, (price+price2)/2.0, path.weight, idx)
    else
      add_sample!(mcmodel.sampleAccumulator, price, path.weight, idx)
    end
    idx += 1
  end

  return mcmodel
end


text:
# Model tuning

_always_true(xs...) = true

function gridtune(estfun::Function,                          # model estimation function
                  evalfun::Function,                         # model evaluation function
                  params::@compat(Tuple{AbstractString, Any})...;    # parameters to tune
                  ord::Ordering=Forward,                     # ordering of score
                  verbose::Bool=false)                       # whether to display the progress

    np = length(params)
    pnames = [p[1] for p in params]
    pvals = [p[2] for p in params]

    # main loop
    t = 0
    first = true
    local best_score, best_model, best_cfg
    for cf in product(pvals...)
        t += 1
        m = estfun(cf...)
        if m == nothing
            continue
        end
        v = evalfun(m)
        update_best = (first || lt(ord, best_score, v))
        first = false

        if update_best
            best_model = m
            best_cfg = cf
            best_score = v
        end

        if verbose
            print("[")
            for i = 1:np
                print("$(pnames[i])=$(cf[i])")
                if i < np
                    print(", ")
                end
            end
            println("] => $v")
        end
    end

    if first
        error("None of the configs are valid.")
    end
    return (best_model, best_cfg, best_score)
end


text:
# Patrick Keast, MODERATE-DEGREE TETRAHEDRAL QUADRATURE FORMULAS
function _get_gauss_tetdata(n::Int)
    if n == 1
        a = 1. / 4.
        w = 1. / 6.
        xw = [a a a w]
    elseif n == 2
        a = ( 5. + 3. * √(5.) ) / 20.
        b = ( 5. - √(5.) ) / 20.
        w = 1. / 24.
        xw = [a b b w
              b a b w
              b b a w
              b b b w]
    elseif n == 3
        a1 = 1. / 4.
        a2 = 1. / 2.
        b2 = 1. / 6.
        w1 = -2. / 15.
        w2 = 3. / 40.
        xw = [a1 a1 a1 w1
              a2 b2 b2 w2
              b2 a2 b2 w2
              b2 b2 a2 w2
              b2 b2 b2 w2]
    else
        throw(ArgumentError("unsupported order for tetraheder gauss-legendre integration"))
    end
    return xw
end


text:
# this is the log-likelihood for the linear model
# y = x*b+e, with e~N(0,s^2)
# the parameter theta = [b' s]
# this file is used for demonstration of MLE, to compare to OLS
function normal(theta, y, x)
    b = theta[1:end-1][1]
    s = theta[end][1]
    e = (y - x*b)./s
    logdensity = -log.(sqrt.(2.0*pi)) .- 0.5*log(s.^2) .- 0.5*e.*e
end

max_stars_repo_name: Vkamalvanshi/Econometrics

max_issues_repo_name: Vkamalvanshi/Econometrics

max_forks_repo_name: Vkamalvanshi/Econometrics


text:
# additional constructors

function FlintFiniteField(p::Integer; cached::Bool = true)
  @assert is_prime(p)
  k = GF(p, cached=cached)
  return k, k(1)
end

function FlintFiniteField(p::fmpz; cached::Bool = true)
  @assert is_prime(p)
  k = GF(p, cached=cached)
  return k, k(1)
end

GF(p::Integer, k::Int, s::Union{AbstractString,Symbol}=:o; cached::Bool = true) = FlintFiniteField(p, k, s, cached = cached)[1]
GF(p::fmpz, k::Int, s::Union{AbstractString,Symbol}=:o; cached::Bool = true) = FlintFiniteField(p, k, s, cached = cached)[1]

##
## rand for Flint-Finite fields
##
## fq_nmod has no base ring
function rand(R::FqNmodFiniteField)
  #gen is not a generator for the group!
  r = R(0)
  for i=0:degree(R)
    r *= gen(R)
    r += rand(1:characteristic(R))
  end

  return r
end

function rand(R::FqFiniteField)
  #gen is not a generator for the group!
  r = R(0)
  for i=0:degree(R)
    r *= gen(R)
    r += rand(1:characteristic(R))
  end

  return r
end

################################################################################
#
#  Iterators for finite fields
#
################################################################################

# FqNmodFiniteField

function Base.iterate(F::FqNmodFiniteField)
  return zero(F), zeros(UInt, degree(F))
end

function Base.iterate(F::FqNmodFiniteField, st::Vector{UInt})
  done = true
  for j in 1:length(st)
    if st[j] != F.n - 1
      done = false
    end
  end

  if done
    return nothing
  end

  st[1] = st[1] + 1
  for j in 1:(length(st) - 1)
    if st[j] == F.n
      st[j] = 0
      st[j + 1] = st[j + 1] + 1
    end
  end

  d = F()
  ccall((:fq_nmod_init2, libflint), Nothing,
        (Ref{fq_nmod}, Ref{FqNmodFiniteField}), d, F)
  for j in 1:length(st)
         ccall((:nmod_poly_set_coeff_ui, libflint), Nothing,
              (Ref{fq_nmod}, Int, UInt), d, j - 1, st[j])
  end

  return d, st
end

Base.length(F::FqNmodFiniteField) = BigInt(F.n)^degree(F)

Base.IteratorSize(::Type{FqNmodFiniteField}) = Base.HasLength()

Base.eltype(::FqNmodFiniteField) = fq_nmod

# FqFiniteField

function Base.iterate(F::FqFiniteField)
  return zero(F), zeros(FlintZZ, degree(F))
end

function Base.iterate(F::FqFiniteField, st::Vector{fmpz})
  done = true
  for j in 1:length(st)
    if st[j] != characteristic(F) - 1
      done = false
    end
  end

  if done
    return nothing
  end

  st[1] = st[1] + 1
  for j in 1:(length(st) - 1)
    if st[j] == characteristic(F)
      st[j] = 0
      st[j + 1] = st[j + 1] + 1
    end
  end

  d = F()
  ccall((:fq_init2, libflint), Nothing,
        (Ref{fq}, Ref{FqFiniteField}), d, F)
  g = Hecke.Globals.Zx()
  for j in 1:length(st)
    #ccall((:fmpz_poly_set_coeff_fmpz, libflint), (Ref{fmpz_poly}, Int, Ref{fmpz}), g, j - 1, st[j])
    ccall((:fmpz_poly_set_coeff_fmpz, libflint), Nothing,
               (Ref{fq}, Int, Ref{fmpz}), d, j - 1, st[j])
  end

  return d, st
end

Base.eltype(::FqFiniteField) = fq

Base.length(F::FqFiniteField) = BigInt(characteristic(F)^degree(F))

Base.IteratorSize(::Type{FqFiniteField}) = Base.HasLength()

sub!(z::T, x::T, y::T) where {T} = x - y

function _reduce(a::fq_nmod)
  A = parent(a)
  if a.length < 2*degree(A)
    ccall((:fq_nmod_reduce, libflint), Nothing, (Ref{fq_nmod}, Ref{FqNmodFiniteField}), a, A)
  else
    ccall((:nmod_poly_rem, libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{Nothing}, Ref{Nothing}), a, a, pointer_from_objref(A)+6*sizeof(Int) + 2*sizeof(Ptr{Nothing}), pointer_from_objref(A)+sizeof(fmpz))
  end
end

function _reduce(a::fq)
  A = parent(a)
  #if a.length < 2*degree(A)
    ccall((:fq_reduce, libflint), Nothing, (Ref{fq}, Ref{FqFiniteField}), a, A)
  #else
  #  ccall((:fmpz_mod_poly_rem, libflint), Nothing, (Ref{fq}, Ref{fq}, Ref{Nothing}, Ref{Nothing}), a, a, pointer_from_objref(A)+6*sizeof(Int) + 2*sizeof(Ptr{Nothing}), pointer_from_objref(A)+sizeof(fmpz))
  #end
end

function (R::FqFiniteField)(x::fmpz_mod_poly)
  z = R()
  ccall((:fq_set_fmpz_mod_poly, libflint), Nothing, (Ref{Nemo.fq}, Ref{Nemo.fmpz_mod_poly}, Ref{Nemo.FqFiniteField}), z, x, R)
  #ccall((:fq_reduce, libflint), Nothing, (Ref{Nemo.fq}, Ref{Nemo.FqFiniteField}), z, R)
  return z
end

function (R::FqFiniteField)(x::gfp_fmpz_poly)
  z = R()
  ccall((:fq_set_fmpz_mod_poly, libflint), Nothing, (Ref{Nemo.fq}, Ref{Nemo.gfp_fmpz_poly}, Ref{Nemo.FqFiniteField}), z, x, R)
  ccall((:fq_reduce, libflint), Nothing, (Ref{Nemo.fq}, Ref{Nemo.FqFiniteField}), z, R)
  return z
end

#TODO: move elsewhere - and use. There are more calls to nmod_set/reduce
function (A::FqNmodFiniteField)(x::nmod_poly)
  u = A()
  ccall((:fq_nmod_set, libflint), Nothing,
                     (Ref{fq_nmod}, Ref{nmod_poly}, Ref{FqNmodFiniteField}),
                                     u, x, A)
  _reduce(u)
  return u
end

function (A::FqNmodFiniteField)(x::gfp_poly)
  u = A()
  ccall((:fq_nmod_set, libflint), Nothing,
                     (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}),
                                     u, x, A)
  _reduce(u)
  return u
end

function _nf_to_fq!(a::fq_nmod, b::nf_elem, K::FqNmodFiniteField, a_tmp::nmod_poly)
  nf_elem_to_nmod_poly!(a_tmp, b)
  ccall((:fq_nmod_set, libflint), Nothing,
                     (Ref{fq_nmod}, Ref{nmod_poly}, Ref{FqNmodFiniteField}),
                                     a, a_tmp, K)
  _reduce(a)
end

function _nf_to_fq!(a::fq_nmod, b::nf_elem, K::FqNmodFiniteField, a_tmp::gfp_poly)
  nf_elem_to_gfp_poly!(a_tmp, b)
  ccall((:fq_nmod_set, libflint), Nothing,
                     (Ref{fq_nmod}, Ref{gfp_poly}, Ref{FqNmodFiniteField}),
                                     a, a_tmp, K)
  _reduce(a)
end

function _nf_to_fq!(a::fq, b::nf_elem, K::FqFiniteField, a_tmp::gfp_fmpz_poly)
  nf_elem_to_gfp_fmpz_poly!(a_tmp, b)
  ccall((:fq_set, libflint), Nothing,
                     (Ref{fq}, Ref{gfp_fmpz_poly}, Ref{FqFiniteField}),
                                     a, a_tmp, K)
  _reduce(a)
end

function (A::FqNmodFiniteField)(x::gfp_elem)
  @assert characteristic(A) == characteristic(parent(x))
  return A(lift(x))
end

function (A::FqFiniteField)(x::Generic.ResF{fmpz})
  @assert characteristic(A) == characteristic(parent(x))
  return A(lift(x))
end

Nemo.promote_rule(::Type{fq_nmod}, ::Type{gfp_elem}) = fq_nmod

Nemo.promote_rule(::Type{fq}, ::Type{gfp_fmpz_elem}) = fq

################################################################################
#
#  Primitive roots in finite fields
#
################################################################################

# This is not used anywhere
function has_primitive_root_1(K::Nemo.FqNmodFiniteField, m::Int)
  @assert m > 0
  if m == 1
    return K(1)
  end
  if size(K) % m != 1
    return false, K(1)
  end
  if m == 2
    return K(-1)
  end
  fm = factor(m)
  while true
    g = rand(K)
    if iszero(g)
      continue
    end
    if any(x -> isone(g^div(m, x)), keys(fm))
      continue
    end
    return true, g^div(size(K)-1, m)
  end
end

## Minpoly/ Charpoly

function minpoly(a::fq_nmod)
  return minpoly(PolynomialRing(GF(Int(characteristic(parent(a))), cached = false), cached = false)[1], a)
end

function minpoly(Rx::GFPPolyRing, a::fq_nmod)
  c = [a]
  fa = frobenius(a)
  while !(fa in c)
    push!(c, fa)
    fa = frobenius(fa)
  end
  St = PolynomialRing(parent(a), cached = false)[1]
  f = prod([gen(St) - x for x = c], init = one(St))
  g = Rx()
  for i = 0:degree(f)
    setcoeff!(g, i, coeff(coeff(f, i), 0))
  end
  return g
end

function charpoly(a::fq_nmod)
  return charpoly(PolynomialRing(GF(Int(characteristic(parent(a))), cached = false), cached = false)[1], a)
end

function charpoly(Rx::GFPPolyRing, a::fq_nmod)
  g = minpoly(Rx, a)
  return g^div(degree(parent(a)), degree(g))
end

function minpoly(a::fq)
  return minpoly(PolynomialRing(GF(characteristic(parent(a)), cached = false), cached = false)[1], a)
end

function minpoly(Rx::GFPFmpzPolyRing, a::fq)
  c = [a]
  fa = frobenius(a)
  while !(fa in c)
    push!(c, fa)
    fa = frobenius(fa)
  end
  St = PolynomialRing(parent(a), cached = false)[1]
  f = prod([gen(St) - x for x = c])
  g = Rx()
  for i = 0:degree(f)
    setcoeff!(g, i, coeff(coeff(f, i), 0))
  end
  return g
end

################################################################################
#
#  Generators multiplicative group
#
################################################################################

function primitive_element(F::T; n_quo::Int = -1) where T <: Union{FqFiniteField, FqNmodFiniteField, GaloisField, Nemo.GaloisFmpzField}
  n = size(F)-1
  k = fmpz(1)
  if n_quo != -1
    if !divisible(n, n_quo)
      return F(1)
    end
    n, k = ppio(n, fmpz(n_quo))
  end
  primefactors_n = collect(keys(factor(n).fac))

  x = rand(F)^k
  while iszero(x)
    x = rand(F)^k
  end
  while true
    found = true
    for l in primefactors_n
      if isone(x^div(n, l))
        found = false
        break
      end
    end
    if found
      break
    end
    x = rand(F)^k
    while iszero(x)
      x = rand(F)^k
    end
  end
  return x
end

function unit_group(F::T; n_quo::Int = -1) where T <: FinField

  g = primitive_element(F, n_quo = n_quo)
  k = size(F) - 1
  inv = fmpz(1)
  npart = fmpz(k)
  if n_quo != -1
    k = fmpz(n_quo)
    npart, nnpart = ppio(size(F)-1, k)
    inv = invmod(nnpart, npart)
  end

  G = abelian_group(Int[k])
  ex = div(size(F)-1, npart)
  function disc_log(x)
    @assert typeof(x) == elem_type(F)
    iszero(x) && error("Not invertible!")
    y = x^ex
    isone(y) && return G[0]
    if k < 20
      c = 1
      el = deepcopy(g)
      while el != y
        c += 1
        if c > npart
          error("Something went wrong!")
        end
        el *= g
      end
      return G([mod(c*inv, k)])
    else
      return G([mod(inv*disc_log_bs_gs(g, y, npart), k)])
    end
  end
  mG = FiniteFieldMultGrpMap{T, elem_type(F)}(G, F, g, disc_log)
  return G, mG
end

################################################################################
#
#  Missing ad hoc operations
#
################################################################################

function (R::FqFiniteField)(x::Nemo.gfp_fmpz_elem)
  return R(lift(x))
end

function *(a::Nemo.fq, b::Nemo.gfp_fmpz_elem)
  return a * parent(a)(b)
end

function *(a::Nemo.gfp_fmpz_elem, b::Nemo.fq)
  return parent(b)(a) * b
end
function Hecke.preimage(phi::Nemo.FinFieldMorphism, x::FinFieldElem)
  return preimage_map(phi)(x)
end

function (R::Nemo.NmodRing)(a::Nemo.gfp_elem)
  @assert modulus(R) == characteristic(parent(a))
  return R(data(a))
end

function (k::Nemo.FqNmodFiniteField)(a::fmpq)
  return k(numerator(a))//k(denominator(a))
end

function (k::Nemo.GaloisFmpzField)(a::fmpq)
  return k(numerator(a))//k(denominator(a))
end

function (k::Nemo.FqFiniteField)(a::fmpq)
  return k(numerator(a))//k(denominator(a))
end

(F::Nemo.FqNmodFiniteField)(a::Nemo.nmod) = F(a.data)

#to avoid using embed - which is (more me) still broken..
# it accumulates fields until the machine dies
function find_morphism(k::Nemo.NmodRing, K::FqNmodFiniteField)
  return x->K(x.data)
end

function find_morphism(k::FqNmodFiniteField, K::FqNmodFiniteField)
   if degree(k) > 1
    phi = Nemo.find_morphism(k, K) #avoids embed - which stores the info
  else
    phi = MapFromFunc(x->K((coeff(x, 0))), y->k((coeff(y, 0))), k, K)
  end
  return phi
end

function frobenius_matrix(F::FqNmodFiniteField, n::Int=1)
  a = frobenius(gen(F), n)
  k = quo(ZZ, Int(characteristic(F)))[1]
  m = zero_matrix(k, degree(F), degree(F))
  ccall((:fq_nmod_embed_composition_matrix_sub, libflint), Nothing, (Ref{nmod_mat}, Ref{fq_nmod}, Ref{FqNmodFiniteField}, Clong), m, a, F, degree(F))
  ccall((:nmod_mat_transpose, libflint), Nothing, (Ref{nmod_mat}, Ref{nmod_mat}), m, m)
  return m
end

mutable struct VeryBad
  entries::Ptr{Nothing}
  r::Clong
  c::Clong
  rows::Ptr{Nothing}
  n::Culong
  ninv::Culong
  norm::Culong

  function VeryBad(n, ninv, norm)
    r = new()
    r.n = n
    r.ninv = ninv
    r.norm = norm
    r.r = 1
    r.rr = [reinterpret(Ptr{Nothing}, 0)]
    r.rows = Base.unsafe_convert(Ptr{Nothing}, r.rr)
    return r
  end

  rr::Vector{Ptr{Nothing}}
end

function VeryBad!(V::VeryBad, a::fq_nmod)
  V.c = a.length
  V.entries = a.coeffs
  V.rr[1] = a.coeffs
#  V.rows = Base.unsafe_convert(Ptr{Nothing}, [a.coeffs])
end

function clear!(V::VeryBad)
  V.entries = reinterpret(Ptr{Nothing}, 0)
#  V.rows = reinterpret(Ptr{Nothing}, 0)
end

struct FrobeniusCtx
  m::nmod_mat
  fa::VeryBad
  fb::VeryBad
  K::FqNmodFiniteField
  i::Int

  function FrobeniusCtx(K::FqNmodFiniteField, i::Int = 1)
    m = frobenius_matrix(K, i)
    return new(m, VeryBad(m.n, m.ninv, m.norm), VeryBad(m.n, m.ninv, m.norm), K, i)
  end
end

function show(io::IO, F::FrobeniusCtx)
  println(io, "$(F.i)-th Frobenius data for $(F.K)")
end

function apply!(b::fq_nmod, a::fq_nmod, F::FrobeniusCtx)
  n = degree(parent(a))
  ccall((:nmod_poly_fit_length, libflint), Nothing, (Ref{fq_nmod}, Clong), b, n)
  VeryBad!(F.fa, a)
  VeryBad!(F.fb, b)
  ccall((:nmod_mat_mul, libflint), Nothing, (Ref{VeryBad}, Ref{VeryBad}, Ref{nmod_mat}), F.fb, F.fa, F.m)
  b.length = n
  clear!(F.fa)
  clear!(F.fb)
  ccall((:_nmod_poly_normalise, libflint), Nothing, (Ref{fq_nmod}, ), b)
  return b
end

function frobenius!(a::fq_nmod, b::fq_nmod, i::Int = 1)
    ccall((:fq_nmod_frobenius, libflint), Nothing,
         (Ref{fq_nmod}, Ref{fq_nmod}, Int, Ref{FqNmodFiniteField}),
                                               a, b, i, a.parent)
end


text:
function forward_pass!(policy::PolicyData, problem::ProblemData, data::SolverData;
    line_search=:armijo,
    min_step_size=1.0e-5,
    c1=1.0e-4,
    c2=0.9,
    max_iterations=25,
    verbose=false)

    # reset solver status
    data.status[1] = false

    # previous cost
    J_prev = data.objective[1]

    # gradient of Lagrangian
    lagrangian_gradient!(data, policy, problem)

    if line_search == :armijo
        trajectory_sensitivities(problem, policy, data)
        delta_grad_product = data.gradient' * problem.trajectory
    else
        delta_grad_product = 0.0
    end

    # line search with rollout
    data.step_size[1] = 1.0
    iteration = 1
    while data.step_size[1] >= min_step_size
        iteration > max_iterations && (verbose && (@warn "forward pass failure"), break)

        J = Inf
        #TODO: remove try-catch
        # try
        rollout!(policy, problem,
            step_size=data.step_size[1])
        J = cost!(data, problem,
            mode=:current)[1]
        # catch
        #     if verbose
        #         @warn "rollout failure"
        #         @show norm(data.gradient)
        #     end
        # end
        if (J <= J_prev + c1 * data.step_size[1] * delta_grad_product)
            # update nominal
            update_nominal_trajectory!(problem)
            data.objective[1] = J
            data.status[1] = true
            break
        else
            data.step_size[1] *= 0.5
            iteration += 1
        end
    end
    data.step_size[1] < min_step_size && (verbose && (@warn "line search failure"))
end


text:
function read(filename)
    f = open(filename, "r")
    # head = readline(f)
    # n, d, l = split(head, " ")
    # n = parse(Int, n)
    # d = parse(Int, d)
    # l = parse(Int, l)

    row_ft = zeros(Int, 0)
    col_ft = zeros(Int, 0)
    val_ft = zeros(Float64, 0)

    row_lbl = zeros(Int, 0)
    col_lbl = zeros(Int, 0)
    val_lbl = zeros(Float64, 0)

    cc = 1
    for line in eachline(f)
        # println(cc)
        tmp = split(line, " ")
        if length(tmp) <= 1
            println("weird at line: ", cc)
            continue
        end

        if tmp[1] != ""
            lbls = split(tmp[1], ",")
        else
            lbls = []
        end

        lbls = [ parse(Int, x) for x in lbls ]

        for lbl in lbls
            push!(row_lbl, lbl)
            push!(col_lbl, cc)
            push!(val_lbl, 1)
        end

        for i = 2:length(tmp)
            if tmp[i] == ""
                continue
            end
            idx, val = split(tmp[i], ":")
            idx = parse(Int, idx)
            val = parse(Float64, val)
            push!(col_ft, idx )
            push!(val_ft, val)
            push!(row_ft, cc)
        end
        cc += 1
    end
    close(f)

    n = cc-1

    ft_mat = sparse( row_ft, col_ft, val_ft )
    lbl_mat = sparse( row_lbl, col_lbl, val_lbl )

    return ft_mat, lbl_mat

end


text:
############################################################################################
"""
$(TYPEDEF)

Stores in diagnostics.

# Fields
$(TYPEDFIELDS)
"""
struct AcceptStatistic
    "Acceptance rate"
    rate::Float64
    "Step accepted or rejected"
    accepted::Bool
    function AcceptStatistic(rate::T, accepted::Bool) where {T<:AbstractFloat}
        ArgCheck.@argcheck 0.0 <= rate <= 1.0 "Acceptance rate is out of boundaries"
        return new(rate, accepted)
    end
end
function AcceptStatistic(_rng::Random.AbstractRNG, ℓacceptrate::F) where {F<:AbstractFloat}
    rate = exp(min(zero(F), ℓacceptrate))
    accepted = rand(_rng) < rate
    return AcceptStatistic(rate, accepted)
end

############################################################################################
# Export
export AcceptStatistic


text:
using ITensors,
      HDF5,
      Test

@testset "HDF5 Read and Write" begin

  i = Index(2,"i")
  j = Index(3,"j")
  k = Index(4,"k")

  @testset "TagSet" begin
    ts = TagSet("A,Site,n=2")
    fo = h5open("data.h5","w")
    write(fo,"tags",ts)
    close(fo)

    fi = h5open("data.h5","r")
    rts = read(fi,"tags",TagSet)
    close(fi)
    @test rts == ts
  end

  @testset "Index" begin
    i = Index(3,"Site,S=1")
    fo = h5open("data.h5","w")
    write(fo,"index",i)
    close(fo)

    fi = h5open("data.h5","r")
    ri = read(fi,"index",Index)
    close(fi)
    @test ri == i

    # primed Index
    i = Index(3,"Site,S=1")
    i = prime(i,2)
    fo = h5open("data.h5","w")
    write(fo,"index",i)
    close(fo)

    fi = h5open("data.h5","r")
    ri = read(fi,"index",Index)
    close(fi)
    @test ri == i
  end

  @testset "IndexSet" begin
    is = IndexSet(i,j,k)

    fo = h5open("data.h5","w")
    write(fo,"inds",is)
    close(fo)

    fi = h5open("data.h5","r")
    ris = read(fi,"inds",IndexSet)
    close(fi)
    @test ris == is
  end

  @testset "ITensor" begin

    # default constructed case
    T = ITensor()

    fo = h5open("data.h5","w")
    write(fo,"defaultT",T)
    close(fo)

    fi = h5open("data.h5","r")
    rT = read(fi,"defaultT",ITensor)
    close(fi)
    @test typeof(store(T)) == typeof(store(ITensor()))

    # real case
    T = randomITensor(i,j,k)

    fo = h5open("data.h5","w")
    write(fo,"T",T)
    close(fo)

    fi = h5open("data.h5","r")
    rT = read(fi,"T",ITensor)
    close(fi)
    @test norm(rT-T)/norm(T) < 1E-10

    # complex case
    T = randomITensor(ComplexF64,i,j,k)

    fo = h5open("data.h5","w")
    write(fo,"complexT",T)
    close(fo)

    fi = h5open("data.h5","r")
    rT = read(fi,"complexT",ITensor)
    close(fi)
    @test norm(rT-T)/norm(T) < 1E-10
  end

  #
  # Clean up the test hdf5 file
  #
  rm("data.h5",force=true)

end

nothing
nothing


text:
"""
arquivo.rsd

SINTAXE:

definir: # SINTAXE=> definir:
    material: "CONCRETO 20MPA" # SINTAXE=> material: {ID Material} ...
    no: 0.0 0.0 0.0 "A" # SINTAXE=> no: {x} {y} {z} {ID Nó}
    no: 0.0 1.0 0.0 "B" # SINTAXE=> no: {x} {y} {z} {ID Nó}
    elemento: "A" "B" "PILAR P1" # elemento: {ID Nó inicial} {ID Nó final} {ID Elemento}

com elemento "PILAR P1": # SINTAXE=> com elemento {ID Elemento}:
    propriedades: material="CONCRETO 20MPA" # SINTAXE=> propriedades: material={ID Material}  ...

com no "A": # SINTAXE=> com no {ID Nó}:
    carregar: 0.0 -12000.00 0.0 # SINTAXE=> carregar: {x} {y} {z} ...
    restringir: xyz # SINTAXE=> restringir: {[z,y,z]} ...

analisar # SINTAXE=> analisar
relatorio # SINTAXE=> relatorio
dimensionar # SINTAXE=> dimensionar

"""
module RSDLang
   function load()

   end

   function tokenize()

   end

   function evaluate()

   end

   function intepret()

   end
end

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
push!(LOAD_PATH, "./")
using CAS_QMDP

jld_full_filename = "../policies/fine_maxdist_2km_colli_150m_pen_action_5.0_pen_conflict_50.0_pen_closeness_0.1.jld"
alphas = load_alphas(jld_full_filename)
csv_filename = "fine_maxdist_2km_colli_150m_pen_action_5.0_pen_conflict_50.0_pen_closeness_0.1"
write_alphas_to_csv(alphas, csv_filename)

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
module DistributedData

using Distributed
using Serialization

include("structs.jl")
export Dinfo

include("base.jl")
export save_at,
    get_from,
    get_val_from,
    remove_from,
    scatter_array,
    unscatter,
    dexec,
    dtransform,
    dmapreduce,
    dmap,
    dpmap,
    gather_array,
    tmp_symbol,
    @remote

include("io.jl")
export dstore, dload, dunlink

include("tools.jl")
export dcopy,
    dselect,
    dapply_cols,
    dapply_rows,
    dstat,
    dstat_buckets,
    dcount,
    dcount_buckets,
    dscale,
    dmedian,
    dmedian_buckets

end # module