text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

# Newton method

# What is this guarded_mid for? Shouldn't it be checking if f(m)==0?
doc"""Returns the midpoint of the interval x, slightly shifted in case
it is zero"""
function guarded_mid{T}(x::Interval{T})
    m = mid(x)
    if m == zero(T)                     # midpoint exactly 0
        α = convert(T, 0.46875)      # close to 0.5, but exactly representable as a floating point
        m = α*x.lo + (one(T)-α)*x.hi   # displace to another point in the interval
    end

    @assert m ∈ x

    m
end

function N{T}(f::Function, x::Interval{T}, deriv::Interval{T})
    m = guarded_mid(x)
    m = Interval(m)

    m - (f(m) / deriv)
end

doc"""If a root is known to be inside an interval,
`newton_refine` iterates the interval Newton method until that root is found."""
function newton_refine{T}(f::Function, f_prime::Function, x::Interval{T};
                          tolerance=eps(T), debug=false)

    debug && (print("Entering newton_refine:"); @show x)

    while diam(x) > tolerance  # avoid problem with tiny floating-point numbers if 0 is a root
        deriv = f_prime(x)
        Nx = N(f, x, deriv)

        debug && @show(x, Nx)

        Nx = Nx ∩ x
        Nx == x && break

        x = Nx
    end

    debug && @show "Refined root", x

    return [Root(x, :unique)]
end

doc"""`newton` performs the interval Newton method on the given function `f`
with its optional derivative `f_prime` and initial interval `x`.
Optional keyword arguments give the `tolerance`, `maxlevel` at which to stop
subdividing, and a `debug` boolean argument that prints out diagnostic information."""

function newton{T}(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
                   tolerance=eps(T), debug=false, maxlevel=30)

    debug && (print("Entering newton:"); @show(level); @show(x))

    isempty(x) && return Root{T}[]  #[(x, :none)]

    z = zero(x.lo)
    tolerance = abs(tolerance)

    # Tolerance reached or maximum level of bisection
    (diam(x) < tolerance || level >= maxlevel) && return [Root(x, :unknown)]

    deriv = f_prime(x)
    debug && @show(deriv)

    if !(z in deriv)
        Nx = N(f, x, deriv)
        debug && @show(Nx, Nx ⊆ x, Nx ∩ x)

        isempty(Nx ∩ x) && return Root{T}[]

        if Nx ⊆ x
            debug && (print("Refining "); @show(x))
            return newton_refine(f, f_prime, Nx, tolerance=tolerance, debug=debug)
        end

        m = guarded_mid(x)

        debug && @show(x,m)

        # bisect:
        roots = vcat(
            newton(f, f_prime, Interval(x.lo, m), level+1,
                   tolerance=tolerance, debug=debug, maxlevel=maxlevel),
            newton(f, f_prime, Interval(m, x.hi), level+1,
                   tolerance=tolerance, debug=debug, maxlevel=maxlevel)
            # note the nextfloat here to prevent repetition
            )

        debug && @show(roots)

    else  # 0 in deriv; this does extended interval division by hand

        y1 = Interval(deriv.lo, -z)
        y2 = Interval(z, deriv.hi)

        if debug
            @show (y1, y2)
            @show N(f, x, y1)
            @show N(f, x, y2)
        end

        y1 = N(f, x, y1) ∩ x
        y2 = N(f, x, y2) ∩ x

        debug && @show(y1, y2)

        roots = vcat(
            newton(f, f_prime, y1, level+1,
                   tolerance=tolerance, debug=debug, maxlevel=maxlevel),
            newton(f, f_prime, y2, level+1,
                   tolerance=tolerance, debug=debug, maxlevel=maxlevel)
            )

        debug && @show(roots)

    end

    roots = clean_roots(f, roots)

    return roots

end

# use automatic differentiation if no derivative function given:
newton{T}(f::Function, x::Interval{T};  args...) =
    newton(f, x->D(f,x), x; args...)

# newton for vector of intervals:
newton{T}(f::Function, f_prime::Function, xx::Vector{Interval{T}}; args...) =
    vcat([newton(f, f_prime, @interval(x); args...) for x in xx]...)

newton{T}(f::Function,  xx::Vector{Interval{T}}, level; args...) =
    newton(f, x->D(f,x), xx, 0, args...)

newton{T}(f::Function,  xx::Vector{Root{T}}; args...) =
    newton(f, x->D(f,x), [x.interval for x in xx], args...)


text:
"""
# References
- Controller
[1] G. P. Falconi and F. Holzapfel,
“Adaptive Fault Tolerant Control Allocation for a Hexacopter System,”
Proc. Am. Control Conf., vol. 2016-July, pp. 6760–6766, 2016.
- Reference model (e.g., xd, vd, ad, ad_dot, ad_ddot)
[2] S. J. Su, Y. Y. Zhu, H. R. Wang, and C. Yun, “A Method to Construct a Reference Model for Model Reference Adaptive Control,” Adv. Mech. Eng., vol. 11, no. 11, pp. 1–9, 2019.
"""
struct AdaptiveAllocator <: AbstractEnv
    B
    B_pinv
    γ
    function AdaptiveAllocator(B, γ=1e-2)
        B_pinv = pinv(B)
        new(B, B_pinv, γ)
    end
end

function State(allocator::AdaptiveAllocator)
    return function (Θ̂=zeros(6, 4))
        Θ̂
    end
end

function Dynamics!(allocator::AdaptiveAllocator)
    @Loggable function dynamics!(dΘ̂, Θ̂, p, t; Θ̂̇)
        @log Θ̂
        dΘ̂ .= Θ̂̇
    end
end

"""
Several functions are exported from `utils.jl`, e.g., T_u_inv(T).
"""
function Command(allocator::AdaptiveAllocator)
    @unpack B, γ = allocator
    return function (νd, e, zB, T, Θ̂, R, J,
                     Ap, Bp, P, Kp, Kt, Kω)
        P̄ = [          P         zeros(6, 6);
             zeros(6, 6) 0.5*Matrix(I, 6, 6)]
        θ_ėp = Bp
        θ_u̇1 = Kp * θ_ėp
        θ_ėt = θ_u̇1
        θ_ëp = Ap * θ_ėp + Bp * θ_ėt
        θ_ü1 = Kp * θ_ëp
        θ_u̇2 = T_u_inv(T) * R * (2*Bp' * P * θ_ėp
                                 + θ_ü1 + Kt * θ_ėt)
        θ_ω̇d = [1 0 0;
                0 1 0;
                0 0 0] * θ_u̇2
        B̄ = [-(θ_ėp*zB) zeros(6, 3);
             -(θ_u̇1*zB) zeros(3, 3);
             -(θ_ω̇d*zB) inv(J)]
        Θ̂̇ = γ * Proj_R(
                       Θ̂,
                       (νd * e' * P̄ * B̄ * B)'
                      )
    end
end

function Proj_R(C, Y)
    proj_R = zeros(size(Y))
    for i in 1:size(Y)[1]
        ci = C[i, :]
        yi = Y[i, :]
        proj_R[i, :] = Proj(ci, yi)
    end
    proj_R
end

function Proj(θ, y; ϵ=1e-2, θ_max=1e5)
    f = ((1+ϵ) * norm(θ)^2 - θ_max^2) / (ϵ * θ_max^2)
    ∇f = (2*(1+ϵ) / (ϵ * θ_max^2)) * θ
    proj = nothing
    if f > 0 && ∇f' * y > 0
        ∇f_unit = ∇f / norm(∇f)
        proj = y - dot(∇f, y) * ∇f_unit
    else
        proj = y
    end
    proj
end

function (allocator::AdaptiveAllocator)(ν, Θ̂)
    @unpack B_pinv = allocator
    u = (B_pinv + Θ̂) * ν
end

ext: jl


text:
imagedir = joinpath(dirname(@__FILE__), "..", "images")
if !isdir(imagedir)
    mkdir(imagedir)
end

REPO_URL = "https://github.com/timholy/TestImages.jl/blob/gh-pages/images/"

stdfiles = [
    "cameraman.tif" ,
    "house.tif" ,
    "jetplane.tif" ,
    "lake_gray.tif" ,
    "lake_color.tif" ,
    "lena_color_256.tif" ,
    "lena_color_512.tif" ,
    "lena_gray_256.tif" ,
    "lena_gray_512.tif" ,
    "livingroom.tif" ,
    "mandril_color.tif" ,
    "mandril_gray.tif" ,
    "mandrill.tiff" ,
    "peppers_color.tif" ,
    "peppers_gray.tif" ,
    "pirate.tif" ,
    "walkbridge.tif" ,
    "woman_blonde.tif" ,
    "woman_darkhair.tif" ,
]

@info "Downloading standard test images"
for f in stdfiles
    fn = joinpath(imagedir, f)
    if !isfile(fn)
        @info "Downloading $fn"
        download(REPO_URL*f*"?raw=true", joinpath(imagedir, f))
    end
end

@info "Download Completed."


text:
module CVHighGUI

export imshow, waitKey, destroyWIndow, destroyAllWindows, namedWindow,
    createTrackbar, getTrackbarPos

using LibOpenCV
using CVCore
using Cxx

libopencv_highgui = LibOpenCV.find_library_e("libopencv_highgui")
try
    Libdl.dlopen(libopencv_highgui, Libdl.RTLD_GLOBAL)
catch e
    warn("You might need to set DYLD_LIBRARY_PATH to load dependencies proeprty.")
    rethrow(e)
end

cxx"""
#include <opencv2/highgui.hpp>
"""

for name in [
    #Flags for cv::namedWindow
    :WINDOW_NORMAL,
    :WINDOW_AUTOSIZE,
    :WINDOW_OPENGL,
    :WINDOW_FULLSCREEN,
    :WINDOW_FREERATIO,
    :WINDOW_KEEPRATIO,
    ]
    ex = Expr(:macrocall, Symbol("@icxx_str"), string("cv::", name, ";"))
    @eval global const $name = $ex
end

imshow(winname::AbstractString, mat::AbstractCvMat) =
    icxx"cv::imshow($(pointer(winname)), $(mat.handle));"
imshow(winname::AbstractString, arr::Array) = imshow(winname, Mat(arr))
imshow(winname::AbstractString, expr::MatExpr) = imshow(winname, Mat(expr))

destroyWindow(winname::AbstractString) =
    icxx"cv::destroyWindow($(pointer(winname)));"
destroyAllWindows() = icxx"cv::destroyAllWindows();"

waitKey(delay) = icxx"cv::waitKey($delay);"
waitKey(;delay::Int=0) = waitKey(delay)

function createTrackbar(trackbarname::String, winname::String,
    value::Ptr{Cint}, count, onChange::Ptr{Void}=C_NULL,
    userdata::Ptr{Void}=C_NULL)
    icxx"""
    cv::createTrackbar($(pointer(trackbarname)), $(pointer(winname)),
        $value, $count, (void(*)(int,void*))$onChange, $userdata);
    """
end
function createTrackbar(trackbarname::String, winname::String,
    value::Vector{Cint}, count;
    onChange::Ptr{Void}=C_NULL, userdata::Ptr{Void}=C_NULL)
    createTrackbar(trackbarname, winname, pointer(value), count,
        onChange, userdata)
end
function createTrackbar(trackbarname::String, winname::String,
    value::Number, count;
    onChange::Ptr{Void}=C_NULL, userdata::Ptr{Void}=C_NULL)
    v = Cint[value]
    createTrackbar(trackbarname, winname, pointer(v), count,
        onChange, userdata)
end

function getTrackbarPos(trackbarname::String, winname::String)
    icxx"cv::getTrackbarPos($(pointer(trackbarname)), $(pointer(winname)));"
end

function namedWindow(winname::String, flags=WINDOW_AUTOSIZE)
    icxx"cv::namedWindow($(pointer(winname)), $flags);"
end

end # module


text:
import Random
using Random

@testset "sort_biosystem_compounds" begin

    compounds = readcompounds("data/compound")

    @testset "sortkey: exact_mass" begin

        biosystem_compounds = readids(["C00068","C00069","C00379","C00380","C00381"])
        sortkey = :exact_mass
        zero_mass_behavior = "end"

        expected_exact_mass = [("C00380",111.0433),
                               ("C00379",152.0685),
                               ("C00068",425.045),
                               ("C00381",0.0),
                               ("C00069",0.0)]

        Random.seed!(1234);
        @test expected_exact_mass == BioXP.sort_biosystem_compounds(compounds,
                                                    biosystem_compounds,
                                                    sortkey,
                                                    zero_mass_behavior)

        for i in 1:1000
            @test expected_exact_mass[1:3] == BioXP.sort_biosystem_compounds(compounds,
                                                    biosystem_compounds,
                                                    sortkey,
                                                    zero_mass_behavior)[1:3]

        end

    end

    @testset "sortkey: mol_weight" begin

        biosystem_compounds = readids(["C00068","C00069","C00379","C00380","C00381"])
        sortkey = :mol_weight
        zero_mass_behavior = "end"

        expected_mol_weight = [("C00380",111.102),
                               ("C00379",152.1458),
                               ("C00068",425.3144),
                               ("C00381",0.0),
                               ("C00069",0.0)]

        Random.seed!(1234);
        @test expected_mol_weight == BioXP.sort_biosystem_compounds(compounds,
                                                    biosystem_compounds,
                                                    sortkey,
                                                    zero_mass_behavior)

        for i in 1:1000
            @test expected_mol_weight[1:3] == BioXP.sort_biosystem_compounds(compounds,
                                                    biosystem_compounds,
                                                    sortkey,
                                                    zero_mass_behavior)[1:3]

        end
    end

    @testset "zero_mass_behavior: random" begin

        biosystem_compounds = readids(["C00068","C00069","C00379","C00380","C00381"])
        sortkey = :exact_mass
        zero_mass_behavior = "random"

        expected_exact_mass = [("C00069",0.0),
                            ("C00380",111.0433),
                            ("C00379",152.0685),
                            ("C00381",0.0),
                            ("C00068",425.045)]

        expected_order = [("C00380",111.0433),
                          ("C00379",152.0685),
                          ("C00068",425.045)]

        Random.seed!(1234);
        # @test expected_exact_mass == BioXP.sort_biosystem_compounds(compounds,
        #                                             biosystem_compounds,
        #                                             sortkey,
        #                                             zero_mass_behavior)

        for i in 1:1000
            sorted_cpd_masses = BioXP.sort_biosystem_compounds(compounds,
                                            biosystem_compounds,
                                            sortkey,
                                            zero_mass_behavior)

            ## Test that the compounds with masses are always in order of inc mass
            pos = 0
            for cm in expected_order
                newpos = findfirst(isequal(cm),sorted_cpd_masses)
                @test newpos>pos
                pos = newpos
            end

        end

    end

    @testset "flip probability" begin

        beta = 20

        @test .5 == BioXP.getflipprobability(("C00069",0.0),("C00066",0.0),beta)
        @test .5 == BioXP.getflipprobability(("C00069",102.0),("C00066",0.0),beta)
        @test .5 == BioXP.getflipprobability(("C00069",0.000),("C00066",6663.0),beta)
        @test .12857 ≈ BioXP.getflipprobability(("C00379",152.0685),("C00380",111.0433),beta) atol=0.00001
        @test 1.0 == BioXP.getflipprobability(("C00380",111.0433),("C00379",152.0685),beta)

    end

    @testset "randomizecompounds" begin
        biosystem_compounds = readids(["C00068","C00069","C00379","C00380","C00381"])
        n_runs = 3

        # println(randomizecompounds(biosystem_compounds,compounds,n_runs))
        runs = randomizecompounds(biosystem_compounds,compounds,n_runs)
        for i in runs
            @test length(biosystem_compounds) == length(i)
        end
    end

end

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
function ψtcNewtonSolve!(F, ∇ₓF, x, nrm, τstop)
  # TODO x_hist = x
    @printf "Solving F(x) = 0 (Using Pseudotransient Continuation)\n"
    nit_NS = 200
    i = 0
    s = copy(x)
    FNorm = nrm(F(x))
    FNorm₀ = copy(FNorm)
    xNorm = nrm(x)
    τ = xNorm / FNorm
    sNorm = nrm(s)
  # evaluate F at the initial iterate
  # compute the stop tolerance
  # main iteration loop
    @printf "  │ i = %d: |F(x)| = %.2g\n" i FNorm
    while (i == 0) || ((τ < τstop) && (i < nit_NS))
        i += 1
        @printf "  │ i = %d: " i
    # compute the ψtc step
        τ = xNorm / FNorm
        s .= (I / τ - ∇ₓF(x)) \ F(x)
        x .= x + s # ψtc step
        FNorm = nrm(F(x))
        xNorm = nrm(x)
        sNorm = nrm(s)
        @printf "|F(x)| = %.2g" FNorm
        @printf "; |x| = %.2g" xNorm
        @printf ", |dx|/|x| = %.2g\n" sNorm / xNorm
    end
  # on failure, set the error flag
    if (nrm(x) / nrm(F(x)) < τstop)
        @printf "  └────> Newton has reached max iterations, "
        @printf "|x|/|F(x)| = %.2g years\n" ((xNorm / FNorm) / (365 * 24 * 60 * 60))
    else
        @printf "  └────> Newton has converged, "
        @printf "|x|/|F(x)| = %.2g years\n" ((xNorm / FNorm) / (365 * 24 * 60 * 60))
    end
  # TODO return sol, x_hist
    return x
end

export ψtcNewtonSolve!


text:
@testset "891.sum-of-subsequence-widths.jl" begin
    @test sum_subseq_widths([2, 1, 3]) == 6
    @test sum_subseq_widths([2, 1, 1]) == 3
    @test sum_subseq_widths([1, 2, 3, 4, 5]) == 72
    @test sum_subseq_widths([3858, 3747, 3753, 2643, 4625, 1331, 4903, 1403, 1594, 3603,
                             728, 2472, 2769, 3113, 2691, 1587, 313, 4913, 2865, 530, 1190,
                             2843, 238, 2244, 1787, 2598, 4863, 1054, 2984, 4907, 1812,
                             1398, 3797, 4548, 2535, 4286, 1937, 3602, 3150, 3808, 1731,
                             4794, 276, 4195, 1527, 1809, 1001, 4429, 2268, 468, 1631, 1710,
                             3279, 1841, 3824, 3636, 730, 3491, 3084, 134, 1188, 3174, 4655,
                             2662, 1277, 185, 2088, 4613, 3557, 1357, 3489, 3738, 3620,
                             2257, 3295, 1194, 109, 4355, 3759, 2552, 2537, 760, 4296, 2522,
                             1962, 648, 828, 1295, 2618, 1425, 4671, 522, 390, 1620, 1028,
                             1608, 2734, 3383, 1559, 2249, 2551, 2495, 3181, 2909, 3430,
                             2955, 3050, 72, 1949, 4984, 1732, 4402, 546, 3819, 1915, 1817,
                             2128, 3175, 153, 273, 2993, 4194, 40, 827, 4046, 2265, 3181,
                             310, 514, 4146, 4132, 4060, 3363, 4509, 4847, 3249, 4626, 2907,
                             1685, 3153, 2730, 3044, 2724, 4658, 1651, 3815, 2385, 2319,
                             3036, 978, 743, 991, 3147, 3610, 2915, 3520, 904, 2292, 4321,
                             2591, 3806, 534, 4872, 2116, 4170, 3283, 4682, 2942, 2117,
                             2871, 3972, 3711, 177, 1525, 3658, 2317, 4027, 1148, 4250,
                             3592, 1252, 1715, 3884, 1157, 3405, 3818, 3603, 4086, 4879, 1,
                             2649, 1528, 4879, 3988, 3104, 779, 2050, 2725, 11, 1452, 4513,
                             3789, 3592, 2847, 7, 3911, 589, 2944, 2805, 2932, 2919, 3812,
                             2372, 1778, 4858, 1739, 1189, 98, 487, 3583, 644, 2437, 4889,
                             1717, 1726, 2459, 1087, 3171, 851, 4395, 1939, 873, 2579, 2741,
                             3622, 2719, 3925, 1700, 2801, 2428, 2718, 2029, 574, 758, 3753,
                             3244, 4313, 3786, 4101, 1423, 4213, 3525, 2152, 1289, 4793,
                             590, 717, 885, 4823, 1506, 1361, 1093, 4215, 173, 1725, 1509,
                             3656, 1520, 2294, 697, 657, 300, 4349, 2773, 745, 3788, 2018,
                             583, 1620, 2620, 802, 3782, 4603, 3455, 2942, 3596, 2536, 45,
                             1695, 4866, 2414, 2748, 1299, 3229, 1330, 1998, 1544, 934,
                             3657, 1918, 401, 235, 213, 612, 4160, 1758, 3547, 1321, 2991,
                             2829, 2690, 1648, 2552, 328, 2210, 4397, 3853, 4861, 3003,
                             4148, 4434, 2863, 1170, 1407, 3204, 2362, 870, 760, 2020, 60,
                             100, 4948, 3501, 4257, 597, 1537, 1616, 189, 1576, 976, 3190,
                             3020, 1584, 4952, 13, 4051, 3797, 3055, 3261, 4908, 3933, 2240,
                             4402, 1317, 2972, 2063, 1956, 4421, 3909, 4244, 2489, 2222,
                             4602, 3862, 371, 3867, 3056, 2210, 965, 534, 4698, 4369, 3552,
                             999, 2193, 4689, 3504, 4026, 3850, 4995, 886, 3914, 2681, 4316,
                             2633, 1229, 2160, 4186, 3505, 2570, 2068, 2315, 4858, 3564,
                             3173, 2739, 4467, 902, 3395, 314]) == 11231607
end


text:
module OSMPBF_pb
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  include("fileformat_pb.jl")
  include("osmformat_pb.jl")
end


text:
function returnFileLength(file_name::String)
    data1 = read(file_name, String)
    splitdata = split(data1, "\n")
    return length(splitdata) - 1
end

function findTotalPoints(folder_name::String)
    point_sum = 0
    for iter in 1:length(workers())
        if iter - 1 < 10
            filename = folder_name * "/" * "partGrid000" * string(iter-1)
        elseif iter - 1 < 100
            filename = folder_name * "/" * "partGrid00" * string(iter-1)
        elseif iter - 1 < 1000
            filename = folder_name * "/" * "partGrid0" * string(iter-1)
        else
            filename = folder_name * "/" * "partGrid" * string(iter-1)
        end
        data = read(filename, String)
        splitdata = @view split(data, "\n")[1:end-1]
        itmdata = split(splitdata[1], " ")
        local_point_count = parse(Int,itmdata[3])
        point_sum += local_point_count
    end
    return point_sum
end

function createGlobalLocalMapIndex(global_local_map_index, global_local_direct_index, folder_name::String, files_length, actual_files_length)
    index_flag = 1
    for iter in 1:length(workers())
        if iter - 1 < 10
            filename = folder_name * "/" * "partGrid000" * string(iter-1)
        elseif iter - 1 < 100
            filename = folder_name * "/" * "partGrid00" * string(iter-1)
        elseif iter - 1 < 1000
            filename = folder_name * "/" * "partGrid0" * string(iter-1)
        else
            filename = folder_name * "/" * "partGrid" * string(iter-1)
        end
        data = read(filename, String)
        splitdata = @view split(data, "\n")[1:end-1]
        itmdata = split(splitdata[1], " ")
        local_point_count = parse(Int,itmdata[3])
        ghost_point_count = parse(Int,itmdata[4])
        files_length[iter] = local_point_count
        actual_files_length[iter] = local_point_count + ghost_point_count
        for (idx, itm) in enumerate(splitdata)
            if idx == 1
                continue
            elseif idx <= local_point_count + 1
                itmdata = split(itm, " ")
                global_local_direct_index[parse(Int,itmdata[1])] = index_flag
                global_local_map_index[(parse(Float64,itmdata[2]), parse(Float64, itmdata[3]))] = index_flag
                index_flag += 1
            end
        end
    end
    return nothing
end

function readDistribuedFile(folder_name::String, defprimal, p, global_local_map_index)
    println("Reading multiple files")
    println(folder_name)
    iter = p - 1
    if iter - 1 < 10
        filename = folder_name * "/" * "partGrid000" * string(iter-1)
    elseif iter - 1 < 100
        filename = folder_name * "/" * "partGrid00" * string(iter-1)
    elseif iter - 1 < 1000
        filename = folder_name * "/" * "partGrid0" * string(iter-1)
    else
        filename = folder_name * "/" * "partGrid" * string(iter-1)
    end
    println(filename)
    data = read(filename, String)
    splitdata = @view split(data, "\n")[1:end-1]
    itmdata = split(splitdata[1], " ")
    local_point_count = parse(Int,itmdata[3])
    ghost_point_count = parse(Int,itmdata[4])
    local_points_holder = Array{Point,1}(undef, local_point_count)
    for (idx, itm) in enumerate(splitdata)
        if idx == 1
            continue
        elseif idx <= local_point_count + 1
            itmdata = split(itm, " ")
            globalID = global_local_map_index[(parse(Float64,itmdata[2]), parse(Float64, itmdata[3]))]
            local_points_holder[idx-1] = Point(parse(Int,itmdata[1]),
                parse(Float64,itmdata[2]),
                parse(Float64, itmdata[3]),
                parse(Int, itmdata[4]),
                parse(Int, itmdata[5]),
                parse(Int8,itmdata[6]),
                parse(Int8,itmdata[7]),
                parse(Float64,itmdata[8]),
                parse(Int8,itmdata[9]),
                parse.(Int, itmdata[10:end-1]),
                0.0,
                0.0,
                copy(defprimal),
                zeros(Float64, 4),
                Array{Array{Float64,1},1}(undef, 2), 0.0, 0, 0, 0, 0, Array{Int32,1}(undef, 0), Array{Int32,1}(undef, 0),
                Array{Int32,1}(undef, 0), Array{Int32,1}(undef, 0), 0.0, zeros(Float64, 4), zeros(Float64, 4), zeros(Float64, 4),
                globalID)
        end
    end
    return local_points_holder
end

function readDistribuedFileQ(folder_name::String, defprimal, p, global_local_map_index, files_length)
    println("Setting Up Q")
    # println(folder_name)
    iter = p - 1
    local_points_count = files_length[iter]

    local_points_holder = Array{TempQ,1}(undef, local_points_count)
    for idx in 1:local_points_count
        local_points_holder[idx] = TempQ(zeros(Float64, 4))
    end
    return local_points_holder
end

function readDistribuedFileDQ(folder_name::String, defprimal, p, global_local_map_index, files_length)
    println("Setting Up DQ")
    # println(folder_name)
    iter = p - 1
    local_points_count = files_length[iter]

    local_points_holder = Array{TempDQ,1}(undef, local_points_count)

    for idx in 1:local_points_count
        local_points_holder[idx] = TempDQ(Array{Array{Float64,1},1}(undef, 2))
    end

    return local_points_holder
end

function readGhostFile(folder_name::String, ghost_holder, global_local_map_index, dist_globaldata)
    # println(ghost_folder_name)
    for iter in 1:length(workers())
        if iter - 1 < 10
            filename = folder_name * "/" * "partGrid000" * string(iter-1)
        elseif iter - 1 < 100
            filename = folder_name * "/" * "partGrid00" * string(iter-1)
        elseif iter - 1 < 1000
            filename = folder_name * "/" * "partGrid0" * string(iter-1)
        else
            filename = folder_name * "/" * "partGrid" * string(iter-1)
        end
        println(filename)
        data = read(filename, String)
        splitdata = @view split(data, "\n")[1:end-1]
        ghost_holder[iter] = Dict{Int32,Point}()
        itmdata = split(splitdata[1], " ")
        local_point_count = parse(Int,itmdata[3])
        ghost_point_count = parse(Int,itmdata[4])

        for (idx, itm) in enumerate(splitdata)
            if idx > local_point_count + 1
                itmdata = split(itm, " ")
                ghost_holder[iter][idx-1] = dist_globaldata[global_local_map_index[(parse(Float64,itmdata[2]), parse(Float64, itmdata[3]))]]
            end
        end
    end
end

function readDistribuedFileMutables(folder_name::String, p, actual_files_length)
    println("Setting Up Mutables")
    # println(folder_name)
    iter = p - 1
    local_points_count = actual_files_length[iter]

    local_points_holder = zeros(Float64, 28, local_points_count)
    return local_points_holder
end

function writeToFile(loc_globaldata, numPoints)
    file = open("results/split/primvalscuda"* "_" * string(getConfig()["core"]["max_iters"]) *
     "_" * string(numPoints) * "_" * string(myid()) * ".txt", "w")
    @showprogress 1 "This takes time" for (idx, _) in enumerate(loc_globaldata)
        primtowrite = loc_globaldata[idx].prim
        for element in primtowrite
            @printf(file,"%0.17f", element)
            @printf(file, " ")
        end
        print(file, "\n")
    end
    close(file)
end

# function readGhostFileMutables(folder_name::String, ghost_holder_mutable)
#     # println(ghost_folder_name)
#     for iter in 1:length(workers())
#         if iter - 1 < 10
#             filename = folder_name * "/" * "partGrid000" * string(iter-1)
#         elseif iter - 1 < 100
#             filename = folder_name * "/" * "partGrid00" * string(iter-1)
#         elseif iter - 1 < 1000
#             filename = folder_name * "/" * "partGrid0" * string(iter-1)
#         else
#             filename = folder_name * "/" * "partGrid" * string(iter-1)
#         end
#         println(filename)
#         data = read(filename, String)
#         splitdata = @view split(data, "\n")[1:end-1]
#         ghost_holder_mutable[iter] = Dict{Int32,Array{Float64,1}}()
#         itmdata = split(splitdata[1], " ")
#         local_point_count = parse(Int,itmdata[3])
#         ghost_point_count = parse(Int,itmdata[4])
#         for (idx, itm) in enumerate(splitdata)
#             if idx > local_point_count + 1
#                 ghost_holder_mutable[iter][idx-1] = zeros(Float64, 20)
#             end
#         end
#     end
# end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
using Flux: Chain, Dense, σ, softmax, params, Conv, BatchNorm, relu
using Test
using Yama: flatten, unflatten_set

@testset "Yama.jl" begin
    @testset "Trajectory" begin
        include("test_trajectory.jl")
    end
end


text:
# if VERSION < v"0.5.0-dev+3294"
#     include("map0_4.jl")
# else
using Base: ith_all
if VERSION < v"0.5.0-dev+3294"
    include("map0_4.jl")
else
    using Base: collect_similar, Generator
end

function _return_type(f, Xs...)
    rtypes = Base.return_types(f, tuple([ inner_eltype(X) for X in Xs ]...))
    T = isempty(rtypes) ? Union{} : rtypes[1]
end

macro nullcheck(Xs, nargs)
    res = :($(Xs)[1].isnull[i])
    for i = 2:nargs
        e = :($(Xs)[$i].isnull[i])
        res = Expr(:||, res, e)
    end
    return res
end

macro fcall(Xs, nargs)
    res = Expr(:call, :f)
    for i in 1:nargs
        push!(res.args, :($(Xs)[$i].values[i]))
    end
    return res
end

inner_eltype{T}(X::NullableArray{T}) = T

# Base.map!

Base.map!{F}(f::F, X::NullableArray; lift=false) = map!(f, X, X; lift=lift)
function Base.map!{F}(f::F, dest::NullableArray, X::NullableArray; lift=false)
    if lift
        for (i, j) in zip(eachindex(dest), eachindex(X))
            if X.isnull[j]
                dest.isnull[i] = true
            else
                dest.isnull[i] = false
                dest.values[i] = f(X.values[j])
            end
        end
    else
        for (i, j) in zip(eachindex(dest), eachindex(X))
            dest[i] = f(X[j])
        end
    end
    return dest
end

function Base.map!{F}(f::F, dest::NullableArray, X1::NullableArray,
                      X2::NullableArray; lift=false)
    if lift
        for (i, j, k) in zip(eachindex(dest), eachindex(X1), eachindex(X2))
            if X1.isnull[j] | X2.isnull[k]
                dest.isnull[i] = true
            else
                dest.isnull[i] = false
                dest.values[i] = f(X1.values[j], X2.values[k])
            end
        end
    else
        for (i, j, k) in zip(eachindex(dest), eachindex(X1), eachindex(X2))
            dest[i] = f(X1[j], X2[k])
        end
    end
    return dest
end

function Base.map!{F}(f::F, dest::NullableArray, Xs::NullableArray...; lift=false)
    _mapn!(f, dest, Xs, lift)
end

@generated function _mapn!{F, N}(f::F, dest::NullableArray, Xs::NTuple{N, NullableArray}, lift)
    return quote
        if lift
            for i in eachindex(dest)
                if @nullcheck Xs $N
                    dest.isnull[i] = true
                else
                    dest.isnull[i] = false
                    dest.values[i] = @fcall Xs $N
                end
            end
        else
            for i in eachindex(dest)
                dest[i] = f(ith_all(i, Xs)...)
            end
        end
        return dest
    end
end

# Base.map

if VERSION < v"0.5.0-dev+3294"
    function Base.map(f, X::NullableArray; lift=false)
        lift ? _liftedmap(f, X) : _map(f, X)
    end
    function Base.map(f, X1::NullableArray, X2::NullableArray; lift=false)
        lift ? _liftedmap(f, X1, X2) : _map(f, X1, X2)
    end
    function Base.map(f, Xs::NullableArray...; lift=false)
        lift ? _liftedmap(f, Xs) : _map(f, Xs...)
    end
else
    function Base.map(f, X::NullableArray; lift=false)
        lift ? _liftedmap(f, X) : collect_similar(X, Generator(f, X))
    end
    function Base.map(f, X1::NullableArray, X2::NullableArray; lift=false)
        lift ? _liftedmap(f, X1, X2) : collect(Generator(f, X1, X2))
    end
    function Base.map(f, Xs::NullableArray...; lift=false)
        lift ? _liftedmap(f, Xs) : collect(Generator(f, Xs...))
    end
end

function _liftedmap(f, X::NullableArray)
    len = length(X)
    # if X is empty, fall back on type inference
    len > 0 || return NullableArray(_return_type(f, X), 0)
    i = 1
    while X.isnull[i] & (i < len)
        i += 1
    end
    # if X is all null, fall back on type inference
    if X.isnull[i]
        T = _return_type(f, X)
        return similar(X, T)
    end
    # otherwise, initialize and map to destination array
    v = f(X.values[i])
    dest = similar(X, typeof(v))
    dest[i] = v
    _liftedmap_to!(f, dest, X, i+1, len)
end

function _liftedmap(f, X1::NullableArray, X2::NullableArray)
    len = prod(promote_shape(size(X1), size(X2)))
    len > 0 || return NullableArray(_return_type(f, X1, X2), 0)
    i = 1
    while (X1.isnull[i] | X2.isnull[i]) & (i < len)
        i += 1
    end
    if X1.isnull[i] | X2.isnull[i]
        T = _return_type(f, X1, X2)
        return similar(X1, T)
    end
    v = f(X1.values[i], X2.values[i])
    dest = similar(X1, typeof(v))
    dest[i] = v
    _liftedmap_to!(f, dest, X1, X2, i+1, len)
end

@generated function _liftedmap{N}(f, Xs::NTuple{N, NullableArray})
    return quote
        shp = mapreduce(size, promote_shape, Xs)
        len = prod(shp)
        len > 0 || return NullableArray(_return_type(f, Xs...), 0)
        i = 1
        while (@nullcheck Xs $N) & (i < len)
            i += 1
        end
        if @nullcheck Xs $N
            T = _return_type(f, Xs...)
            return similar(Xs[1], T)
        end
        v = @fcall Xs $N
        dest = similar(Xs[1], typeof(v))
        dest[i] = v
        _liftedmap_to!(f, dest, Xs, i+1, len)
    end
end

function _liftedmap_to!{T}(f, dest::NullableArray{T}, X, offs, len)
    # map to dest array, checking the type of each result. if a result does not
    # match, widen the result type and re-dispatch.
    i = offs
    while i <= len
        @inbounds if X.isnull[i]
            i += 1; continue
        end
        @inbounds el = f(X.values[i])
        S = typeof(el)
        if S === T || S <: T
            @inbounds dest[i] = el::T
            i += 1
        else
            R = typejoin(T, S)
            new = similar(dest, R)
            copy!(new, 1, dest, 1, i-1)
            @inbounds new[i] = el
            return map_to!(f, new, X, i+1, len)
        end
    end
    return dest
end

function _liftedmap_to!{T}(f, dest::NullableArray{T}, X1, X2, offs, len)
    i = offs
    while i <= len
        @inbounds if X1.isnull[i] | X2.isnull[i]
            i += 1; continue
        end
        @inbounds el = f(X1.values[i], X2.values[i])
        S = typeof(el)
        if S === T || S <: T
            @inbounds dest[i] = el::T
            i += 1
        else
            R = typejoin(T, S)
            new = similar(dest, R)
            copy!(new, 1, dest, 1, i-1)
            @inbounds new[i] = el
            return map_to!(f, new, X1, X2, i+1, len)
        end
    end
    return dest
end

@generated function _liftedmap_to!{T, N}(f, dest::NullableArray{T}, Xs::NTuple{N,NullableArray}, offs, len)
    return quote
        i = offs
        while i <= len
            @inbounds if @nullcheck Xs $N
                i += 1; continue
            end
            @inbounds el = @fcall Xs $N
            S = typeof(el)
            if S === T || S <: T
                @inbounds dest[i] = el::T
                i += 1
            else
                R = typejoin(T, S)
                new = similar(dest, R)
                copy!(new, 1, dest, 1, i-1)
                @inbounds new[i] = el
                return map_to!(f, new, Xs, i+1, len)
            end
        end
        return dest
    end
end


text:
"""
credits to: @stevengj professor steven G. Johnson of
the MIT 18.335J / 6.337J -  Numerical Analysis 2019 by
18.335J / 6.337J
"""
#stableRNG
rng = StableRNG(123)
A = randn(rng, 10, 10)
num1 = randn(rng, 1)
num2 = randn(rng, 1) # instead of randn(10, 10)
@test inv(inv(A)) ≈ A
@test inv(randn(A))

struct ratio()
    numerator
    denominator
    return numerator / denominator
end
global maxiterations = 10
@testset genericTest begin

    function half(n)
        enumerate(n)
           x = 1/a /2
        end
    end

    function iterate_sqrt(a, n, x)
        """
        a = numerator
        x = denominator
        """
        for i = 1:n
            x = (x + a / x) / 2  #parabolic shape
        end
        return x
    end

    for (i, value) in enumerate(a)
        println("$i $value")
    end

    function fun1(ϵ, maxiterations, x = 2, min = 1, fun = iterate_sqrt(x, iterations, min))
# i doesn't matter
        for iterations = min:maxiterations
            val = iterate_sqrt(x, iterations, min)
            ϵ = val - sqrt(x) # 1 form of difference (euclidean distance )[the Geodesic ] (optimally. should be absolute)
            @show iterations, val, ϵ

            return iterations, val, ϵ
        end
    end

function fun1(ϵ, maxiterations, x = 2, min = 1, fun = iterate_sqrt(x, iterations, min))
    # i doesn't matter
    for i in enumerate(maxiterations)
        val = iterate_sqrt(x, i, min)
        ϵ = val - sqrt(x) # 1 form of difference (euclidean distance )[the Geodesic ] (optimally. should be absolute)
        @show i, val, ϵ

        return i, val, ϵ
    end
#end

    """
    in julia, defaults to a Double-precision, floating-point values, or like 15 decimal digits)
    """
function error()

end
    function printIterations(BigFloat, maxiterations = 16) # 400 binary digits, about 120 decimal places
        #enumerate(maxiterations)
        for iterations = 1:10
            val = iterate_sqrt(big"2.0", iterations, big"1.0")
            err = Float64(val - sqrt(big"2.0"))
            println("iterations = $iterations:\n    val = $val\n    err = $err")
        end
    end

     err = Float64(iterate_sqrt(big"2.0", n, big"1.0") - sqrt(big"2.0"))
    function setprecision()
        #enumerate(maxiterations)
        n= 10

            enumerate(n)
            val = iterate_sqrt(big"2.0", n, big"1.0")
            err = Float64(val - sqrt(big"2.0"))
            println("iterations = $iterations:\n    val = $val\n    err = $err")
        end
    end

    # output x as HTML, with digits matching x0 printed in bold
    function matchdigits(message, pattern)
        s = string(message)
        pattern = string(pattern)
        buf = IOBuffer()
        matches = true
        i = 0
        print(buf, "")
        while (i += 1) <= length(s)
            i % 70 == 0 && print(buf, "
")
            if matches && i <= length(pattern) && isdigit(message[i])
                if message[i] == pattern[i]
                    print(buf, message[i])
                    continue
                end
                print(buf, "")
                matches = false
            end
            print(buf, message[i])
        end
        matches && print(buf, "")
        HTML(String(take!(buf)))
    end
    using Interact
    range = 0:12
    function setprecision(BigFloat, fun = iterate_sqrt(x, range, iterations)) # 4000 binary digits, about 1200 decimal places

        @manipulate for n = 0:12
            print(typeof(n))
            matchdigits(iterate_sqrt(big(2), n, big(1)), sqrt(big(2)))
        end
    end

    #

    setprecision(BigFloat, maxiterations) # =16
    for range = 0:10
        display(matchdigits(iterate_sqrt(big(2), n, big(1)), sqrt(big(2))))
    end
end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
module HistogramThresholding

using Images
using LinearAlgebra

abstract type ThresholdAlgorithm end
struct Otsu <: ThresholdAlgorithm end
struct MinimumError <: ThresholdAlgorithm end
struct MinimumIntermodes <: ThresholdAlgorithm end
struct Intermodes <: ThresholdAlgorithm end
struct UnimodalRosin <: ThresholdAlgorithm end
struct Moments <: ThresholdAlgorithm end

include("common.jl")
include("otsu.jl")
include("minimum_error.jl")
include("unimodal.jl")
include("moments.jl")
include("minimum.jl")
include("intermodes.jl")

export
  # main functions
  find_threshold,
  Otsu,
  UnimodalRosin,
  Moments,
  MinimumIntermodes,
  Intermodes,
  MinimumError
end # module


text:
module Saving

using Dates: TimeType

import DrWatson
import ..ControlProblem

struct OptimizationConfig
    problem::ControlProblem
    method::Symbol
    kwargs::Dict{Symbol,Any}
end

function allaccess(c)
    no_anon_funcs((k, v)) = !((v isa Function) && contains(String(nameof(v)), "#"))
    all_keys = String["method"]
    for dict in (c.kwargs, c.problem.kwargs)
        for (k, v) ∈ dict
            key = String(k)
            if !((v isa Function) && contains(String(nameof(v)), "#"))
                push!(all_keys, key)
            end
        end
    end
    return all_keys
end

DrWatson.allaccess(c::OptimizationConfig) = allaccess(c)

function access(c, key::AbstractString)
    (key == "method") && (return String(c.method))
    key_sym = Symbol(key)
    if key_sym ∈ keys(c.kwargs)
        val = c.kwargs[key_sym]
    else
        val = c.problem.kwargs[key_sym]
    end
    return val
end

DrWatson.access(c::OptimizationConfig, key::AbstractString) = access(c, key)

const ALLOWED_TYPES = [Real, String, Symbol, TimeType, Function]

function default_allowed(c)
    return ALLOWED_TYPES
end

DrWatson.default_allowed(c::OptimizationConfig) = default_allowed(c)

end


text:
using BEAST
using CompScienceMeshes
#using LinearForms
using MATLAB
#include(Pkg.dir("BEAST","src","lusolver.jl"))

function patch_mat(Γ,fcr=nothing)
  fcr != nothing ?  C = norm.(fcr) : C=0
  V = vertexarray(Γ)
  F = cellarray(Γ)
  mat"clf"
  @mput V F C
  @matlab begin
      hold("on")
      patch("Vertices", V, "Faces", F, "FaceVertexCData", C, "FaceColor", "flat")
  end

  mat"axis equal"
  #mat"colormap jet"
  mat"xlabel('x-axis')"
  mat"ylabel('y-axis')"
  mat"colorbar"
  mat"title('Plot of Surface currents J')"
  mat"view(3)"
end

# define geometry & material parameters
f₀= 3.0e6 # 3 GHz
 c = 3.0e8
 ϵ₀ = 8.85418782e-12
 μ₀ = 1.25663706e-6
 ϵᵣ=1; μᵣ=1
 l = w = 1.0 #Length and width of capacitor plates
 d = 0.1    #seperation of plates
 h = 1/12   #size of meshes
 ω = 2π*f₀
 vₑ = c/sqrt(ϵᵣ)
 λₑ = vₑ/f₀
 η = sqrt(μ₀*μᵣ/ϵ₀*ϵᵣ)
 κ = ω/vₑ

println("wavenumber = ", κ)
println("freq = ", f₀, " Hz")
println("λ = ", λₑ, "m; l = ", l ,"m")

# Build and mesh the structure
Γ₁ = meshrectangle(l,w,h);
 translate!(Γ₁, point(0.0,0.0,d))
Γ₀ = meshrectangle(l,w,h)

γ₁ = meshsegment(l, h, 3)
  translate!(γ₁, point(0.0,0.0,d))
γ₀ = meshsegment(l, h, 3)

Γ = weld(Γ₁, Γ₀)

# define the excitation
V₀ = 1.0
f = ScalarTrace(p -> V₀ )

#define basis function
RT = rt_ports(Γ,[γ₁,γ₀])

# define the equations
j, = hilbertspace( :j)
k, = hilbertspace( :m)
T = MWSingleLayer3D(im*κ)
EFIE = @varform η*T[k,j] == f[ntrace(k,γ₁)]

# discretise & solve the equation
efie = @discretise EFIE j∈RT k∈RT
u = solve(efie)

#Calculate & plot face currents
fcr,geo = facecurrents(u, RT); println("Face currents calculated")
patch_mat(geo, fcr); println("Facecurrents plotted")

# Compute the Scalar Potential across the ports
pts1 = [point(0.5,0.5,a) for a in linspace(-0.5,0.5,100) ]
volt = η*potential(MWSingleLayerPotential3D(κ), pts1, u, RT)

# Plot the Scalar potential
using Plots
  plotlyjs()
Φ = [real(v[1]) for v in volt]
Plots.plot(Φ)
println("Potential plotted")

# Calculate total charge Q
idx = getindex_rtg(RT) #get index of global RWG
Q = norm(-u[idx]/(im*ω))

#Check capacitance with
C = (ϵ₀*ϵᵣ*(l*w))/d


text:
using InfrastructureSystemsPlots
using DataFrames
using Dates
const ISP = InfrastructureSystemsPlots

variables = Dict()
variables[:P_ThermalStandard] = DataFrames.DataFrame(
    :one => [1, 2, 3, 2, 1],
    :two => [3, 2, 1, 2, 3],
    :three => [1, 2, 3, 2, 1],
)
variables[:P_RenewableDispatch] = DataFrames.DataFrame(
    :one => [3, 2, 3, 2, 3],
    :two => [1, 2, 1, 2, 1],
    :three => [3, 2, 3, 2, 3],
)
optimizer_log = Dict()
objective_value = Dict()
right_now = round(Dates.now(), Dates.Hour)
time_stamp =
    DataFrames.DataFrame(:Range => right_now:Dates.Hour(1):(right_now + Dates.Hour(4)))
res = ISP.Results(variables, optimizer_log, objective_value, time_stamp)

generators = Dict("Coal" => [:one; :two], "Wind" => [:three])