text:
#=
This file attempts to implement the "imporoved dimension estimation" by Judd:

Judd, K. (1992). An improved estimator of dimension and some comments on providing
confidence intervals. Physica D: Nonlinear Phenomena, 56(2–3), 216–228. https://doi.org/10.1016/0167-2789(92)90025-I

Judd, K. (1994). Estimating dimension from small samples. Physica D: Nonlinear Phenomena,
71(4), 421–429. https://doi.org/10.1016/0167-2789(94)90008-6

So far this attempt has failed.
=#
using Distances, BigCombinatorics

"""
    interpoint_distances(X, N, norm = Euclidean()) -> ds

returns N interpoint distances of randomly chosen pairs of points from time series X.
N must be smaller than the free memory. N ≈ 5000 is actually already sufficient for the histograms. N ≈ 45000 sort of the upper limit of what you should use, this makes for a smoother hist later
"""
function interpoint_distances(X, N, norm = Euclidean())
    ds = eltype(X)[]

    @assert N < Int(Sys.free_memory()) ÷ 8

    # choosing N random points from X.
    sizehint!(ds, Int(N*(N-1)/2))
    while length(ds) < Int(N*(N-1)/2)
        i,j = rand(1:length(X),2)
        dist = evaluate(norm, X[i], X[j])
        (i != j && !any(ds .== dist)) && push!(ds, dist)
    end
    return filter(x-> x!= 0.0,ds)
end

"""
    logspace_histogram(distances::Vector, λ = exp(-w0), ϵ₀ = λ^2 * maximum(x)) -> bins, counts
Calculate a logarithmically spaced histogram of the values in `x` (must be positive
definite, i.e. distances) and return the bins and count of values in the bins.

The `i`-th bin is `[ϵ₀*λ^i, ϵ₀*λ^(i+1))`. Elements of `x` greater than `ϵ₀` are disregarded.
"""
function logspace_histogram(distances::AbstractVector, ε₀, λ = exp(-w0(x)))
    @assert λ < 1
    εs = [ε₀]
    λmin = minimum(distances)*λ
    εᵢ = ε₀
    while εᵢ > λmin
        εᵢ = εs[end]*λ
        push!(εs, εᵢ)
    end

    reverse!(εs)
    L = length(εs)
    bs = zeros(Int, L)
    # Fill bins with values
    for v in distances
        i = searchsortedfirst(εs, v)
        i == L + 1 && continue
        bs[i] += 1
    end
    return reverse!(εs), reverse!(bs)
end
# TODO: test the algorithm with including the bin [∞, ε0)
w0(x) = log(maximum(x)/minimum(x))/√length(x)

"""
    bin_judd(x::AbstractVector) -> εs, bs, ε₀

returns a logspace histogram with optimal bin width.
Estimation of the latter could probably be done more efficient

"""
function bin_judd(x::AbstractVector)

    opt_λ =  optimal_λ(x) # This is probably not very efficient
    εs, bs = logspace_histogram(x,maximum(x),opt_λ)

    ε₀ = εs[findmax(bs)[2]+1]
    new_εs = reverse(εs[findmax(bs)[2]+1:end])
    new_bs = reverse(bs[findmax(bs)[2]+1:end])
    push!(new_bs, sum(bs[1:findmax(bs)[2]]))
    push!(new_εs, εs[findmax(bs)[2]])
    reverse!(new_bs)
    reverse!(new_εs)
    return new_εs, new_bs, ε₀
end

function P(εᵢ, d,a,ε₀)
    polynomial = length(a) == 0 ? 1 : 0
    for t in 1:length(a)
        polynomial += a[t] * (εᵢ/ε₀)^(t-1)
    end
    return (εᵢ/ε₀)^d * polynomial
end

""" log_l(εs, bs, P, d, a,ε₀)

returns the log-likelihood function according to Judd
"""
log_l(εs, bs, P, d, a,ε₀) = -sum(bs[i]*log( P(εs[i],d,a,ε₀) - P(εs[i+1],d,a,ε₀) ) for i in 1:length(εs)-1 )

""" optimal_λ(dists)

searches for the optimal bin width w = log(1/λ) of the judd histogram.
The search is currently not how Judd proposed, because I didn't get anywhere with that. Right now, this is a brute force search, this should definitely be improved if this is actuallly going to be used.

"""
function optimal_λ(dists)
    n = length(dists)
    w₀ = w0(dists)
    ws = collect(w₀:0.001:.5)
    f_ws = []

    for w in ws
        εs, bs = logspace_histogram(dists, maximum(dists), exp(-w))
        m = findlast(x->x!=0,bs)
        f = log(w) + log(Multinomial(bs[1:m])) + log(Binomial(n+m+1,m))
        push!(f_ws,f)
    end

    w = ws[findmin(f_ws)[2]]
    λ = exp(-w)
    return λ
end

problargerzero(εs, bs, d₀,a₀,ε₀) = all(P(εs[i], d₀, a₀,ε₀)-P(εs[i+1],d₀,a₀,ε₀) > 0 for i in 1:length(εs)-1)

""" nested_optimizer(εs, bs; ε₀, d₀, a₀, stp = 1e-5, η = 1e-3, tol=1e-4, decay=1e-4, maxiter = 10000)

    uses numerical (center-diff) gradient descent with adaptive (decaying) learning rate to find
    a minimum of the log-liklelihood function.
    Seems to produce reasonable results if `deg(a) < 2`, but still sensitive to initial
    conditions.
    The key was to use the nested version that Judd suggested in a footnote of the 1994 paper:
    nest the optimizations so that first `d` is optimized with fixed `a`, then the coefficients
    of `a` are optimized with fixed `d`. Repeat until total step size is smaller than `tol`.

    There are a lot of parameters to this optimization problem that themselves can probably be
    optimized and/ or made more adaptive so that they fit the specific problem.

"""
function nested_optimizer(εs, bs; ε₀, d₀, a₀, stp = 1e-5, η = 1e-3, γ = 1e-3, tol=1e-4, decay=1e-4, maxiter = 10000)

    iter = 0

    deg_a = length(a₀) # Technically, this is not the *degree*, but deg(a)+1, but it's more convenient to write it like this

    cur_x = [d₀, a₀...]
    cur_d = cur_x[1]
    cur_a = cur_x[2:end]

    step_size = 10

    while iter < maxiter && step_size > tol

        prev_v = cur_x

        # First we will optimize d with fixed polynomial coefficients
        cur_d = cur_x[1]
        v = 0
        iter_d = 0
        d_step = 1
        while iter_d < maxiter÷10 && d_step > tol/10
            prev_d = cur_d
            g = (log_l(εs, bs, P, cur_d+stp, a₀,ε₀)-log_l(εs, bs, P, cur_d-stp, a₀,ε₀))/(2stp)
            v = η * g + γ * v
            cur_d -= v
            d_step = abs(cur_d-prev_d)
            iter_d += 1
        end

        # now we will optimize the coefficients of a with the previous best d
        grad = zeros(deg_a)
        cur_a = cur_x[2:end]
        v = zeros(deg_a)
        iter_a = 0
        a_step = 1
        while iter_a < maxiter÷10 && a_step > tol/10
            prev_a = cur_a
            for k in 1:deg_a
                a_stp = zeros(deg_a)
                a_stp[k] = stp
                grad[k] = (log_l(εs, bs, P, d₀, a₀ + a_stp,ε₀)-log_l(εs, bs, P, d₀, a₀-a_stp,ε₀))/(2stp)
            end
            v = γ * v + η * grad
            cur_a -= v
            a_step = evaluate(Euclidean(), cur_a, prev_a)
            iter_a += 1
        end

        if problargerzero(εs, bs, cur_d,cur_a,ε₀)
            cur_x = [cur_d, cur_a...]
        else
            break
        end

        # check total step size of parameter vector
        step_size = evaluate(Euclidean(), cur_x, prev_v)

        # update learning rate
        η /= (1+decay*iter)

        iter += 1
    end
    cur_d, cur_a
end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
#Generate Parameters of WNNM-algorithm
#input the nSig  the noise Level

mutable struct PAR{T<:Int64, F<:Float32}
    nSig        ::T                         #level of Gaussian noise
    SearchWin   ::T                         #the search range for Neighbour Patches
    delta       ::F                         #parameter of iterative method
    c           ::F
    Innerloop   ::T

    patsize     ::T                         #the size of Patches is: patsize * patsize
    patnum      ::T
    Iter        ::T                         #parameter of iterative method
    lamada      ::F
    step        ::T
end

function autopar(nSig)

    SearchWin =   30
    delta     =   Float32(0.1)
    c         =   Float32(3*√(2))
    Innerloop =   2

    if          nSig <= 20
        patsize       =   6
        patnum        =   70
        Iter          =   8
        lamada        =   Float32(0.54)
    elseif      nSig <= 40
        patsize       =   7
        patnum        =   90
        Iter          =   12
        lamada        =   Float32(0.56)
    elseif      nSig <= 60
        patsize       =   8
        patnum        =   120
        Iter          =   14
        lamada        =   Float32(0.58)
    else
        patsize       =   9
        patnum        =   140
        Iter          =   14
        lamada        =   Float32(0.58)
    end
    step              =   patsize ÷ 2 -1

    PAR(nSig, SearchWin, delta, c, Innerloop, patsize, patnum , Iter, lamada, step)
end

ext: jl


text:
module TEST13H_examples
using FinEtools
using FinEtools.AlgoDeforLinearModule: ssit
using Gaston

function TEST13H_hva()
    # Harmonic forced vibration problem is solved for a homogeneous square plate,
    # simply-supported on the circumference.
    # This is the TEST 13H from the Abaqus v 6.12 Benchmarks manual.
    # The test is recommended by the National Agency for Finite Element Methods and Standards (U.K.):
    # Test 13 from NAFEMS “Selected Benchmarks for Forced Vibration,” R0016, March 1993.
    #
    #
    # The plate is discretized with hexahedral solid elements. The simple support
    # condition is approximated by distributed rollers on the boundary.
    # Because only the out of plane displacements are prevented, the structure
    # has three rigid body modes in the plane of the plate.
    #
    #
    # The nonzero benchmark frequencies are (in hertz): 2.377, 5.961, 5.961,
    # 9.483, 12.133, 12.133, 15.468, 15.468 [Hz].

    println("""
    Homogeneous square plate, simply-supported on the circumference from
    the test 13 from NAFEMS “Selected Benchmarks for Forced Vibration,” R0016, March 1993.
    The nonzero benchmark frequencies are (in hertz): 2.377, 5.961, 5.961,
    9.483, 12.133, 12.133, 15.468, 15.468 [Hz].
    """)

    # t0 = time()

    E = 200*phun("GPa");# Young's modulus
    nu = 0.3;# Poisson ratio
    rho = 8000*phun("KG*M^-3");# mass density
    qmagn = 100.0*phun("Pa")
    L = 10.0*phun("M"); # side of the square plate
    t = 0.05*phun("M"); # thickness of the square plate
    nL = 16; nt = 4;
    tolerance = t/nt/100;
    # neigvs = 11;
    # OmegaShift = (2*pi*0.5) ^ 2; # to resolve rigid body modes
    frequencies = vcat(linearspace(0.0,2.377,20), linearspace(2.377,15.0,70))

    # Compute the parameters of Rayleigh damping. For the two selected
    # frequencies we have the relationship between the damping ratio and
    # the Rayleigh parameters
    # $\xi_m=a_0/\omega_m+a_1\omega_m$
    # where $m=1,2$.  Solving for the Rayleigh parameters $a_0,a_1$ yields:
    zeta1= 0.02; zeta2  =0.02;
    o1 =2*pi*2.377;  o2 =2*pi*15.468;
    Rayleigh_mass = 2*(o1*o2)/(o2^2-o1^2)*(o2*zeta1-o1*zeta2);# a0
    Rayleigh_stiffness = 2*(o1*o2)/(o2^2-o1^2)*(-1/o2*zeta1+1/o1*zeta2);# a1

    Rayleigh_mass = Rayleigh_mass;
    Rayleigh_stiffness = Rayleigh_stiffness;

    MR = DeforModelRed3D
    fens,fes  = H8block(L, L, t, nL, nL, nt)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(FCplxFlt, size(fens.xyz,1), 3)) # displacement field
    nl = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, 3)
    nl = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, 3)
    nl = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, 3)
    nl = selectnode(fens, box=[-Inf Inf L L -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, 3)
    applyebc!(u)
    numberdofs!(u)
    println("nfreedofs = $(u.nfreedofs)")

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearMSH8(MR, IntegData(fes, GaussRule(3,2)), material)
    femm = associategeometry!(femm, geom)
    K = stiffness(femm, geom, u)
    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,3)), material)
    M = mass(femm, geom, u)
    C = Rayleigh_mass*M + Rayleigh_stiffness*K

    # if true
    #     t0 = time()
    #     d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    #     d = d - OmegaShift;
    #     fs = real(sqrt.(complex(d)))/(2*pi)
    #     println("Reference Eigenvalues: $fs [Hz]")
    #     println("eigs solution ($(time() - t0) sec)")
    # end

    bdryfes = meshboundary(fes)
    topbfl = selectelem(fens, bdryfes, facing=true, direction=[0.0 0.0 1.0])
    el1femm =  FEMMBase(IntegData(subset(bdryfes,topbfl), GaussRule(2,2)))
    function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
        forceout .=  [0.0, 0.0, -qmagn]
        return forceout
    end
    fi = ForceIntensity(FFlt, 3, pfun);
    F = distribloads(el1femm, geom, u, fi, 2);

    U1 = zeros(FCplxFlt, u.nfreedofs, length(frequencies))
    for k = 1:length(frequencies)
        frequency = frequencies[k];
        omega = 2*pi*frequency;
        U1[:, k] = (-omega^2*M + 1im*omega*C + K)\F;
    end

    midpoint = selectnode(fens, box=[L/2 L/2 L/2 L/2 0 0], inflate=tolerance);
    midpointdof = u.dofnums[midpoint, 3]

    umidAmpl = abs.(U1[midpointdof, :])/phun("MM")
    set(axis="semilogx", plotstyle="linespoints", linewidth=4, pointsize = 2, color = "blue", xlabel = "Frequency [Hz]", ylabel = "Midpoint  displacement amplitude [mm]", grid="on", title = "Thin plate midpoint Amplitude FRF")
    f = figure()
    plot(vec(frequencies), vec(umidAmpl), legend = "", marker = "edmd")
    figure(f)

    umidReal = real.(U1[midpointdof, :])/phun("MM")
    umidImag = imag.(U1[midpointdof, :])/phun("MM")
    set(axis="semilogx", plotstyle="linespoints", linewidth=4, pointsize = 2, xlabel = "Frequency [Hz]", ylabel = "Displacement FRF [mm]", grid="on", title = "Thin plate midpoint Real/Imag FRF")
    f = figure()
    plot(vec(frequencies), vec(umidReal), color = "red", legend = "real", marker = "ecircle")
    plot!(vec(frequencies), vec(umidImag), color = "blue", legend = "imag", marker = "x")
    figure(f)

    umidPhase = atan2.(umidImag,umidReal)/pi*180
    set(axis="semilogx", plotstyle="linespoints", linewidth=4, pointsize = 2, xlabel = "Frequency [Hz]", ylabel = "Phase shift [deg]", grid="on", title = "Thin plate midpoint FRF Phase")
    f = figure()
    plot(vec(frequencies), vec(umidPhase), color = "blue", legend = "", marker = "x")
    figure(f)
    true
end # TEST13H_hva

function allrun()
    println("#####################################################")
    println("# TEST13H_hva ")
    TEST13H_hva()
    return true
end # function allrun

end # module TEST13H_examples


text:
function time_sec()
    time_ns() * 1e-9
end

function timed_result(action::Function, seed::Int64)
    start = time_sec()
    result = action(seed)
    delta = time_sec() - start
    return (delta, result)
end

function time_clock(niter::Integer)
    start = time_sec()
    for _ in 1:niter
        time_sec()
    end
    time_sec() - start
end

function timed_noresult(action::Function, seed::Int64)
    start = time_sec()
    for _ in 1:seed
        action()
    end
    time_sec() - start
end

function run_for_atleast(howlong::Float64,
			 seed::Int64,
			 action::Function)
    #ensure that time_sec is compiled
    time_sec()
    init_time = time_sec()
    iters = 0
    while true
        now = time_sec()
        if ((now - init_time) > 10.0 * howlong)
           throw(error("took too long to run: seed $seed, iters $iters"))
        end
	elapsed, result = timed_result(action, seed)
	if elapsed < howlong
	    seed  *= 2
     	    iters += 1
   	else
 	    return elapsed, seed, result
	end
     end
end

function timed_func(func::Function, n::Integer)
    start = time_ns()
    local ret
    for _ in 1:n
	ret = func()
    end
    (int(time_ns() - start), 1.0)
end

function collect_samples(sample_count::Integer,
			 exec_count::Integer,
			 func::Function,
			 gc_before::Bool)
    @assert sample_count > 0
    results = zeros(sample_count)
    times = zeros(sample_count)
    for i in 1:sample_count
        if gc_before
	    Base.gc()
	end
        t, r = timed_func(func, exec_count)
        times[i] = t;
        results[i] = r
    end
    (times, results)
end

function estimate_execution_count(period, func,
			          gc_before_sample,
				  est_run_time)
    print("Estimating execution count...")
    n = int(max(period / est_run_time / 5, 1))
    while true
        times = collect_samples(1, n, func, gc_before_sample)[1]
        t = max(1.0, mean(times)) # prevent zero times
   	print("....")
        if t >= period
            print("\n")
            return n
        end
        n = min(2 * n, n * (period / t) + 1)
    end
    @printf("%d\n", n)
    return n
end

function time_str(k::Float64)

    function fmt(time::Float64, unit::String)
        if     time >= 1e9 return @sprintf("%.4f %s", time, unit)
        elseif time >= 1e6 return @sprintf("%.0f %s", time, unit)
        elseif time >= 1e5 return @sprintf("%.1f %s", time, unit)
        elseif time >= 1e4 return @sprintf("%.2f %s", time, unit)
        elseif time >= 1e3 return @sprintf("%.3f %s", time, unit)
        elseif time >= 1e2 return @sprintf("%.4f %s", time, unit)
        elseif time >= 1e1 return @sprintf("%.5f %s", time, unit)
        else return @sprintf("%.6f %s", time, unit) end
    end

    if k < 0          return "-" + secs(-k)
    elseif k >= 60    return fmt(k / 60.0, "min")
    elseif k >= 1     return fmt(k, "s")
    elseif k >= 1e-03 return fmt(k * 1e03, "ms")
    elseif k >= 1e-06 return fmt(k * 1e06, "us")
    elseif k >= 1e-09 return fmt(k * 1e09, "ns")
    elseif k >= 1e-12 return fmt(k * 1e12, "ps")
    else return @sprintf("%.6f %s", k * 1e12, "ps")
    end
end


text:
module AnalyticDiffEq

using ModelingToolkit

abstract type AbstractEquation end
abstract type LaguerreAbstractEquation <: AbstractEquation end
abstract type HermiteAbstractEquation <: AbstractEquation end
abstract type LegendreAbstractEquation <: AbstractEquation end

struct LaguerreEquation{α,n} <: LaguerreAbstractEquation end
struct HermiteEquation{n} <: HermiteAbstractEquation end
struct LegendreEquation{n} <: LegendreAbstractEquation end

# struct HermiteEquation{n} <: AbstractEquation end
# struct LegendreEquation{n} <: AbstractEquation end
# struct ChebyshevFirstKindEquation{n} <: AbstractEquation end
# struct ChebyshevSecondKindEquation{n} <: AbstractEquation end
# struct GegenbauerEquation{α,n} <: AbstractEquation end
# struct JacobiEquation{α, β, n} <: AbstractEquation end

# Laguerre
# xy'' + (α + 1 - x)* y' + ny = 0
# eqtest1 = x*D2(y) + (1 + 1 - x) * D(y) + 1*y ~ 0
function isit(eq :: Equation, y :: Type{T}) where T<:AbstractEquation
        check_lhs(eq, y) && check_rhs(eq, y)
end

check_rhs(eq :: Equation, y :: Type{T}) where T<:AbstractEquation = eq.rhs.value == 0
check_lhs(eq :: Equation, y :: Type{T}) where T<:AbstractEquation = check_term1(eq,y) && check_term2(eq,y) && check_term3(eq,y)

function check_term1(eq :: Equation, y :: Type{LaguerreAbstractEquation})
        term1 = eq.lhs.args[1].args[1]
        flags = Bool[]

        flag1 = isequal(term1.op, *)
        push!(flags, flag1)

        flag2 = length(term1.args) == 2
        push!(flags, flag2)

        flag3 = typeof(term1.args[1]) == Variable
        push!(flags, flag3)

        x = term1.args[1]
        flag4 = x.known == false
        push!(flags, flag4)

        flag5 = x.name == :x # cross check dependents
        push!(flags, flag5)

        flag6 = length(x.dependents) == 1
        push!(flags, flag6)

        DD2 = term1.args[2]
        flag7 = typeof(DD2.op) == Differential
        push!(flags, flag7)

        # Time as a dependent?
        y = DD2.args[1].args[1]
        flag8 = isequal(y.dependents, x.dependents)
        push!(flags, flag8)

        flag9 = typeof(DD2.args[1].op) == Differential
        push!(flags, flag9)

        flag10 = y.name == :y
        push!(flags, flag10)

        flag11 = length(y.dependents) >= 1
        push!(flags, flag11)

        return all(flags)
end
function check_term2(eq :: Equation, y :: Type{LaguerreAbstractEquation})
        term2 = eq.lhs.args[1].args[2]
        flags = Bool[]

        flag1 = length(term2.args) == 2
        push!(flags, flag1)

        # check that it is indeed the derivative we want
        flag2 = typeof(term2.args[2].op) == Differential
        push!(flags, flag2)

        y = term2.args[2].args[1] # check names and dependents
        flag3 = length(y.dependents) == 1

        flag4 = typeof(y) == Variable
        push!(flags,flag4)

        # check it is an algebraic form we want
        flag5 = length(term2.args[1].args) == 2
        push!(flags, flag5)

        x = term2.args[1].args[2] # check names and dependents
        flag6 = typeof(x) == Variable
        push!(flags, flag6)

        # alpha + 1 TODO assume alpha is more than 0
        a = term2.args[1].args[1].value
        α = a - 1
        flag7 = α isa Real
        push!(flags, flag7)

        flag8 = term2.args[1].op == -
        push!(flags, flag8)

        return all(flags)
end
function check_term3(eq :: Equation, y :: Type{LaguerreAbstractEquation})
        term3 = eq.lhs.args[2]
        flags = Bool[]

        flag1 = isequal(term3.op, *)
        push!(flags, flag1)

        flag2 = length(term3.args) == 2
        push!(flags, flag2)

        y = term3.args[2]
        # TODO add check across terms so that all dependents are equal, all names are equal
        flag3 = typeof(y) == Variable
        push!(flags, flag3)

        flag4 = length(y.dependents) >= 1
        push!(flags, flag4)

        # Add logic for order 0
        n = term3.args[1].value
        flag5 = n isa Int && n >= 1
        push!(flags, flag5)

        return all(flags)
end

# Todo - cross check dependent variable is always the same
# This lets us use different dependent variabls as long as they are consistent
function check_dependents(eq :: Equation, y ::Type{LaguerreAbstractEquation})

end

# todo - cross check dependent names is always the same
function check_names(eq :: Equation, y :: Type{LaguerreAbstractEquation})

end

# Hermite
# y'' - 2xy' + λy = 0
# eq = D2(y) - 2*x*D(y) + λ*y ~ 0
function check_term1(eq :: Equation, y :: Type{HermiteAbstractEquation})
        term1 = eq.lhs.args[1].args[1]
        flags = Bool[]

        DD2 = term1
        flag1 = typeof(DD2.op) == Differential
        push!(flags, flag1)

        flag2 = typeof(DD2.args[1].op) == Differential
        push!(flags, flag2)

        y = term1.args[1].args[1]
        flag3 = length(y.dependents) >= 1
        push!(flags, flag3)

        flag4 = typeof(y) == Variable
        push!(flags, flag4)

        return all(flags)
end
function check_term2(eq :: Equation, y :: Type{HermiteAbstractEquation})
        term2 = eq.lhs.args[1].args[2]
        flags = Bool[]

        flag1 = length(term2.args) == 2
        push!(flags, flag1)

        # check that it is indeed the derivative we want
        flag2 = typeof(term2.args[2].op) == Differential
        push!(flags, flag2)

        y = term2.args[2].args[1] # check names and dependents
        flag3 = length(y.dependents) == 1

        flag4 = typeof(y) == Variable
        push!(flags,flag4)

        # check it is an algebraic form we want
        flag5 = length(term2.args[1].args) == 2
        push!(flags, flag5)

        x = term2.args[1].args[2] # check names and dependents
        flag6 = typeof(x) == Variable
        push!(flags, flag6)

        # alpha + 1 TODO assume alpha is more than 0
        flag7 = term2.args[1].args[1].value == 2
        push!(flags, flag7)

        flag8 = eq.lhs.args[1].op == -
        push!(flags, flag8)

        return all(flags)
end
function check_term3(eq :: Equation, y :: Type{HermiteAbstractEquation})
        term3 = eq.lhs.args[2]
        flags = Bool[]

        flag1 = isequal(term3.op, *)
        push!(flags, flag1)

        flag2 = length(term3.args) == 2
        push!(flags, flag2)

        y = term3.args[2]
        # TODO add check across terms so that all dependents are equal, all names are equal
        flag3 = typeof(y) == Variable
        push!(flags, flag3)

        flag4 = length(y.dependents) >= 1
        push!(flags, flag4)

        # Add logic for order 0
        n = term3.args[1].value
        flag5 = n isa Int && n >= 1
        push!(flags, flag5)

        return all(flags)
end

# TODO - cross check dependent variable is always the same
# This lets us use different dependent variabls as long as they are consistent
function check_dependents(eq :: Equation, y ::Type{HermiteAbstractEquation})

end

# TODO - cross check dependent names is always the same
function check_names(eq :: Equation, y :: Type{HermiteAbstractEquation})

end

# Legendre
# (1 - x^2) * y'' -2xy' + l(l+1)y = 0
# eqtest1 = (1-x^2) * D2(y) - 2*x*D(y) + l*(l+1)*y ~ 0
function check_term1(eq :: Equation, y :: Type{LegendreAbstractEquation})
        term1 = eq.lhs.args[1].args[1]
        flags = Bool[]

        flag1 = isequal(term1.op, *)
        push!(flags, flag1)

        flag2 = length(term1.args) == 2
        push!(flags, flag2)

        flag3 = typeof(term1.args[1]) == Operation
        push!(flags, flag3)

        parens = term1.args[1]
        flag4 = isequal(parens.op, -)
        push!(flags, flag4)

        flag5 = parens.args[1].value == 1
        push!(flags, flag5)

        xsq = parens.args[2]
        flag6 = length(xsq.args) == 2
        push!(flags, flag6)

        flag7 = isequal(xsq.op, ^)
        push!(flags, flag7)

        flag8 = xsq.args[2].value == 2
        push!(flags, flag8)

        DD2 = term1.args[2]
        flag9 = typeof(DD2.op) == Differential
        push!(flags, flag9)

        # Time as a dependent?
        y = DD2.args[1].args[1]
        flag10 = typeof(DD2.args[1].op) == Differential
        push!(flags, flag10)

        flag11 = length(y.dependents) >= 1
        push!(flags, flag11)

        return all(flags)
end
function check_term2(eq :: Equation, y :: Type{LegendreAbstractEquation})
        term2 = eq.lhs.args[1].args[2]
        flags = Bool[]

        flag1 = length(term2.args) == 2
        push!(flags, flag1)

        # check that it is indeed the derivative we want
        flag2 = typeof(term2.args[2].op) == Differential
        push!(flags, flag2)

        y = term2.args[2].args[1] # check names and dependents
        flag3 = length(y.dependents) == 1

        flag4 = typeof(y) == Variable
        push!(flags,flag4)

        # check it is an algebraic form we want
        flag5 = length(term2.args[1].args) == 2
        push!(flags, flag5)

        x = term2.args[1].args[2] # check names and dependents
        flag6 = typeof(x) == Variable
        push!(flags, flag6)

        flag8 = eq.lhs.args[1].op == -
        push!(flags, flag8)

        return all(flags)
end
function check_term3(eq :: Equation, y :: Type{LegendreAbstractEquation})
        term3 = eq.lhs.args[2]
        flags = Bool[]

        flag1 = isequal(term3.op, *)
        push!(flags, flag1)

        flag2 = length(term3.args) == 2
        push!(flags, flag2)

        y = term3.args[2]
        # TODO add check across terms so that all dependents are equal, all names are equal
        flag3 = typeof(y) == Variable
        push!(flags, flag3)

        flag4 = length(y.dependents) >= 1
        push!(flags, flag4)

        # Add logic for order 0
        l = term3.args[1].value
        flag5 = l isa Int && l >= 1 && llplus1checker(l)
        push!(flags, flag5)

        return all(flags)
end

function llplus1checker(l)
        l1 = isqrt(l)
        l2 = l1 + 1
        l1 * l2 == l ? true : false
end

# TODO - cross check dependent variable is always the same
# This lets us use different dependent variabls as long as they are consistent
function check_dependents(eq :: Equation, y ::Type{LegendreAbstractEquation})

end

# TODO - cross check dependent names is always the same
function check_names(eq :: Equation, y :: Type{LegendreAbstractEquation})

end

# Analytic Extraction

# Laguerre
# xy'' + (α + 1 - x)* y' + ny = 0
# eqlag = x*D2(y) + (1 + 1 - x) * D(y) + 1*y ~ 0
function extract(eq :: Equation, y :: Type{LaguerreAbstractEquation})
        n = eq.lhs.args[2].args[1].value
        α = eq.lhs.args[1].args[2].args[1].args[1].value - 1
        LaguerreEquation{α, n}()
end

# Hermite
# y'' - 2xy' + λy = 0
# eqher = D2(y) - 2*x*D(y) + λ*y ~ 0
function extract(eq :: Equation, y :: Type{HermiteAbstractEquation})
        n = eq.lhs.args[2].args[1].value
        HermiteEquation{n}()
end

# Legendre
# (1 - x^2) * y'' -2xy' + l(l+1)y = 0
# eqleg = (1-x^2) * D2(y) - 2*x*D(y) + l*(l+1)*y ~ 0
function extract(eq :: Equation, y :: Type{LegendreAbstractEquation})
        lp = eq.lhs.args[2].args[1].value
        n = isqrt(lp)
        LegendreEquation{n}()
end

function build(x :: LaguerreEquation{α,n}) where {α,n}
        return Ls[n + 1]
end

function build(x :: HermiteEquation{n}) where n
        return Hs[n + 1]
end

function build(x :: LegendreEquation{n}) where n
        Ps[n + 1]
end

function transform(x :: Equation, y ::Type{T}) where T<:AbstractEquation
        typ = extract(x,y)
        build(typ)
end

function analyze(eq :: Equation)
        try
                if isit(eq, LaguerreAbstractEquation)
                        return transform(eq, LaguerreAbstractEquation)
                end
        catch
        end

        try
                if isit(eq, HermiteAbstractEquation)
                        return transform(eq, HermiteAbstractEquation)
                end
        catch
        end

        try
                if isit(eq, LegendreAbstractEquation)
                        return transform(eq, LegendreAbstractEquation)
                end
        catch
        end

        throw("No analytic solution detected.")
end

export isit, extract, build, transform, analyze
export LaguerreAbstractEquation, LaguerreEquation
export HermiteAbstractEquation, HermiteEquation
export LegendreAbstractEquation, LegendreEquation

L0(x) = 1.0
    L1(x) = -x + 1.0
    L2(x) = .5 * (x^2 - 4x + 2)
    L3(x) = (1/6)*(-x^3 +9x^2 -18x + 6)
    L4(x) = (1/24)*(x^4 - 16x^3 + 72x^2 - 96x + 24)
    L5(x) = (1/120)*(-x^5 + 25x^4 -200x^3 + 600x^2 - 600x + 120)
    L6(x) = (1/720)*(x^6 - 36x^5 + 450x^4 - 2400x^3 + 5400x^2 - 4320x + 720)
    L7(x) = (1/5040)*(-x^7 + 49x^6 - 882x^5 + 7350x^4 - 29_400x^3 + 52_920x^2 - 35_280x + 5040)
    L8(x) = (1/40320)*(x^8 - 64x^7 + 1568x^6 - 18_816x^5 + 117_600x^4 - 376_320x^3 + 564_480x^2 - 322_560x + 40320)
    L9(x) = (1/362880)*(-x^9 + 81x^8 - 2592x^7 + 42_336x^6 - 381_024x^5 + 1_905_120x^4 - 5_080_320x^3 + 6_531_840x^2 - 3_265_920x + 362_880)
    L10(x) = (1/3628800)*(x^10 - 100x^9 + 4_050x^8 - 86_400x^7 + 1_058_400x^6 - 7_620_480x^5 + 31_752_000x^4 - 72_576_000x^3  + 81_648_000x^2 - 36_288_000x + 3_628_800)
    L11(x) = (1/39916800)*(-x^11 + 121x^10 - 6_050x^9 + 163_350x^8 - 2_613_600x^7 + 25_613_280x^6 - 153_679_680x^5 + 548_856_000x^4 - 1_097_712_000x^3 + 1_097_712_000x^2 - 439_084_800x + 39_916_800)
    L12(x) = (1/479001600)*(x^12 - 144x^11 + 8_712x^10 - 290_400x^9 + 5_880_600x^8 - 75_271_680x^7 + 614_718_720x^6 - 3_161_410_560x^5 + 9_879_408_000x^4 - 17_563_392_000x^3 + 15_807_052_800x^2 - 5_748_019_200x + 479_001_600)

    Ls = [L0, L1, L2, L3, L4, L5, L6, L7, L8, L9, L10, L11, L12]

    # Credit to @yingboma
    function laguerre(n, x)

        n == 0 ? (return 1) : nothing

        l0, l1 = 1, 1-x
        for k in 1:n-1
            l0, l1 = l1, ( (2k+1-x)*l1 - k*l0 )/(k+1)
        end
        l1
    end

H0(x) = 1.0
    H1(x) = 2x
    H2(x) = -2 + 4x^2
    H3(x) = -12x + 8x^3
    H4(x) = 12 - 48x^2 + 16x^4
    H5(x) = 120x -160x^3 + 32x^5
    H6(x) = -120 + 720x^2 - 480x^4 + 64x^6
    H7(x) = -1680x + 3360x^3 - 1344x^5 + 128x^7
    H8(x) =1680 - 13440x^2 + 13440x^4 - 3584x^6 + 256x^8
    H9(x) =30240x - 80640x^3 + 48384x^5 - 9216x^7 + 512x^9
    H10(x) =-30240 + 302400x^2 - 403200x^4 + 161280x^6 - 23040x^8 + 1024x^10
    H11(x) =-665280x + 2217600x^3 - 1774080x^5 + 506880x^7 - 56320x^9 + 2048x^11
    H12(x) =665280 - 7983360x^2 + 13305600x^4 - 7096320x^6 + 1520640x^8 - 135168x^10 + 4096x^12

    Hs = [H0, H1, H2, H3, H4, H5, H6, H7, H8, H9, H10, H11, H12]

    function hermite(x,n)
        if n == 0
            return 1
        elseif n == 1
            return 2x
        else
            return 2x*hermite(x,n-1) - 2*(n-1)*hermite(x,n-2)
        end
    end

P0(x) = 1
    P1(x) =  x
    P2(x) = .5*(3x^2 -1)
    P3(x) = .5*(5x^3 -3x)
    P4(x) = (1/8)*(35x^4 - 30x^2 + 3)
    P5(x) = (1/8)*(63x^5 - 70x^3 + 15x)
    P6(x) = (1/16)*(231x^6 - 315x^4 + 105x^2 - 5)
    P7(x) = (1/16)*(429x^7 - 693x^5 + 315x^3 - 35x)
    P8(x) = (1/128)*(6435x^8 - 12012x^6 + 6930x^4 - 1260x^2 + 35)
    P9(x) = (1/128)*(12155x^9 - 25740x^7 + 18018x^5 - 4620x^3 + 315x)
    P10(x) = (1/256)*(46189x^10 - 109395x^8 + 90090x^6- 30030x^4 + 3465x^2 - 63)
    P11(x) = (1/256)*(2^11)*(705432x^11 - 1847560x^9 + 1750320x^7 - 720720x^5 + 120120x^3 - 5544x)
    P12(x) = (1/512)*(2^12)*(2704156x^12 - 7759752x^10 + 8314020x^8 - 4084080x^6 + 900900x^4 - 72072x^2 + 924)

    Ps = [P0, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12]
    testPs = [P2, P3, P9, P10]

    function legendre(x, n)
        if n == 0
            return 1
        elseif n == 1
            return x
        else
            return ((2n-1)*x*legendre(x,n-1) - (n - 1)*legendre(x, n-2))/n
        end
    end

export Ls, Hs, Ps
end # module


text:
exact_energies = Dict(pc"He" => [-2.8616800, [o"1s" => -0.91795555]],
                      pc"Be" => [-14.573023, [o"1s" => -4.7326698, o"2s" => -0.3092695]],

                      pc"Ne" => [-128.54710, [o"1s" => -32.7724455, o"2s" => -1.93039095, o"2p" => -0.85040965]],
                      PseudoPotentials.NeonHF => [-34.709465, [o"2s" => -1.93039095, o"2p" => -0.85040965]],
                      PseudoPotentials.NeonWB => [-34.708785, [o"2s" => -1.93039095, o"2p" => -0.85040965]],

                      pc"Mg" => [-199.61463, [o"1s" => -49.0317255, o"2s" => -3.767718, o"2p" => -2.2822236, o"3s" => -0.25305275]],

                      pc"Ar" => [-526.81751, [o"1s" => -237.22070/2, o"2s" => -24.644306/2, o"2p" => -19.142932/2,
                                             o"3s" => -2.5547063/2, o"3p" => -1.1820348/2]],
                      PseudoPotentials.ArgonHF => [-20.849867, [o"3s" => -2.5547063/2, o"3p" => -1.1820348/2]],
                      PseudoPotentials.ArgonWB => [-20.884584, [o"3s" => -2.5547063/2, o"3p" => -1.1820348/2]],

                      pc"Ca" => [-676.75818, [o"1s" => -298.72744/2, o"2s" => -33.645481/2, o"2p" => -27.258531/2,
                                          o"3s" => -4.4907488/2, o"3p" => -2.6814114/2, o"4s" => -0.3910594/2]],
                      pc"Zn" => [-1777.8481, [o"1s" => -706.60909/2, o"2s" => -88.723452/2, o"2p" => -77.849691/2,
                                             o"3s" => -11.275642/2, o"3p" => -7.6787581/2, o"3d" => -1.5650843/2,
                                             o"4s" => -0.5850141/2]],
                      pc"Kr" => [-2752.0550, [o"1s" => -1040.3309/2, o"2s" => -139.80617/2, o"2p" => -126.01957/2,
                                             o"3s" => -21.698934/2, o"3p" => -16.663004/2, o"3d" => -7.6504697/2,
                                             o"4s" => -2.3058703/2, o"4p" => -1.0483734/2]],
                      pc"Xe" => [-7232.1384, [o"1s" => -2448.7956/2, o"2s" => -378.68024/2, o"2p" => -355.56490/2,
                                             o"3s" => -80.351324/2, o"3p" => -70.443321/2, o"3d" => -52.237736/2,
                                             o"4s" => -15.712602/2, o"4p" => -12.016674, o"4d" => -5.5557597,
                                             o"5s" => -1.8888279/2, o"5p" => -0.9145793/2]],

                      PseudoPotentials.XenonHF => [-14.989100, [o"5s" => -1.8888279/2, o"5p" => -0.9145793/2]],
                      PseudoPotentials.XenonWB => [-15.277055, [o"5s" => -1.8888279/2, o"5p" => -0.9145793/2]],
                      PseudoPotentials.XenonDF2c => [-328.74543, [o"5s" => -1.0097, ro"5p-" => -0.4915, ro"5p" => -0.4398]])

function shift_unit(u::U, d) where {U<:Unitful.Unit}
    d = Int(3floor(Int, d/3))
    iszero(d) && return u
    for tens = Unitful.tens(u) .+ (d:(-3*sign(d)):0)
        haskey(Unitful.prefixdict, tens) && return U(tens, u.power)
    end
    u
end

function shift_unit(u::Unitful.FreeUnits, d)
    tu = typeof(u)
    us,ds,a = tu.parameters

    uu = shift_unit(us[1], d)

    Unitful.FreeUnits{(uu,us[2:end]...), ds, a}()
end

function si_round(q::Quantity; fspec="{1:+9.4f} {2:s}")
    v,u = ustrip(q), unit(q)
    if !iszero(v)
        u = shift_unit(u, log10(abs(v)))
        q = u(q)
    end
    format(fspec, ustrip(q), unit(q))
end

function energy_errors(fock, exact_energies, Δ, δ)
    atom = fock.quantum_system

    Eexact,orbExact = exact_energies
    exact_energies = if first(atom.orbitals) isa SpinOrbital
        vcat(Eexact, [repeat([E], degeneracy(o)) for (o,E) in orbExact]...)
    else
        vcat(Eexact, [E for (o,E) in orbExact]...)
    end

    orbital_refs = unique(nonrelorbital.(first.(orbExact)))

    H = zeros(1,1)
    Etot = SCF.energy_matrix!(H, fock)[1,1]
    energies = vcat(Etot, collect(SCF.energy(fock.equations.equations[i])/degeneracy(o)
                                  for (i,o) in enumerate(atom.orbitals)
                                  if nonrelorbital(Atoms.getspatialorb(o)) ∈ orbital_refs))
    errors = energies - exact_energies

    labels = vcat("Total", string.(collect(o for o in atom.orbitals
                                           if nonrelorbital(Atoms.getspatialorb(o)) ∈ orbital_refs)))

    Ha_formatter = (v,i) -> si_round(v*u"hartree")
    eV_formatter = (v,i) -> si_round(v*u"eV")

    pretty_table([labels exact_energies energies errors 27.211energies 27.211errors errors./abs.(exact_energies)],
                 ["", "HF limit", "Energy", "Δ", "Energy", "Δ", "Δrel"],
                 formatter=Dict(
                     2 => Ha_formatter,
                     3 => Ha_formatter,
                     4 => Ha_formatter,
                     5 => eV_formatter,
                     6 => eV_formatter,
                     7 => (v,i) -> si_round(100v*u"percent")
                 ),
                 highlighters=(Highlighter((v,i,j) -> abs(v[i,7])>0.2, foreground=:red, bold=true),
                               Highlighter((v,i,j) -> abs(v[i,7])>0.01, foreground=:yellow, bold=true),))

    @test abs(errors[1]) < Δ
    @test all(abs.(errors[2:end]) .< δ)
end

function atom_calc(nucleus::AbstractPotential, grid_type::Symbol, rₘₐₓ, ρ,
                   Δ, δ;
                   config_transform=identity, kwargs...)
    # If we're using a pseudopotential, we don't really need higher
    # order of the FEDVR basis functions in the first finite element,
    # since the orbitals almost vanish there.
    R,r = get_atom_grid(grid_type, rₘₐₓ, ρ, nucleus,
                        amend_order=nucleus isa PointCharge)

    atom = Atom(R, [spin_configurations(config_transform(ground_state(nucleus)))[1]],
                nucleus, eltype(R))

    fock = Fock(atom)

    optimize!(fock; kwargs...)
    energy_errors(fock, exact_energies[nucleus], Δ, δ)
end

@testset "Calculation accuracy" begin
    @testset "Helium" begin
        @testset "$(orb_type) orbitals" for (orb_type,config_transform) in [
            ("Non-relativistic", identity),
            ("Relativistic", relconfigurations)
        ]
            @testset "$(grid_type)" for (grid_type,Δ,δ) in [(:fedvr,6e-9,7e-9),
                                                            (:fd,4e-3,4e-3)]
                atom_calc(pc"He", grid_type, 10.0, 0.1, Δ, δ, ω=0.9,
                          config_transform=config_transform)
            end
        end
    end

    @testset "Beryllium" begin
        @testset "$(grid_type)" for (grid_type,ρ,Δ,δ) in [(:fedvr,0.2,6e-7,2e-6),
                                                          (:fd,0.05,0.02,0.009)]
            atom_calc(pc"Be", grid_type, 15.0, ρ, Δ, δ, ω=0.9, scf_method=:arnoldi)
        end
    end

    @testset "Neon" begin
        @testset "$nucleus" for (nucleus,Δ,δ) in [(pc"Ne",0.005,0.005),
                                                  (PseudoPotentials.NeonHF,0.2,0.05)]
            atom_calc(nucleus, :fedvr, 10, 0.2, Δ, δ,
                      ω=0.999, ωmax=1.0-1e-3, scf_method=:arnoldi)
        end
    end

    @testset "Xenon" begin
        @testset "$nucleus" for (nucleus,grid_type,ρ,Δ,δ,config_transform) in [
            (PseudoPotentials.XenonHF,:fd,0.1,4e-3,4e-3,identity),
            (PseudoPotentials.XenonDF2c,:fd,0.1,0.2,1e-2,relconfigurations)
        ]
            atom_calc(nucleus, grid_type, 7.0, ρ, Δ, δ,
                      ω=0.999, ωmax=1.0-1e-3,
                      config_transform=config_transform, scf_method=:arnoldi)
        end
    end
end


text:
"""
genetrap.jl

Fit GRSM models to live cell and smFISH data
"""

"""
Gene information
"""
const genes = ["CANX";"DNAJC5";"ERRFI1";"KPNB1";"MYH9";"Rab7a";"RHOA";"RPAP3";"Sec16A";"SLC2A1"]
const genelength = Dict([("Sec16A", 42960);("SLC2A1", 33802);("ERRFI1", 14615);("RHOA", 52948);("KPNB1", 33730);("MYH9", 106741);("DNAJC5", 40930);("CANX", 32710);("Rab7a", 88663);("RPAP3", 44130)])
const MS2end = Dict([("Sec16A", 5220);("SLC2A1", 26001);("ERRFI1", 5324);("RHOA", 51109);("KPNB1", 24000);("MYH9", 71998);("DNAJC5", 14857);("CANX", 4861);("Rab7a", 83257);("RPAP3", 38610)])
const halflife = Dict([("CANX", 50.),("DNAJC5", 5.),("ERRFI1", 1.35),("KPNB1", 9.),("MYH9", 10.),("Rab7a", 50.),("RHOA", 50.),("RPAP3", 7.5),("Sec16A", 8.),("SLC2A1", 5.)])

"""
genetrap()

Load data, model and option structures for metropolis-hastings fit
root is the folder containing data and results

FISH counts is divided by tempfish to adjust relative weights of data
set tempfish = 0 to equalize FISH and live cell counts
"""
function genetrap(root::String,infolder::String,rinchar::String,label::String,gene::String,G::Int,R::Int,nalleles::Int,type::String,maxtime::Float64,samplesteps::Int,method=1,temp=1.,tempfish=1.,warmupsteps=0,annealsteps=0)
    r = readrates_genetrap(root,infolder,rinchar,gene,"$G$R",type)
    println(r)
    genetrap(root,r,label,gene,G,R,nalleles,type,maxtime,samplesteps,method,temp,tempfish,annealsteps,warmupsteps)
end

function genetrap(root,r,label::String,gene::String,G::Int,R::Int,nalleles::Int,type::String,maxtime::Float64,samplesteps::Int,method,temp=1.,tempfish=1.,warmupsteps=0,annealsteps=0)
    data = data_genetrap(root,label,gene,tempfish)
    model = model_genetrap(r,gene,G,R,nalleles,type,method)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp)
    return data,model,options
end

function data_genetrap(root,label,gene,tempfish=1.)
    LC = readLCPDF_genetrap(root,gene)
    if tempfish == 0
        counts = Int(div(sum(LC[:,2]+LC[:,3]),2))
        println(counts)
        # set FISH counts to mean of live cell counts
        histFISH = readFISH_genetrap(root,gene,counts)
    else
        histFISH = readFISH_genetrap(root,gene,tempfish)
    end
    RNALiveCellData(label,gene,LC[:,1],LC[:,3],LC[:,2],length(histFISH),histFISH)
end

"""
model_genetrap

load model structure
"""
function model_genetrap(r,gene::String,G::Int,R::Int,nalleles::Int,type::String,method,propcv=[.02*ones(2*(G-1));1e-4;1e-4*ones(R);.02*ones(R)],fittedparam=collect(1:num_rates(G,R)-1),Gprior=(.01,.5),Sprior=(.1,.5),Rcv=1e-2)
    rm,rcv = prior_rates_genetrap(G,R,gene,Gprior,Sprior,Rcv)
    # d = priordistributionLogNormal_genetrap(rm,rcv,G,R)
    d = priordistributionGamma_genetrap(rm,rcv,G,R)
    if r == 0
        r = rm
    else
        r = r[:,1]
    end
    GRSMmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,R,nalleles,type,r,d,propcv,fittedparam,method)
    # GRSMmodel(G,R,nalleles,type,r,d,prop,fittedparam,method)
end

"""
setpriorrate(G,R,gene,Gprior::Tuple,Sprior::Tuple,Rcv::Float64)
Set prior distribution for mean and cv of rates
"""
function prior_rates_genetrap(G,R,gene,Gprior::Tuple,Sprior::Tuple,Rcv::Float64)
    n = G-1
    rm = [fill(Gprior[1],2*n);2000/(genelength[gene]-MS2end[gene]);fill(2000/MS2end[gene]*R,R);fill(Sprior[1], R);log(2.)/(60 .* halflife[gene])]
    rcv = [fill(Gprior[2],2*n);Rcv;Rcv;fill(Sprior[2], R);Rcv]
    return rm,rcv
end

"""
priordistributionLogNormal_genetrap(r,cv,G,R)

LogNormal Prior distribution
"""
function priordistributionLogNormal_genetrap(r,cv,G,R)
    sigma = sigmalognormal(cv)
    d = []
    j = 1
    #  G rates
    for i in Grange(G)
        push!(d,Distributions.LogNormal(log(r[i]),sigma[j]))
        j += 1
    end
    # initiation rate
    i = initiation(G)
    push!(d,Distributions.LogNormal(log(r[i]),sigma[j]))
    j += 1
    # Step rates are bounded by length of insertion site to end of gene, i.e sum of reciprocal rates is bounded
    t = sum(1 ./ r[Rrange(G,R)])
    push!(d,Distributions.LogNormal(-log(t),sigma[j]))
    j += 1
    # priors for splice rates
    rs = 0.
    for i in Srange(G,R)
        rs += r[i]
        push!(d,Distributions.LogNormal(log(rs),sigma[j]))
        j += 1
    end
    return d
end

"""
priordistributionGamma_genetrap(r,cv,G,R)

Gamma Prior distribution
"""
function priordistributionGamma_genetrap(rm,cv,G,R)
    n = G-1
    zet = R
    d = []
    rsig = cv .^ 2
    j = 1
    # priors for G rates
    for i in Grange(G)
        theta = rsig[j]*rm[i]
        alpha = 1/rsig[j]
        push!(d,Distributions.Gamma(alpha,theta))
        j += 1
    end
    # prior for nu1 is upper bounded by rm, i.e. distance is bounded by distance to insertion site
    theta = rsig[j]*rm[initiation(G)]
    alpha = 1/rsig[j]
    push!(d,Distributions.Gamma(alpha,theta))
    j += 1
    # priors for nu rates are bounded by length of insertion site to end of gene, i.e sum of reciprocal rates is bounded
    tm = sum(1 ./ rm[Rrange(G,R)])
    # for i in 2*n + 2 : 2*n + zet + 1
    #     tm += 1/rm[i]
    # end
    theta = rsig[j]/tm
    alpha = 1/rsig[j]
    push!(d,Distributions.Gamma(alpha,theta))
    j += 1
    # priors for ejection rates
    rs = 0
    for i in Srange(G,R)
        rs += rm[i]
        theta = rsig[j]*rs
        alpha = 1/rsig[j]
        push!(d,Distributions.Gamma(alpha,theta))
        j += 1
    end
    return d
end

"""
readLCPDF_genetrap(root,gene)

Read in dwell time PDF
"""
function readLCPDF_genetrap(root,gene)
    infile = joinpath(root,"DwellTimePDF/$(gene)_PDF.csv")
    if isfile(infile)
        LC = readdlm(infile,',')
        x = truncate_histogram(LC[:,2],.999,1000)
        LC[1:length(x),:]
    else
        println("No data for gene: ", gene)
    end
end

"""
readFISH_genetrap(root,gene,temp=1.,clone=true)

Read and assemble smFISH data
"""
function readFISH_genetrap(root::String,gene::String,temp::Float64=1.,clone=true)
    counts = countsFISH_genetrap(root,gene,clone)
    readFISH_genetrap(root,gene,Int(div(counts,temp)),clone)
end

function readFISH_genetrap(root::String,gene::String,counts::Int,clone=true)
    fishfile=joinpath(root,"Fish_4_Genes/$(gene)_steady_state_mRNA.csv")
    col = clone ? 3 : 2
    # Get smFISH RNA histograms
    if isfile(fishfile)
        x = readdlm(fishfile,',')[3:end,col]
        x = x[x .!= ""]
        x = truncate_histogram(x,.99,1000)
        x /= sum(x)
        x *= counts
        return x
    else
        println("No smFISH data for gene: ", gene)
    end
    nothing
end

"""
countsFISH_genetrap(root,gene,x,counts,clone=true)

read in smFISH counts
"""
function countsFISH_genetrap(root,gene,clone=true)
    countfile = joinpath(root,"counts/total_counts.csv")
    wttext = "WT_"
    clonetext = "_clone"
    if isfile(countfile)
        countsdata = readdlm(countfile,',')[:,:]
        if clone
            counts=(countsdata[findall(uppercase("$(gene)$clonetext") .== uppercase.(countsdata[:,1])),2])[1]
        else
            counts=(countsdata[findall("$wttext$(gene)" .== countsdata[:,1]),2])[1]
        end
    else
        counts = 1000
    end
    return counts
end

"""
readrates_genetrap(infolder::String,rinchar::String,gene::String,model::String,type::String)

Read in initial rates from previous runs
"""
function readrates_genetrap(root::String,infolder::String,rinchar::String,gene::String,model::String,type::String)
    if rinchar == "rml" || rinchar == "rlast"
        rskip = 1
        rstart = 1
    elseif rinchar == "rmean"
        rskip = 2
        rstart = 1
    elseif rinchar == "rquant"
        rskip = 3
        rstart = 2
    end
    if type == "offeject" || type == "on"
        type = ""
    end
    # infile = "$resultpath/$infolder/$(gene)/$rinchar$model$type$txtstr"
    infile = getratefolder_genetrap(root,infolder,rinchar,gene,model,type)
    print(gene," ",model," ",type)
    readrates_genetrap(infile,rstart,rskip)
end

function readrates_genetrap(infile::String,rstart::Int,rskip::Int)
    if isfile(infile) && ~isempty(read(infile))
        r = readdlm(infile)
        r = r[:,rstart:rskip:end]
        println("")
        return r
    else
        println(" no prior")
        return 0
    end
    nothing
end

function getratefolder_genetrap(root::String,infolder::String,rinchar::String,gene::String,model::String,type::String)
    results = joinpath(root,"Results")
    infolder = joinpath(infolder,gene)
    infolder = joinpath(infolder,rinchar * model * type * txtstr)
    joinpath(results,infolder)

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:
export runOEDTypeA

"""
function runOEDTypeA

Input:

	param
	Xtrue
	L
	n
	theta
	Alphas
	model
	IPmode
"""
function runOEDTypeA(OEDparam::Array{OEDTypeAParam},Xtrue::Array{Float64},L::AbstractArray,Alphas,modFun::Function,IPmode::Symbol,m0::Vector{Float64};Ce=[],ce=[],Ci=[],ci=[],pcgMaxIter=10,maxIter=20,maxStep=1.0,exName="")

	ids = mod.(0:size(Xtrue,2)-1,nworkers())+1

	println("\n\n=== setup forward problem ===")
	if IPmode==:unconstr
		pFor = Array{UnconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
			pFor[j] = UnconstrainedOEDParam(OEDparam[j],L,zeros(0,0),[])
		end
	elseif IPmode==:eqconstr
		if isempty(Ce) || isempty(ce)
			error("equality constraints must be provided. To run unconstrained OED, use IPmode=:unconstr.")
		end
		pFor = Array{EqconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
			pFor[j] = EqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],zeros(0,0),[])
		end
	elseif IPmode==:nonnegconstr
		I = speye(size(Xtrue,1))
		pFor = Array{IneqconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
				pFor[j] = getIneqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],I,0*Xtrue[:,ids.==j],
				          ftolIPQP=1e-3,mutolIPQP=1e-3,outIPQP=0)
		end

	elseif IPmode==:boxconstr
		if isempty(Ci) || isempty(ci)
			error("inequality constraints must be provided. To run equality constrained OED, use IPmode=:eqconstr.")
		end
		pFor = Array{IneqconstrainedOEDParam}(nworkers())

		for j=1:nworkers()
		pFor[j] = 	getIneqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],Ci,ci[:,ids.==j],
		ftolIPQP=1e-3,mutolIPQP=1e-3,outIPQP=0)
		end
	else
		error("Invalid IPmode = $IPmode")
	end
	display(pFor[1])

	np = length(m0)
	A  = copy(pFor[1].Design.A)
	D  = zeros(size(A,1),size(Xtrue,2))
	for k=1:nworkers()
		D[:,ids.==k] = pFor[k].Design.D
	end

	for k=1:length(Alphas)
		alpha = Alphas[k]
		println("\n\n==== processing alpha=($alpha) which is $(k) out of $(length(Alphas)) ====\n\n")
		println("\n\n=== setup first optimization to identify non-zeros ===")
		if nworkers()>1
			pMis  = Array{RemoteChannel}(nworkers())
			for j=1:nworkers()
				pFor[j].Design.A = A
				pFor[j].Design.D = D[:,ids.==j]
				pMis[j] = initRemoteChannel(getMisfitParam,workers()[j], pFor[j],ones(size(Xtrue[:,ids.==j])),Xtrue[:,ids.==j],SSDFun)
			end
		else
			pFor[1].Design.A = A
			pFor[1].Design.D = D
			pMis        = getMisfitParam(pFor[1],ones(size(Xtrue)),Xtrue,SSDFun);
		end

		# make sure the full size problem is in there
		Minv        = getRegularMesh([0  1 0 1],[12 12]);
		pInv        = getInverseParam(Minv,modFun,l1Reg,alpha,zeros(np),zeros(np),Inf*ones(np));
		pInv.pcgMaxIter = pcgMaxIter
		pInv.maxIter    = maxIter
		pInv.maxStep    = maxStep
		display(pInv)

		println("\n\n=== first inversion started $(now()) ===")
		tic()
		mOpt1,Xo,flag1,His1 = projGN(copy(m0),pInv,pMis);
		time1 = toc()
		wOpt1, = modFun(mOpt1)

		Xopt1 = 0*Xtrue
		if nworkers() == 1
			Xopt1 = Xo
		else
			for jj=1:nworkers()
				Xopt1[:,ids.==jj] = fetch(Xo[jj])
			end
		end

		println("model:   nnz=$(countnz(mOpt1)), sparsity=$(countnz(mOpt1)/length(mOpt1))")
		println("weights: nnz=$(countnz(wOpt1)), sparsity=$(countnz(wOpt1)/length(wOpt1))")

		println("\n\n=== setup second optimization to tune nonzero weights ===")
		if nworkers()>1
			pMis  = Array{RemoteChannel}(nworkers())
			for j=1:nworkers()
				pFor[j].Design.A = A[wOpt1.>0,:]
				pFor[j].Design.D = D[wOpt1.>0,ids.==j]
				pMis[j] = initRemoteChannel(getMisfitParam,workers()[j], pFor[j],ones(size(Xtrue[:,ids.==j])),Xtrue[:,ids.==j],SSDFun)
			end
		else
			pFor[1].Design.A = A[wOpt1.>0,:]
			pFor[1].Design.D = D[wOpt1.>0,:]
			pMis        = getMisfitParam(pFor[1],ones(size(Xtrue)),Xtrue,SSDFun);
		end

		nnzM = countnz(mOpt1)
		pInv.boundsLow  = zeros(nnzM)
		pInv.boundsHigh = Inf*ones(nnzM)
		pInv.mref       = zeros(nnzM)
		pInv.alpha      = 0.0
		pInv.HesPrec.applyPrec = (a,b,x,c) -> x
		display(pInv)

		println("\n\n=== second inversion started $(now()) ===")
		tic()
		mOpt2,Xo,flag2,His2 = projGN(mOpt1[mOpt1.>0],pInv,pMis);
		time2 = toc()
		wOpt2, = modFun(mOpt2)
		Xopt2 = 0*Xtrue
		if nworkers() == 1
			Xopt2 = Xo
		else
			for jj=1:nworkers()
				Xopt2[:,ids.==jj] = fetch(Xo[jj])
			end
		end

		if !isempty(exName)
			fileName = @sprintf "%s-%s-alpha-%2.2e-nnz-%d.mat" exName IPmode alpha countnz(mOpt2)
			logName = @sprintf "%s-%s-alpha-%2.2e-nnz-%d.log" exName IPmode alpha countnz(mOpt2)
			mv("jInv.out",logName,remove_destination=true)
			println("=== save results in $fileName===")
			His1.Dc = []; His2.Dc = []
			results = Dict("mOpt1"=>mOpt1, "Xopt1"=>Xopt1,"flag1"=>flag1,"His1"=>His1,"time1"=>time1,"wOpt1"=>wOpt1,
			               "mOpt2"=>mOpt2, "Xopt2"=>Xopt2,"flag2"=>flag2,"His2"=>His2,"time2"=>time2,"wOpt2"=>wOpt2)
			matwrite(fileName,results)
		end
		m0 = mOpt1;
	end
	return m0
end


text:
# A basic example of using the BayesianMixtures package.

using BayesianMixtures
B = BayesianMixtures

# Simulate some data
x = randn(500)

# Specify model, data, and MCMC options
n_total = 1000  # total number of MCMC sweeps to run
options = B.options("Normal","MFM",x,n_total, x_full=x_full)  # MFM model with univariate Normal components

# Run MCMC sampler
result = B.run_sampler(options)

# Display the posterior on k (number of components)
println("Posterior on k:")
for (k,pk) in enumerate(B.k_posterior(result; upto=10))
    @printf "%d: %.3g\n" k pk
end


text:
function make_number(data, start, nbytes)
    sum = 0
    for ii = 1:nbytes
        sum *= 256
        sum += data[start+ii]
    end
    return sum
end

function get_xref_record(data, start, w)
    v = Vector{Int}()
    for tw in w
        n = make_number(data, start, tw)
        start += tw
        push!(v, n)
    end
    return v
end

createIfObjectStream(stm::CosStream) =
    get(stm, cn"Type") === cn"ObjStm" ? CosObjectStream(stm) : stm

function read_xref_stream(xrefstm::IDD{CosStream},
                          xref::Dict{CosIndirectObjectRef, CosObjectLoc})
    @assert get(xrefstm, cn"Type") === cn"XRef"
    size = get(xrefstm, cn"Size")
    @assert size !== CosNull

    w = get(xrefstm, cn"W")
    @assert w != CosNull
    @assert length(w) == 3

    index = get(xrefstm, cn"Index")

    index === CosNull &&
        (index = CosArray(CosObject[CosInt(0), size]))

    # The xref stream may be accessed later.
    # There is no point encrypting this data
    # Ideal will be to remove the filter.
    cosStreamRemoveFilters(xrefstm)

    input = get(xrefstm)
    data = read(input)
    util_close(input)
    datasize = length(data)

    w_n = get(w,true) #This size is 3
    recsize = sum(w_n)

    lenidx = length(index)
    @assert rem(lenidx,2) == 0
    idx_int=get(index,true)

    it = 0 #iterator for data
    for i = 1:div(lenidx,2)
        for j = 0:idx_int[2*i] - 1
            oid = idx_int[2*i-1]+j
            record = get_xref_record(data, it, w_n)
            @assert length(record) == 3
            @assert record[1] in 0:2

            loc = (record[1] == 1) ? record[2] :
                  (record[1] == 2) ? record[3] : 0
            stm = (record[1] == 2) ? CosIndirectObjectRef(record[2],0) : CosNull
            ref = (record[1] == 1) ? CosIndirectObjectRef(oid, record[3]) :
                  (record[1] == 2) ? CosIndirectObjectRef(oid, 0) :
                                     CosIndirectObjectRef(0,0)
            it += recsize
            record[1] != 0 && !haskey(xref,ref) && (xref[ref] = CosObjectLoc(loc,stm))
            it >= datasize && break
        end
        it >= datasize && break
    end
    return xref
end

function read_object_info_from_stm(stm::CosStream,
                                   oids::Vector{Int},
                                   oloc::Vector{Int},
                                   n::Int, first::Int)
    filename = get(stm, cn"F")
    io = util_open(String(filename),"r")
    try
        for i = 1:n
            oids[i] = convert(Int, parse_unsignednumber(io))
            oloc[i] = convert(Int, parse_unsignednumber(io)) + first
        end
    finally
        util_close(io)
    end
end

cosObjectStreamGetObject(stm::CosIndirectObject{CosObjectStream},
                         ref::CosIndirectObjectRef, loc::Int) =
                             cosObjectStreamGetObject(stm.obj, ref, loc)

function cosObjectStreamGetObject(stm::CosObjectStream,
                                  ref::CosIndirectObjectRef,
                                  loc::Int)
    objtuple = get(ref)
    stm.oids[loc+1] != objtuple[1] && return CosNull
    dirobj = cosObjectStreamGetObject(stm, CosNull, loc)
    return CosIndirectObject(objtuple[1], objtuple[2], dirobj)
end

function cosObjectStreamGetObject(stm::CosObjectStream,
                                  ref::CosNullType,
                                  loc::Int)
    filename = get(stm, cn"F")
    io = util_open(String(filename),"r")
    ps = io
    try
        seek(ps, stm.oloc[loc+1])
        obj = parse_value(ps)
        return obj
    finally
        util_close(ps)
    end
end


text:
mutable struct Bounds
	min
	max
end

function Bounds(vertices)
	b = Bounds(Vector(0., 0.), Vector(0., 0.))
	update!(b, vertices)
	return b
end

# TODO velocity
function update!(b::Bounds, vertices)
	b.min.x = Inf
	b.max.x = -Inf
	b.min.y = Inf
	b.max.y = -Inf

	for vertex in vertices
		if (vertex.x > b.max.x) b.max.x = vertex.x end
		if (vertex.x < b.min.x) b.min.x = vertex.x end
		if (vertex.y > b.max.y) b.max.y = vertex.y end
		if (vertex.y < b.min.y) b.min.y = vertex.y end
	end
end

function contains!(b::Bounds, point::Vector)
	return (
			point.x >= b.min.x && point.x <= bounds.max.x &&
			point.y >= b.min.y && point.y <= bounds.max.y
			)
end

function overlaps!(boundsA::Bounds, boundsB::Bounds)
	return (boundsA.min.x <= boundsB.max.x && boundsA.max.x >= boundsB.min.x
			&& boundsA.max.y >= boundsB.min.y && boundsA.min.y <= boundsB.max.y)
end

function translate!(b::Bounds, v::Vector)
	b.min += v
	b.max += v
end

function shift!(b::Bounds, pos::Vector)
	delta = b.max - b.min
	b.min = position
	b.max = position + delta
end


text:
using Test
import Pluto: Notebook, ServerSession, ClientSession, Cell, updated_topology, is_just_text

@testset "Analysis" begin
    notebook = Notebook([
        Cell(""),
        Cell("md\"a\""),
        Cell("html\"a\""),
        Cell("md\"a \$b\$\""),
        Cell("md\"a ``b``\""),
        Cell("""
        let
            x = md"a"
            md"r \$x"
        end
        """),
        Cell("html\"a 7 \$b\""),

        Cell("md\"a 8 \$b\""),
        Cell("@a md\"asdf 9\""),
        Cell("x()"),
        Cell("x() = y()"),
        Cell("12 + 12"),
        Cell("import Dates"),
        Cell("import Dates"),
        Cell("while false end"),
        Cell("for i in [16]; end"),
        Cell("[i for i in [17]]"),
        Cell("module x18 end"),
        Cell("""
        module x19
            exit()
        end
        """),
        Cell("""quote end"""),
        Cell("""quote x = 21 end"""),
        Cell("""quote \$(x = 22) end"""),
    ])

    old = notebook.topology
    new = notebook.topology = updated_topology(old, notebook, notebook.cells)

    @testset "Only-text detection" begin
        @test is_just_text(new, notebook.cells[1])
        @test is_just_text(new, notebook.cells[2])
        @test is_just_text(new, notebook.cells[3])
        @test is_just_text(new, notebook.cells[4])
        @test is_just_text(new, notebook.cells[5])
        @test is_just_text(new, notebook.cells[6])
        @test is_just_text(new, notebook.cells[7])

        @test !is_just_text(new, notebook.cells[8])
        @test !is_just_text(new, notebook.cells[9])
        @test !is_just_text(new, notebook.cells[10])
        @test !is_just_text(new, notebook.cells[11])
        @test !is_just_text(new, notebook.cells[12])
        @test !is_just_text(new, notebook.cells[13])
        @test !is_just_text(new, notebook.cells[14])
        @test !is_just_text(new, notebook.cells[15])
        @test !is_just_text(new, notebook.cells[16])
        @test !is_just_text(new, notebook.cells[17])
        @test !is_just_text(new, notebook.cells[18])
        @test !is_just_text(new, notebook.cells[19])
        @test !is_just_text(new, notebook.cells[20])
        @test !is_just_text(new, notebook.cells[21])
        @test !is_just_text(new, notebook.cells[22])
    end
end


text:
using SQLite
using Base.Test, CSV, DataStreams, DataFrames, NullableArrays

import Base: +, ==

if !isdefined(Core, :String)
    typealias String UTF8String
end

a = SQLite.DB()

temp = tempname()
SQLite.DB(temp)

dbfile = joinpath(dirname(@__FILE__),"Chinook_Sqlite.sqlite")
dbfile2 = joinpath(dirname(@__FILE__),"test.sqlite")
# dbfile = joinpath("/Users/jacobquinn/.julia/v0.5/SQLite/test","Chinook_Sqlite.sqlite")
# dbfile2 = joinpath("/Users/jacobquinn/.julia/v0.5/SQLite/test","test.sqlite")
cp(dbfile, dbfile2; remove_destination=true)
db = SQLite.DB(dbfile2)

# DataStreams interface tests
source_table = "employee"
sink_table = "employee2"
sink_file = "employee2.csv"
selall(x) = "select * from " * x

# SQLite.Source ==> DataFrame
ds = SQLite.query(db, selall(source_table))
@test size(ds) == (8,15)
# SQLite.Source ==> SQLite.Sink
SQLite.query(db, selall(source_table), SQLite.Sink, db, sink_table)
ds2 = SQLite.query(db, selall(sink_table))
@test size(ds2) == (8,15)
@test Data.types(ds) == Data.types(ds2) && Data.header(ds) == Data.header(ds2)
# SQLite.Source ==> CSV.Sink
SQLite.query(db, selall(source_table), CSV.Sink, sink_file)
ds4 = CSV.read(sink_file; dateformat="yyyy-mm-dd HH:MM:SS")
@test size(ds4) == (8,15)
@test Data.header(ds) == Data.header(ds4)

SQLite.query(db, selall(source_table), SQLite.Sink, db, sink_table; append=true)
ds5 = SQLite.query(db, selall(sink_table))
@test size(ds5) == (16,15)
@test Data.types(ds) == Data.types(ds2) && Data.header(ds) == Data.header(ds2)

# constructed SQLite.Sink
sink = SQLite.Sink(db, Data.schema(ds5); name=sink_table)
SQLite.query(db, selall(source_table), sink)
ds2 = SQLite.query(db, selall(sink_table))
@test size(ds2) == (8,15)
@test Data.types(ds) == Data.types(ds2) && Data.header(ds) == Data.header(ds2)
SQLite.query(db, selall(source_table), sink; append=true)
ds3 = SQLite.query(db, selall(sink_table))
@test size(ds3) == (16,15)

# constructed SQLite.Source
source = SQLite.Source(db, selall(source_table))
ds = SQLite.query(source)
@test size(ds) == (8,15)
source = SQLite.Source(db, selall(source_table))
SQLite.query(source, SQLite.Sink, db, sink_table)
ds2 = SQLite.query(db, selall(sink_table))
@test size(ds2) == (8,15)
source = SQLite.Source(db, selall(source_table))
SQLite.query(source, SQLite.Sink, db, sink_table; append=true)
ds3 = SQLite.query(db, selall(sink_table))
@test size(ds3) == (16,15)

sink = SQLite.Sink(db, Data.schema(ds); name=sink_table)
source = SQLite.Source(db, selall(source_table))
SQLite.query(source, sink)
ds = SQLite.query(db, selall(sink_table))
@test size(ds) == (8,15)
source = SQLite.Source(db, selall(source_table))
SQLite.query(source, sink; append=true)
ds2 = SQLite.query(db, selall(sink_table))
@test size(ds2) == (16,15)

# SQLite.load
si = SQLite.load(db, sink_table, SQLite.Source, db, selall(source_table))
ds = SQLite.query(db, selall(sink_table))
@test size(ds) == (8,15)
si = SQLite.load(db, sink_table, SQLite.Source, db, selall(source_table); append=true)
ds2 = SQLite.query(db, selall(sink_table))
@test size(ds2) == (16,15)

source = SQLite.Source(db, selall(source_table))
si = SQLite.load(db, sink_table, source)
ds = SQLite.query(db, selall(sink_table))
@test size(ds) == (8,15)
source = SQLite.Source(db, selall(source_table))
si = SQLite.load(db, sink_table, source; append=true)
ds2 = SQLite.query(db, selall(sink_table))
@test size(ds2) == (16,15)

sink = SQLite.Sink(db, Data.schema(ds); name=sink_table)
si = SQLite.load(sink, SQLite.Source, db, selall(source_table))
ds = SQLite.query(db, selall(sink_table))
@test size(ds) == (8,15)
sink = SQLite.Sink(db, Data.schema(ds); name=sink_table, append=true)
si = SQLite.load(sink, SQLite.Source, db, selall(source_table); append=true)
ds2 = SQLite.query(db, selall(sink_table))
@test size(ds2) == (16,15)

source = SQLite.Source(db, selall(source_table))
sink = SQLite.Sink(db, Data.schema(ds); name=sink_table)
si = SQLite.load(sink, source)
ds = SQLite.query(db, selall(sink_table))
@test size(ds) == (8,15)
source = SQLite.Source(db, selall(source_table))
sink = SQLite.Sink(db, Data.schema(ds); name=sink_table, append=true)
si = SQLite.load(sink, source; append=true)
ds2 = SQLite.query(db, selall(sink_table))
@test size(ds2) == (16,15)
SQLite.drop!(db, sink_table)
SQLite.drop!(db, "sqlite_stat1")
rm(sink_file)

# regular SQLite tests
so = SQLite.Source(db,"SELECT name FROM sqlite_master WHERE type='table';")
ds = Data.stream!(so,DataFrame)
@test length(ds.columns) == 1
@test Data.header(ds)[1] == "name"
@test size(ds) == (11,1)

results1 = SQLite.tables(db)
@test Data.types(ds) == Data.types(results1) && Data.header(ds) == Data.header(results1)
@test ds.columns[1].values == results1.columns[1].values

results = SQLite.query(db,"SELECT * FROM Employee;")
@test length(results.columns) == 15
@test size(results) == (8,15)
@test typeof(results[1,1]) == Nullable{Int}
@test typeof(results[1,2]) == Nullable{String}
@test isnull(results[1,5])

SQLite.query(db,"SELECT * FROM Album;")
SQLite.query(db,"SELECT a.*, b.AlbumId
	FROM Artist a
	LEFT OUTER JOIN Album b ON b.ArtistId = a.ArtistId
	ORDER BY name;")

r = SQLite.query(db,"create table temp as select * from album")
@test length(r.columns) == 0
r = SQLite.query(db,"select * from temp limit 10")
@test length(r.columns) == 3
@test size(r) == (10,3)
@test length(SQLite.query(db,"alter table temp add column colyear int").columns) == 0
@test length(SQLite.query(db,"update temp set colyear = 2014").columns) == 0
r = SQLite.query(db,"select * from temp limit 10")
@test length(r.columns) == 4
@test size(r) == (10,4)
@test all(Bool[get(x) == 2014 for x in r[:,4]])
@test length(SQLite.query(db,"alter table temp add column dates blob").columns) == 0
stmt = SQLite.Stmt(db,"update temp set dates = ?")
SQLite.bind!(stmt,1,Date(2014,1,1))
SQLite.execute!(stmt)
r = SQLite.query(db,"select * from temp limit 10")
@test length(r.columns) == 5
@test size(r) == (10,5)
@test typeof(r[1,5]) == Nullable{Date}
@test all(Bool[get(x) == Date(2014,1,1) for x in r[:,5]])
@test length(SQLite.query(db,"drop table temp").columns) == 0

dt = DataFrame(Data.Schema([Float64,Float64,Float64,Float64,Float64],5))
sink = SQLite.Sink(db,dt)
Data.stream!(dt,sink)
Data.close!(sink)
r = SQLite.query(db,"select * from $(sink.tablename)")
@test size(r) == (5,5)
@test Data.header(r) == ["Column1","Column2","Column3","Column4","Column5"]
SQLite.drop!(db,"$(sink.tablename)")

dt = DataFrame(zeros(5,5))
sink = SQLite.Sink(db,dt)
Data.stream!(dt,sink)
Data.close!(sink)
r = SQLite.query(db,"select * from $(sink.tablename)")
@test size(r) == (5,5)
@test all([get(i) for i in r.columns[1]] .== 0.0)
@test all([eltype(i) for i in r.columns[1]] .== Float64)
SQLite.drop!(db,"$(sink.tablename)")

dt = DataFrame(zeros(Int,5,5))
sink = SQLite.Sink(db,dt)
Data.stream!(dt,sink)
Data.close!(sink)
r = SQLite.query(db,"select * from $(sink.tablename)")
@test size(r) == (5,5)
@test all([get(i) for i in r.columns[1]] .== 0)
@test all([eltype(i) for i in r.columns[1]] .== Int)

dt = DataFrame(ones(Int,5,5))
Data.stream!(dt,sink,true) # stream to an existing Sink
Data.close!(sink)
r = SQLite.query(db,"select * from $(sink.tablename)")
@test size(r) == (10,5)
@test [get(i) for i in r.columns[1]] == [0,0,0,0,0,1,1,1,1,1]
@test all([eltype(i) for i in r.columns[1]] .== Int)
SQLite.drop!(db,"$(sink.tablename)")

rng = Date(2013):Date(2013,1,5)
dt = DataFrame([i for i = rng, j = rng])
sink = SQLite.Sink(db,dt)
Data.stream!(dt,sink)
Data.close!(sink)
r = SQLite.query(db,"select * from $(sink.tablename)")
@test size(r) == (5,5)
@test all([get(i) for i in r.columns[1]] .== rng)
@test all([eltype(i) for i in r.columns[1]] .== Date)
SQLite.drop!(db,"$(sink.tablename)")

SQLite.query(db,"CREATE TABLE temp AS SELECT * FROM Album")
r = SQLite.query(db, "SELECT * FROM temp LIMIT ?"; values=[3])
@test size(r) == (3,3)
r = SQLite.query(db, "SELECT * FROM temp WHERE Title LIKE ?"; values=["%time%"])
@test [get(i) for i in r.columns[1]] == [76, 111, 187]
SQLite.query(db, "INSERT INTO temp VALUES (?1, ?3, ?2)"; values=[0,0,"Test Album"])
r = SQLite.query(db, "SELECT * FROM temp WHERE AlbumId = 0")
@test r[1,1] === Nullable(0)
@test get(r[1,2]) == "Test Album"
@test r[1,3] === Nullable(0)
SQLite.drop!(db, "temp")

binddb = SQLite.DB()
SQLite.query(binddb, "CREATE TABLE temp (n NULL, i6 INT, f REAL, s TEXT, a BLOB)")
SQLite.query(binddb, "INSERT INTO temp VALUES (?1, ?2, ?3, ?4, ?5)"; values=Any[SQLite.NULL, convert(Int64,6), 6.4, "some text", b"bytearray"])
r = SQLite.query(binddb, "SELECT * FROM temp")
@test isa(get(r.columns[1][1],SQLite.NULL),SQLite.NullType)
@test isa(get(r.columns[2][1]),Int)
@test isa(get(r.columns[3][1]),Float64)
@test isa(get(r.columns[4][1]),AbstractString)
@test isa(get(r.columns[5][1]),Vector{UInt8})
SQLite.query(binddb, "CREATE TABLE blobtest (a BLOB, b BLOB)")
SQLite.query(binddb, "INSERT INTO blobtest VALUES (?1, ?2)"; values=Any[b"a", b"b"])
SQLite.query(binddb, "INSERT INTO blobtest VALUES (?1, ?2)"; values=Any[b"a", BigInt(2)])
type Point{T}
    x::T
    y::T
end
==(a::Point, b::Point) = a.x == b.x && a.y == b.y
p1 = Point(1, 2)
p2 = Point(1.3, 2.4)
SQLite.query(binddb, "INSERT INTO blobtest VALUES (?1, ?2)"; values=Any[b"a", p1])
SQLite.query(binddb, "INSERT INTO blobtest VALUES (?1, ?2)"; values=Any[b"a", p2])
r = SQLite.query(binddb, "SELECT * FROM blobtest";stricttypes=false)
for v in r.columns[1]
    @test get(v) == b"a"
end
for (v1, v2) in zip(r.columns[2], Any[b"b", BigInt(2), p1, p2])
    @test get(v1) == v2
end
############################################

SQLite.query(db,"CREATE TABLE temp AS SELECT * FROM Album")
r = SQLite.query(db, "SELECT * FROM temp LIMIT :a"; values=Dict(:a => 3))
@test size(r) == (3,3)
r = SQLite.query(db, "SELECT * FROM temp WHERE Title LIKE @word"; values=Dict(:word => "%time%"))
@test [get(i) for i in r.columns[1]] == [76, 111, 187]
SQLite.query(db, "INSERT INTO temp VALUES (@lid, :title, \$rid)"; values=Dict(:rid => 0, :lid => 0, :title => "Test Album"))
r = SQLite.query(db, "SELECT * FROM temp WHERE AlbumId = 0")
@test r[1,1] === Nullable(0)
@test get(r[1,2]) == "Test Album"
@test r[1,3] === Nullable(0)
SQLite.drop!(db, "temp")

r = SQLite.query(db, SQLite.@sr_str("SELECT LastName FROM Employee WHERE BirthDate REGEXP '^\\d{4}-08'"))
@test get(r.columns[1][1]) == "Peacock"

triple(x) = 3x
@test_throws AssertionError SQLite.register(db, triple, nargs=186)
SQLite.register(db, triple, nargs=1)
r = SQLite.query(db, "SELECT triple(Total) FROM Invoice ORDER BY InvoiceId LIMIT 5")
s = SQLite.query(db, "SELECT Total FROM Invoice ORDER BY InvoiceId LIMIT 5")
for (i, j) in zip(r.columns[1], s.columns[1])
    @test_approx_eq get(i) 3*get(j)
end

SQLite.@register db function add4(q)
    q+4
end
r = SQLite.query(db, "SELECT add4(AlbumId) FROM Album")
s = SQLite.query(db, "SELECT AlbumId FROM Album")
@test get(r[1,1]) == get(s[1,1])+4

SQLite.@register db mult(args...) = *(args...)
r = SQLite.query(db, "SELECT Milliseconds, Bytes FROM Track")
s = SQLite.query(db, "SELECT mult(Milliseconds, Bytes) FROM Track")
@test (get(r[1,1]) * get(r[1,2])) == get(s[1,1])
t = SQLite.query(db, "SELECT mult(Milliseconds, Bytes, 3, 4) FROM Track")
@test (get(r[1,1]) * get(r[1,2]) * 3 * 4) == get(t[1,1])

SQLite.@register db sin
u = SQLite.query(db, "select sin(milliseconds) from track limit 5")
@test all(-1 .< convert(Vector{Float64},u[:,1]) .< 1)

SQLite.register(db, hypot; nargs=2, name="hypotenuse")
v = SQLite.query(db, "select hypotenuse(Milliseconds,bytes) from track limit 5")
@test [round(Int,get(i)) for i in v.columns[1]] == [11175621,5521062,3997652,4339106,6301714]

SQLite.@register db str2arr(s) = convert(Array{UInt8}, s)
r = SQLite.query(db, "SELECT str2arr(LastName) FROM Employee LIMIT 2")
@test [get(i) for i in r.columns[1]] == Any[UInt8[0x41,0x64,0x61,0x6d,0x73],UInt8[0x45,0x64,0x77,0x61,0x72,0x64,0x73]]

SQLite.@register db big
r = SQLite.query(db, "SELECT big(5)")
@test get(r[1,1]) == big(5)

doublesum_step(persist, current) = persist + current
doublesum_final(persist) = 2 * persist
SQLite.register(db, 0, doublesum_step, doublesum_final, name="doublesum")
r = SQLite.query(db, "SELECT doublesum(UnitPrice) FROM Track")
s = SQLite.query(db, "SELECT UnitPrice FROM Track")
@test_approx_eq get(r[1,1]) 2*sum(convert(Vector{Float64},s.columns[1]))

mycount(p, c) = p + 1
SQLite.register(db, 0, mycount)
r = SQLite.query(db, "SELECT mycount(TrackId) FROM PlaylistTrack")
s = SQLite.query(db, "SELECT count(TrackId) FROM PlaylistTrack")
@test get(r[1,1]) == get(s[1,1])

bigsum(p, c) = p + big(c)
SQLite.register(db, big(0), bigsum)
r = SQLite.query(db, "SELECT bigsum(TrackId) FROM PlaylistTrack")
s = SQLite.query(db, "SELECT TrackId FROM PlaylistTrack")
@test get(r[1,1]) == big(sum(convert(Vector{Int},s.columns[1])))

SQLite.query(db, "CREATE TABLE points (x INT, y INT, z INT)")
SQLite.query(db, "INSERT INTO points VALUES (?, ?, ?)"; values=[1, 2, 3])
SQLite.query(db, "INSERT INTO points VALUES (?, ?, ?)"; values=[4, 5, 6])
SQLite.query(db, "INSERT INTO points VALUES (?, ?, ?)"; values=[7, 8, 9])
type Point3D{T<:Number}
    x::T
    y::T
    z::T
end
==(a::Point3D, b::Point3D) = a.x == b.x && a.y == b.y && a.z == b.z
+(a::Point3D, b::Point3D) = Point3D(a.x + b.x, a.y + b.y, a.z + b.z)
sumpoint(p::Point3D, x, y, z) = p + Point3D(x, y, z)
SQLite.register(db, Point3D(0, 0, 0), sumpoint)
r = SQLite.query(db, "SELECT sumpoint(x, y, z) FROM points")
@test get(r[1,1]) == Point3D(12, 15, 18)
SQLite.drop!(db, "points")

db2 = SQLite.DB()
SQLite.query(db2, "CREATE TABLE tab1 (r REAL, s INT)")

@test_throws SQLite.SQLiteException SQLite.drop!(db2, "nonexistant")
# should not throw anything
SQLite.drop!(db2, "nonexistant", ifexists=true)
# should drop "tab2"
SQLite.drop!(db2, "tab2", ifexists=true)
@test !in("tab2", SQLite.tables(db2).columns[1])

SQLite.drop!(db, "sqlite_stat1")
@test size(SQLite.tables(db)) == (11,1)

source = SQLite.Source(db,"select * from album")
temp = tempname()
sink = CSV.Sink(temp)
Data.stream!(source,sink)
Data.close!(sink)
dt = Data.stream!(CSV.Source(sink),DataFrame)
@test get(dt[1,1]) == 1
@test get(dt[1,2]) == "For Those About To Rock We Salute You"
@test get(dt[1,3]) == 1

db = nothing; gc(); gc();

db = SQLite.DB()
source = CSV.Source(temp)
sink = SQLite.Sink(db,source; name="temp")
Data.stream!(source,sink)
Data.close!(sink)
source2 = SQLite.Source(sink)
dt = Data.stream!(source2,DataFrame)
@test get(dt[1,1]) == 1
@test string(get(dt[1,2])) == "For Those About To Rock We Salute You"
@test get(dt[1,3]) == 1

sink = SQLite.Sink(db, Data.schema(dt); name="temp2")
Data.stream!(dt,sink)
Data.close!(sink)
source3 = SQLite.Source(sink)
dt = Data.stream!(source3,DataFrame)
@test get(dt[1,1]) == 1
@test string(get(dt[1,2])) == "For Those About To Rock We Salute You"
@test get(dt[1,3]) == 1

#Make sure we handle undefined values
db = SQLite.DB() #In case the order of tests is changed
arr = Array(String,2)
arr[1] = "1" #Now an array with the second value undefined
nv = NullableArrays.NullableArray(arr, [false, true])
schema = DataStreams.Data.Schema(["nv"], [String],2)
d = Any[nv]
dt = DataFrame(d)
SQLite.drop!(db, "temp", ifexists=true)
sink = SQLite.Sink(db, dt; name="temp")
Data.stream!(dt, sink)
Data.close!(sink)
dt2 = SQLite.query(db, "Select * from temp")
#There might be a better way to check this
@test dt.columns[1][1].value==dt2.columns[1][1].value
@test dt.columns[1][2].isnull==dt2.columns[1][2].isnull

#Test removeduplicates!
db = SQLite.DB() #In case the order of tests is changed
ints = Int64[1,1,2,2,3]
strs = String["A", "A", "B", "C", "C"]
nvInts = NullableArrays.NullableArray(ints)
nvStrs = NullableArrays.NullableArray(strs)
schema = Data.Schema(["ints", "strs"], [Int64, String],5)
d = Any[nvInts, nvStrs]
dt = DataFrame(d,[:ints,:strs])
SQLite.drop!(db, "temp", ifexists=true)
sink = SQLite.Sink(db, dt; name="temp")
Data.stream!(dt, sink)
Data.close!(sink)
SQLite.removeduplicates!(db, "temp", ["ints","strs"]) #New format
dt3 = SQLite.query(db, "Select * from temp")
@test get(dt3[1,1]) == 1
@test get(dt3[1,2]) == "A"
@test get(dt3[2,1]) == 2
@test get(dt3[2,2]) == "B"
@test get(dt3[3,1]) == 2
@test get(dt3[3,2]) == "C"

# issue #104
db = SQLite.DB() #In case the order of tests is changed
SQLite.execute!(db, "CREATE TABLE IF NOT EXISTS tbl(a  INTEGER);")
stmt = SQLite.Stmt(db, "INSERT INTO tbl (a) VALUES (@a);")
SQLite.bind!(stmt, "@a", 1)

rm(dbfile2)


text:
# Eric S. Tellez <eric.tellez@infotec.mx>

#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")

import JSON
import Glob

include("textmodel.jl")
include("io.jl")

function harvesine(p1::Vector{Float64}, p2::Vector{Float64})
    # Tweets use GEOJson format: [longitude, latitude]
    lon1, lat1 = p1
    lon2, lat2 = p2
    2 * 6372.8 * asin(sqrt(sind((lat2-lat1)/2)^2 + cosd(lat1) * cosd(lat2) * sind((lon2 - lon1)/2)^2))
end

function radius(fun, input, key, coordinates)
    center = [x |> Float64 for x in coordinates["coordinates"]]
    distance = coordinates["radius"]
    # @show center, distance
    itertweets(input) do tweet
        p = getkeypath(tweet, key)

        if p["type"] == "Point"
            coord = [x |> Float64 for x in p["coordinates"]]
        else  #  p["type"] == "Polygon"
            A = p["coordinates"][1]
            coord = Float64[0.0, 0.0]
            for x in A
                coord[1] += x[1]
                coord[2] += x[2]
            end
            coord[1] /= length(A)
            coord[2] /= length(A)
        end
        if harvesine(center, coord) <= distance
            fun(tweet)
        end
    end
end

function main()
    config = TextConfig()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--query"
        help = "a GeoJSON-formatted dictionary describing the query, e.g. '{\"coordinates\": [long,lat], \"radius\": radius_km}'"
        arg_type = String
        default = "{}"
        "--key"
        help = "specifies the key pointing to the text"
        arg_type = String
        default = "text"
    end

    args = parse_args(ARGS, s)
    radius(STDIN, args["key"], args["query"]) do tweet
        write(STDOUT, JSON.json(tweet), '\n')
    end
end

main()

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
function _cleanup_old_branches(url)
    git_repo_dir = with_cloned_repo(url)

    cd(git_repo_dir) do
        DELETE_OLDER_THAN_HOURS = 3
        branches = String(read(`git branch -r`))  # Get all remote branches
        branches = strip.(split(branches, "\n "))  # Convert string into an array of branch names

        branches = [
            b for b in branches if
            contains(b, "compathelper/new_version/") || contains(b, "integration/")
        ]

        for branch in branches
            # Check if there is a commit within the last 3 hours
            # If there is, is_old will have the commit information
            # If there has not been, it will be empty
            is_old = String(
                read(`git log -1 --since="$(DELETE_OLDER_THAN_HOURS) hours" -s $(branch)`)
            )

            if isempty(is_old)
                branch = replace(branch, "origin/" => "")
                run(`git push -d origin $(branch)`)  # Delete the branch on origin
            end
        end
    end

    return nothing
end

function _generate_branch_name(name::AbstractString)
    return "integration/$(CompatHelper.get_random_string())/$(name)"
end

function templates(parts...)
    test_directory = abspath(joinpath(@__DIR__, ".."))
    templates_directory = joinpath(test_directory, "templates")
    result = joinpath(templates_directory, parts...)

    return result
end

function with_master_branch(
    f::Function,
    path_to_content::AbstractString,
    repo_url::AbstractString,
    parent_branch::AbstractString,
)
    b = generate_branch("master", path_to_content, repo_url, parent_branch)
    return f(b)
end

function generate_branch(
    name::AbstractString,
    path_to_content::AbstractString,
    repo_url::AbstractString,
    parent_branch::AbstractString="master",
)
    b = _generate_branch_name(name)
    git_repo_dir = with_cloned_repo(repo_url)

    cd(git_repo_dir) do
        CompatHelper.git_checkout(parent_branch)
        CompatHelper.git_branch(b; checkout=true)
        empty_git_repo(git_repo_dir)

        for x in readdir(path_to_content)
            src = joinpath(path_to_content, x)
            dst = joinpath(git_repo_dir, x)
            rm(dst; force=true, recursive=true)
            cp(src, dst; force=true)
        end

        CompatHelper.api_retry(() -> CompatHelper.git_add())
        CompatHelper.api_retry(
            () -> CompatHelper.git_commit(
                "Automatic commit - CompatHelper integration tests"
            ),
        )
        CompatHelper.api_retry(() -> CompatHelper.git_push("origin", b))

        rm(git_repo_dir; force=true, recursive=true)
    end

    return b
end

function empty_git_repo(git_repo_dir::AbstractString)
    cd(git_repo_dir) do
        for x in readdir(git_repo_dir)
            if x != ".git"
                path = joinpath(git_repo_dir, x)
                rm(path; force=true, recursive=true)
            end
        end
    end

    return nothing
end

function with_cloned_repo(repo_url)
    git_repo_dir = joinpath(mktempdir(), CompatHelper.LOCAL_REPO_NAME)
    CompatHelper.api_retry(() -> CompatHelper.git_clone(repo_url, git_repo_dir))

    return git_repo_dir
end


text:
using PedestrianAvoidancePOMDP
using Test

@testset "PedestrianAvoidancePOMDP.jl" begin
    # Write your own tests here.
    @test 1 == 2
end


text:
module myexample
using ForwardDiff
include("func.jl")
export my_f,my_derivative_f

end


text:
using Test
using SKTB, JuLIP, LinearAlgebra
using JuLIP.Testing

COMPAREQUIP = false     # if QUIP and quippy are installed
TESTDEPTH = 1

println("============================================")
println("    SKTB Tests  ")
println("============================================")

# =========== Main tests =================
@testset "SKTB Tests" begin
   include("testtoymodel.jl")
   include("testnrltb.jl")
   include("testcontour.jl")
   include("testsiteE.jl")
   include("testkwon.jl")
   include("testdual.jl")
   include("test0T.jl")
   include("testvirial.jl")
end

# ============= Compare with Atoms.jl and QUIP implementations
if COMPAREQUIP; include("comparequip.jl"); end

# ============= Performance benchmarks
# (uncomment these only for performance tests)
# include("benchmarkEandFrc.jl")
# include("perfsiteE.jl")

# =========== TEMP =======


text:
#=
A builder is looking to build a row of N houses that can be of K different colors. He has a goal of minimizing cost while ensuring that no two neighboring houses are of the same color.

Given an N by K matrix where the nth row and kth column represents the cost to build the nth house with kth color, return the minimum cost which achieves this goal.
=#

function minimum_building_cost(cost_matrix)
    # Create a matrix with the cumulative minimum cost of building the nth house (row) with the kth color (column)
    cumulative_minimum_cost = zeros(Float64, size(cost_matrix))
    # For the first house, it is the first row of the cost_matrix
    cumulative_minimum_cost[1,:] = cost_matrix[1,:]
    number_of_colors = size(cost_matrix, 2)
    number_of_houses = size(cost_matrix, 1)
    # For the second to last houst, for each color, the minimum cumulative cost is the price of building the house in that color + the minimum cumulative cost of building the last house in a different color
    for house in 2:number_of_houses
        for color in 1:number_of_colors
            allowed_prev_colors = setdiff(1:number_of_colors, color)
            cumulative_minimum_cost[house, color] = cost_matrix[house, color] + minimum(cumulative_minimum_cost[house-1,allowed_prev_colors])
        end
    end
    # We want the minimum cost of building all the house
    minimum_cost = minimum(cumulative_minimum_cost[end,:])
    return minimum_cost
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using SchumakerSpline
using Test
tol = 100 * eps()
close(a,b) = abs(a-b) < tol

from = 0.5
to = 10
x1 = collect(range(from, stop=to, length=40))
y1 = (x1).^2
x2 = [0.5,0.75,0.8,0.93,0.9755,1.0,1.1,1.4,2.0]
y2 = sqrt.(x2)
s1 = Schumaker(x1,y1)
s2 = Schumaker(x2,y2)
con = 7.6
x = 0.9

close( (s1 +con)(x) , s1(x) + con)
close( (con + s1)(x) , s1(x) + con)

close( (s1 - con)(x) , s1(x) - con)
close( (con - s1)(x) , con - s1(x) )

close( (s1 *con)(x) , s1(x) * con)
close( (con * s1)(x) , s1(x) * con)

close( (s1 /con)(x) , s1(x) / con)


text:
using Images, ImageIO, FileIO
using ImageInTerminal, ImageShow
using ColorVectorSpace
using ImageFiltering, ImageContrastAdjustment

function diffusion!(img_new, img, N)
    m, n = size(img)
    for i in 1:m
        for j in 1:n
            if img[i, j] == 0
                continue
            else
                img_new[i, j] = img[i, j]
                for k in 1:N
                    i_new = i + rand(-1:1)
                    j_new = j + rand(-1:1)
                    if i_new < 1 || i_new >= m-1 || j_new < 1 || j_new >= n-1
                        continue
                    else
                        img_new[i_new, j_new] = img[i, j]
                    end
                end
            end
        end
    end
end

function diffusions!(imgs, img, N)
    m, n = size(img)
    for i in 1:m
        for j in 1:n
            if img[i, j] == 0
                continue
            else
                for k in 1:N
                    img_new = imgs[k]
                    img_new[i, j] = img[i, j]
                    i_new = i + rand(-1:1)
                    j_new = j + rand(-1:1)
                    if i_new < 1 || i_new >= m-1 || j_new < 1 || j_new >= n-1
                        continue
                    else
                        img_new[i_new, j_new] = img[i, j]
                    end
                end
            end
        end
    end
end

function diffusion(img, N)
    img_new = similar(img)
    diffusion!(img_new, img, N)
    img_new
end

function animate_diffusion(img, N)
    imgs = [similar(img) for _ in 1:N]
    diffusions!(imgs, img, N)
    imgs
end

function animate_filt(img, N)
    imgs = [similar(img) for _ in 1:N]
    for i in 1:N
        # diffusions!(imgs[i], img, N)
        imfilter!(imgs[i], img, Kernel.gaussian(i));
        # imfilter!(imgs[i], img, Kernel.sobel());
    end
    imgs
end

function anim(imgs; sleep_amt=0.1)
    for i in 1:length(imgs)
        display(imgs[i])
        sleep(sleep_amt)
    end
end

img = load("/Users/anand/Pictures/81.jpeg")
img = load("/Users/anand/Pictures/81.jpeg")
img = load("pic.jpg")
img = ColorVectorSpace.complement.(img)

alg = Equalization(nbins = 256)
img = adjust_histogram(img, alg)

imgs = animate_diffusion(img, 1000)
imgs = animate_filt(img, 100)
anim(imgs)

for (i, img) in enumerate(imgs)
    save("outs/$i.jpeg", img)
end

run(`ffmpeg -framerate 60 -start_number 1 -i 'out2/%d.jpeg' -r 60 -y vid/test2.mp4`)

imgs[1] imgs[end] == imgs[1]

imgg = imfilter(img_adjusted, Kernel.gaussian(3));
imgg = imfilter(imgg, Kernel.gaussian(100))

max_stars_repo_name: anandijain/audio

max_issues_repo_name: anandijain/audio

max_forks_repo_name: anandijain/audio


text:
using Interpolations
using Tullio

##
x = 0:0.1:1 |> collect
f = exp.(x)

interp_lin = LinearInterpolation(x, f)

x_new = 0:(1//3):1 |> collect
f_new = interp_lin(x_new)

##
y = copy(x)
g = zeros(size(x, 1), size(y, 1))
@tullio g[i, j] = exp(x[i])

interp_lin_2d = LinearInterpolation((x, y), g)

y_new = copy(x_new)
g_new = interp_lin_2d(x_new, y_new)


text:
using PiecewiseLinearApprox
using Test

@testset "PiecewiseLinearApprox.jl" begin

  function foo(x)

    y = log.(x[1])*log.(x[2])'
    return y

  end

  x1 = piecewise_linear_nodes(11,[3.0,1.5])
  x2 = piecewise_linear_nodes(21,[2.5,0.5])

  x = (x1,x2)
  y = foo(x)

  point = [2.1,0.7]

  y_hat = piecewise_linear_evaluate(y,x,point)

  y_actual = log(point[1])*log(point[2])
  diff = (y_actual - y_hat)
  println(diff)
  @test abs(diff) <= 1e-3

end


text:
module WordNet

using DataDeps

include("init.jl")
include("pointer.jl")
include("lemma.jl")
include("constants.jl")
include("synset.jl")
include("db.jl")
include("operations.jl")
include("sensekeys.jl")

end


text:
function prime_factors(n)
    pf = []
    temp = n
    i = 2
    while i <= temp
        if temp % i == 0
            push!(pf, i)
            temp ÷= i
        else
            i += 1
        end
    end
    pf
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
using ITensors, Test

@testset "README Examples" begin
  @testset "ITensor Basics" begin
    i = Index(3)
    j = Index(5)
    k = Index(2)
    l = Index(7)

    A = ITensor(i, j, k)
    B = ITensor(j, l)

    A[i(1), j(1), k(1)] = 11.1
    A[i(2), j(1), k(2)] = -21.2
    A[k(1), i(3), j(1)] = 31.1  # can provide Index values in any order
    # ...

    # Contract over shared index j
    C = A * B

    @test hasinds(C, i, k, l) == true

    D = randomITensor(k, j, i) # ITensor with random elements

    # Add two ITensors
    # must have same set of indices
    # but can be in any order
    R = A + D
  end

  @testset "SVD of a Matrix" begin
    i = Index(10)
    j = Index(20)
    M = randomITensor(i, j)
    U, S, V = svd(M, i)
    @test norm(M - U * S * V) < 1E-12
  end

  @testset "SVD of a Tensor" begin
    i = Index(4, "i")
    j = Index(4, "j")
    k = Index(4, "k")
    l = Index(4, "l")
    T = randomITensor(i, j, k, l)
    U, S, V = svd(T, i, k)
    @test hasinds(U, i, k)
    @test hasinds(V, j, l)
    @test norm(T - U * S * V) < 1E-12
  end

  @testset "Making Tensor Indices" begin
    i = Index(3)     # Index of dimension 3
    @test dim(i) == 3     # dim(i) = 3

    ci = copy(i)
    @test ci == i    # true

    j = Index(5, "j") # Index with a tag "j"

    @test j != i     # false

    s = Index(2, "n=1,Site") # Index with two tags,
    # "Site" and "n=1"
    @test hastags(s, "Site") # hastags(s,"Site") = true
    @test hastags(s, "n=1")  # hastags(s,"n=1") = true

    i1 = prime(i) # i1 has a "prime level" of 1
    # but otherwise same properties as i
    @test i1 != i # false, prime levels do not match
  end

  @testset "DMRG" begin
    N = 100
    sites = siteinds("S=1", N)

    # Input operator terms which define
    # a Hamiltonian matrix, and convert
    # these terms to an MPO tensor network
    ampo = AutoMPO()
    for j in 1:(N - 1)
      add!(ampo, "Sz", j, "Sz", j + 1)
      add!(ampo, 0.5, "S+", j, "S-", j + 1)
      add!(ampo, 0.5, "S-", j, "S+", j + 1)
    end
    H = MPO(ampo, sites)

    # Create an initial random matrix product state
    psi0 = randomMPS(sites)

    sweeps = Sweeps(2)
    maxdim!(sweeps, 10, 20, 100, 100, 200)
    cutoff!(sweeps, 1E-10)

    # Run the DMRG algorithm, returning energy
    # (dominant eigenvalue) and optimized MPS
    energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)
    #println("Final energy = $energy")
  end
end

nothing