text:
include("../util.jl")

using JuMP, Gurobi, MathOptFormat

function TEAVAR_Star(env,
                edges,
                capacity,
                flows,
                demand,
                beta,
                k,
                T,
                Tf,
                scenarios,
                scenario_probs,
		outputfilename;
                max_concurrent_flow=true,
                explain=false,
                verbose=false,
                utilization=false,
                average=false)

    nedges = length(edges)
    nflows = length(flows)
    ntunnels = length(T)
    nscenarios = length(scenarios)
    p = scenario_probs

    println(Dates.format(now(), "HH:MM:SS"), ": #edges ", nedges, " #flows ", nflows, " #tunnels ", ntunnels, " #scenarios ", nscenarios)

    #CREATE TUNNEL SCENARIO MATRIX
    X  = ones(nscenarios,ntunnels)
    for s in 1:nscenarios
        for t in 1:ntunnels
            if size(T[t],1) == 0
                X[s,t] = 0
            else
                for e in 1:nedges
                    if scenarios[s][e] == 0
                        back_edge = findfirst(x -> x == (edges[e][2],edges[e][1]), edges)
                        if in(e, T[t]) || in(back_edge, T[t])
                        # if in(e, T[t])
                            X[s,t] = 0
                        end
                    end
                end
            end
        end
    end

    println(Dates.format(now(), "HH:MM:SS"), ": created tunnel scenario matrix")

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    println(Dates.format(now(), "HH:MM:SS"), ": created tunnel edge matrix")

    model = Model(solver=GurobiSolver(env, OutputFlag=1, Method=2, Crossover=0))
    # flow per commodity per path variables
    # @variable(model, a[1:nflows, 1:k] >= 0, basename="a", category=:SemiCont)
	@variable(model, a[1:nflows, 1:k] >= 0, basename="a")

    # alpha variable
    # @variable(model, alpha >= 0, basename="alpha", category=:SemiCont)
    @variable(model, alpha >= 0, basename="alpha")

    # maximum flow lost in that scenario
    @variable(model, umax[1:nscenarios] >= 0, basename="umax")
    # flow lost per commod per scenario
    @variable(model, u[1:nscenarios, 1:nflows] >= 0, basename="u")

    # capacity constraints for final flow assigned to "a" variables
    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows, t in 1:size(Tf[f],1)) <= capacity[e])
    end

    println(Dates.format(now(), "HH:MM:SS"), ": some variables and constraints are created")
    total_scenario_prob = sum(p[s] for s in 1:nscenarios)
    println("total_scenario_prob= ", total_scenario_prob)
    total_demand = sum(demand[f] for f in 1:nflows)

    if max_concurrent_flow
      println("Running max concurrent flow")
      # FLOW LEVEL LOSS
      @expression(model, satisfied[s=1:nscenarios, f=1:nflows], sum(a[f,t] * X[s,Tf[f][t]] for t in 1:size(Tf[f],1)) / demand[f])

      for s in 1:nscenarios
          for f in 1:nflows
              @constraint(model, u[s,f] >= 1 - satisfied[s,f])
          end
      end

      for s in 1:nscenarios
          if average
              @constraint(model, umax[s] + alpha >= (sum(u[s,f] for f in 1:nflows)) / nflows)
          else
              for f in 1:nflows
                  @constraint(model, umax[s] + alpha >= u[s,f])
              end
          end
      end

      @objective(model, Min, alpha + (1 / (1 - beta)) * (sum((p[s] * umax[s] for s in 1:nscenarios)) + ((1-total_scenario_prob) * (1-alpha) )))

    else
      # TODO: do these loss formulas hold if demands are not satisfiable?
      println("Running max total flow")
      @expression(model, satisfied[s=1:nscenarios, f=1:nflows], sum(a[f,t] * X[s, Tf[f][t]] for t in 1:size(Tf[f],1)))

      for s in 1:nscenarios
	@constraint(model, umax[s] >= 0)
	for f in 1:nflows
	  # note: satisfied[s,f] can be more than demand[f]
	  @constraint(model, u[s, f] >= demand[f] - satisfied[s,f])
	end
	# @expression(model, sum_u_f[s=1:nscenarios], sum(u[s,f] for f in 1:nflows))
	@constraint(model, umax[s] + alpha >= sum(u[s, f] for f in 1:nflows))
      end

      @objective(model, Min, alpha + (1/ (1-beta)) * (sum((p[s] * umax[s] for s in 1:nscenarios)) + ((1-total_scenario_prob)*(total_demand-alpha))))
    end

    # println(Dates.format(now(), "HH:MM:SS"), ": writing to file")
    # writeMPS(model, "blahblahjl.mps")
    # MOI.write_to_file(model, "blahblahjl.lp")

    println(Dates.format(now(), "HH:MM:SS"), ": ready to solve")
    flush(stdout)

    solve(model)

    println(Dates.format(now(), "HH:MM:SS"), ": solver finished; explaining")
    if (explain)
        println("Runtime: ", getsolvetime(model))

	# some simple debugging information
	println("beta: ", beta)
	println("#flows: ", nflows)
	println("#edges: ", nedges)
	println("#tunnels: ", ntunnels)
	println("#demands: ", nflows, " total demand=", total_demand)
	println("#scenarios: ", nscenarios, " total_prob= ", sum(p[s] for s in 1:nscenarios))

	result_u = getvalue(u)
	result_umax = getvalue(umax)
	result_satisfied = getvalue(satisfied)
	result_alpha = getvalue(alpha)

	for f in 1:nflows
		sat_vector=Array{Float64}(undef, nscenarios)
		for s in 1:nscenarios
			sat_vector[s] = result_satisfied[s,f]
		end
		println("satisfied[*,", f, "]: ", sat_vector)
	end

	per_scenario_losses = Array{Tuple{Float64,Int64,Float64}}(undef, nscenarios)

	if max_concurrent_flow
		loss_s_f = ones(nflows)' .- result_satisfied
	else
		loss_s_f = demand' .- result_satisfied
	end
	loss_s_f[loss_s_f .< 0] .= 0

	for s in 1:nscenarios
	     if max_concurrent_flow
		loss_s = maximum(loss_s_f[s, :])
	     else
		loss_s = sum(loss_s_f[s, :])
	     end
	     # max_loss_s = maximum(loss_s_f[s, :])
	     # max_u_s = maximum(result_u[s, :])
	     # println(s, ": max_loss_s/max_u_s/result_umax_s= ", max_loss_s, " ", max_u_s, " ", result_umax[s])
	     per_scenario_losses[s] = (loss_s, s, p[s])
        end
	println("psl before sort=", per_scenario_losses)
	per_scenario_losses = sort(per_scenario_losses, by = first)
	println("psl after sort=", per_scenario_losses)

	println("Alpha after obj=", result_alpha)

	result_a = getvalue(a)

	total_prob_thus_far = 0
	ell_k_beta = -1
	k_beta = -1
	for ss in 1:nscenarios
		tup = per_scenario_losses[ss]
		total_prob_thus_far += p[tup[2]]
		if total_prob_thus_far > beta
			ell_k_beta = tup[1]
			k_beta = ss
			break;
		end
	end
	if total_prob_thus_far <= beta
		ell_k_beta = result_alpha
	end
	println("Sorting scenarios by loss: beta= ", beta, " k_beta= ", k_beta, " ell_k_beta= ", ell_k_beta)

	if max_concurrent_flow
	else
		# x
		result_alpha /= total_demand
		ell_k_beta /= total_demand
	end

	# compute net flow
	numtunnel_perflow = 0
	for f in 1:nflows
		numtunnel_perflow = max(numtunnel_perflow, size(Tf[f], 1))
	end

	flow_allocation_1 = Array{Float64}(undef, nflows)
	flow_allocation_2 = Array{Float64}(undef, nflows)
	flow_tunnel_weights = Matrix{Float64}(undef, nflows, numtunnel_perflow)
	all_flow_alloc = 0

	for f in 1:nflows
		flow_totalalloc = 0
	 	for t in 1:size(Tf[f], 1)
	  		flow_totalalloc += result_a[f,t]
         	end

	 	all_flow_alloc += flow_totalalloc

	 	flow_allocation_1[f] = (1 - result_alpha) * min(flow_totalalloc, demand[f])
	 	flow_allocation_2[f] = (1 - ell_k_beta) * min(flow_totalalloc, demand[f])
	 	for t in 1:size(Tf[f], 1)
	   		weight = 0
	   		if flow_totalalloc > 0
				weight = result_a[f,t]/ flow_totalalloc
	   		end
	   		flow_tunnel_weights[f,t] = weight
        end
	end

	flow_allocation_3 = Array{Float64}(undef, nflows)
	for f in 1:nflows
		flow_per_scenario_losses = Array{Tuple{Float64,Int64,Float64}}(undef, nscenarios)
		for ss in 1:nscenarios
			flow_per_scenario_losses[ss] = (loss_s_f[ss, f], ss, p[ss])
		end

		flow_per_scenario_losses = sort(flow_per_scenario_losses, by = first)

		flow_total_prob_thus_far = 0
		flow_loss = 0
		flow_loss_cutoff_iter = 0
		if flow_total_prob_thus_far < beta
			for ss in 1:nscenarios
				tup = flow_per_scenario_losses[ss]
				flow_total_prob_thus_far += tup[3]
				if flow_total_prob_thus_far >= beta
					flow_loss = tup[1]
					flow_loss_cutoff_iter = ss
					break
				end
			end
			if flow_loss_cutoff_iter == 0
				println("WARN!!!! did not get flow loss")
			else
				println("Flow ", f, " loss = ", flow_loss, " at cutoffindex= ", flow_loss_cutoff_iter, " totprob= ", flow_total_prob_thus_far)
				println("Flow losses sorted: ", flow_per_scenario_losses)
			end
		end

		if max_concurrent_flow
			flow_allocation_3[f] = demand[f] * (1- flow_loss)
		else
			flow_allocation_3[f] = demand[f] - flow_loss
		end
	end

	# to account for minor off-by-epsilon negative values
	# should be asserting
	flow_allocation_1[flow_allocation_1 .< 0] .= 0
	flow_allocation_2[flow_allocation_2 .< 0] .= 0
	flow_allocation_3[flow_allocation_3 .< 0] .= 0

	println("totalalloc= ", all_flow_alloc, " netalloc1/2/3= ",
		sum(flow_allocation_1[f] for f in 1:nflows), " ",
		sum(flow_allocation_2[f] for f in 1:nflows), " ",
		sum(flow_allocation_3[f] for f in 1:nflows))

	output=open(outputfilename, "w")
	write(output, "---flow_allocation with alpha---\n")
	writedlm(output, flow_allocation_1, ',')
	write(output, "---flow_allocation with ell_k_beta---\n")
	writedlm(output, flow_allocation_2, ',')
	write(output, "---flow_allocation with per_flow_loss---\n")
	writedlm(output, flow_allocation_3, ',')
	write(output, "---weights per flow per tunnel---\n")
	writedlm(output, flow_tunnel_weights, ',')
	flush(output)
	close(output)

        printResults(getobjectivevalue(model), result_alpha, result_a, result_u, result_umax, edges, scenarios, T, Tf, L, capacity, p, demand, verbose=verbose, utilization=utilization)
    end

    return (getvalue(alpha), getobjectivevalue(model), getvalue(a), getvalue(umax), getsolvetime(model))
end


text:
"""If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.
Find the sum of all the multiples of 3 or 5 below 1000."""

function sum_multiples_3_5(n::Int)
    multiples3or5 = union(3:3:n-1 , 5:5:n-1)
    return sum(multiples3or5)
end

@assert sum_multiples_3_5(10) ==23

@show sum_multiples_3_5(1000)

# Oneliner and faster
@show sum([i for i in 1:999 if i%3==0 || i%5 ==0]);

ext: jl


text:
# CUDA CuDevice management

function devcount()
    # Get the number of CUDA-capable CuDevices
    a = Ref{Cint}()
    @cucall(cuDeviceGetCount, (Ptr{Cint},), a)
    return Int(a[])
end

immutable CuDevice
    ordinal::Cint
    handle::Cint

    function CuDevice(i::Int)
        ordinal = convert(Cint, i)
        a = Ref{Cint}()
        @cucall(cuDeviceGet, (Ptr{Cint}, Cint), a, ordinal)
        handle = a[]
        new(ordinal, handle)
    end
end

immutable CuCapability
    major::Int
    minor::Int
end

function name(dev::CuDevice)
    const buflen = 256
    buf = Array(Cchar, buflen)
    @cucall(cuDeviceGetName, (Ptr{Cchar}, Cint, Cint), buf, buflen, dev.handle)
    bytestring(pointer(buf))
end

function totalmem(dev::CuDevice)
    a = Ref{Csize_t}()
    @cucall(cuDeviceTotalMem, (Ptr{Csize_t}, Cint), a, dev.handle)
    return Int(a[])
end

function attribute(dev::CuDevice, attrcode::Integer)
    a = Ref{Cint}()
    @cucall(cuDeviceGetAttribute, (Ptr{Cint}, Cint, Cint), a, attrcode, dev.handle)
    return Int(a[])
end

capability(dev::CuDevice) = CuCapability(attribute(dev, 75), attribute(dev, 76))

function list_devices()
    cnt = devcount()
    if cnt == 0
        println("No CUDA-capable CuDevice found.")
        return
    end

    for i = 0:cnt-1
        dev = CuDevice(i)
        nam = name(dev)
        tmem = iround(totalmem(dev) / (1024^2))
        cap = capability(dev)

        println("device[$i]: $(nam), capability $(cap.major).$(cap.minor), total mem = $tmem MB")
    end
end


text:
# Methods for fitting and querying under the RKHS regularization framework.

# for univariate inputs.
function constructkernelmatrix( X,
                                θ)::Matrix{Float64} where T
    #
    K = Matrix{Float64}(undef,length(X),length(X))
    constructkernelmatrix!(K,X,θ)

    return K
end

function constructkernelmatrix!(K::Matrix{T},
                                X,
                                θ)::Matrix{T} where T

    M = length(X)
    @assert size(K) == (M,M)

    fill!(K,Inf) # debug
    for j = 1:M
        for i = j:M
            K[i,j] = evalkernel(X[i], X[j], θ)
        end
    end

    for j = 2:M
        for i = 1:(j-1)
            K[i,j] = K[j,i]
        end
    end

    return K
end

# experimental.
function constructkernelmatrix2(X::Vector{Vector{T}}, θ_array::Vector)::Matrix{T} where T
    M = length(X)
    @assert length(θ_array) == M

    K = Matrix{T}(undef,M,M)
    fill!(K,Inf) # debug
    for j = 1:M
        #for i = j:M
        for i = 1:M
            K[i,j] = evalkernel(X[i], X[j], θ_array[j])
        end
    end

    # for j = 2:M
    #     for i = 1:(j-1)
    #         K[i,j] = K[j,i]
    #     end
    # end
    # K = (K+K') ./ 2
    #
    # if !isposdef(K)
    #     s,Q = eigen(K)
    #     println("K not posdef! eigenvalues are:")
    #     display(sort(s))
    # end

    return K
end

# v are the derivatives of x.
function constructkernelematrices( X::Vector{T},
                                   θ::ODEKernelType) where T
    N_x = length(X)
    N_v = length(X)

    K = constructkernelmatrix(X, θ)

    K_xv = Matrix{T}(undef, N_x, N_v)
    fill!(K_xv,Inf) # debug
    for j = 1:N_v
        for i = 1:N_x
            K_xv[i,j] = evalkernelderivative(X[i], X[j], θ)
        end
    end

    K_vv = Matrix{T}(undef, N_x, N_v)
    fill!(K_vv,Inf) # debug
    for j = 1:N_v
        for i = 1:N_x
            K_vv[i,j] = evalDkernel𝐷(X[i], X[j], θ)
        end
    end

    return K, K_xv, K_vv
end

# returns K_XZ. Entry i,j is k(x[i],z[j]).
function constructkernelmatrix( X::Vector{Vector{T}},
                                Z::Vector{Vector{T}},
                                θ)::Matrix{T} where T
    Nr = length(X)
    Nc = length(Z)

    K = Matrix{T}(undef,Nr,Nc)
    fill!(K,Inf) # debug
    for j = 1:Nc
        for i = 1:Nr
            K[i,j] = evalkernel(X[i], Z[j], θ)
        end
    end

    return K
end

function constructkernelmatrix!(K,
            η::RKHSProblemType{KT,T,XT}) where {KT,T,XT}

    #
    constructkernelmatrix!(K, η.X, η.θ)

    return nothing
end

function constructkernelmatrix!(K,
            η::RKHSProblemType{FastAdaptiveKernelType{KT,T,D},T,XT}) where {KT,T,XT,D}

    constructkernelmatrix!(K, η.X, η.θ)

    return nothing
end

# for multivariate inputs.
# for use with adaptive kernels.
function constructkernelmatrix!(    K::Matrix{T},
                                    X::Vector{Vector{T}},
                                    θ::FastAdaptiveKernelType{KT,T,D}) where {KT,T,D}
    M = length(X)
    @assert size(K) == (M,M)

    # update warpmap evals.
    for i = 1:length(θ.warpfuncs)
        for n = 1:length(X)
            θ.w_X[n,i] = θ.warpfuncs[i](X[n])
        end
    end

    #sb = one(T) - sum(θ.s)
    #@assert one(T) >= sb > zero(T)
    sb = one(T)

    fill!(K,Inf) # debug
    for j = 1:M
        for i = j:M
            if i == j
                K[i,j] = one(T)
            else
                K[i,j] = evalkernel(X[i], X[j], θ)
            end
        end
    end

    for j = 2:M
        for i = 1:(j-1)
            K[i,j] = K[j,i]
        end
    end

    return nothing
end

function constructkernelmatrix!( X::Vector{Vector{T}},
                                θ::FastAdaptiveKernelType)::Matrix{T} where T
    #
    K = Matrix{T}(undef,length(X),length(X))
    constructkernelmatrix!(K, X, θ)

    return K
end

#### Try separable RKHS via inverse. Assume N_r == D
# Assume η.a is same for all dimensions.

function fitRKHS!(  η,
                    y::Vector{T}) where T
    #
    N = length(η.X)
    K = Matrix{T}(undef, N, N)
    fitRKHS!(K, η, y)

    return nothing
end

"""
U is K + σ²I.
"""
function fitRKHS!(  U::Matrix{T},
                    η,
                    y::Vector{T}) where T

    @assert !isempty(η.X)
    @assert !isempty(y)

    M = length(η.X)
    @assert M == length(y)

    # update kernel matrix
    constructkernelmatrix!(U, η.X, η.θ)

    # add observation model's noise.
    for i = 1:size(U,1)
        U[i,i] += η.σ²[1]
    end

    # solve.
    η.c[:] = U\y

    return nothing
end

# Version where Xq being SharedArray, ! versoin.
function query!(    Yq::Vector{T},
                    Xq::Vector{Vector{T}},
                    η::RKHSProblemType{KT,T}) where {KT, T}

    N = length(η.X)
    @assert !isempty(Xq)

    @assert size(Yq) == size(Xq)

    # Pre-allocate
    kq = Vector{T}(undef,N)
    Nq = length(Xq)

    for iq = 1:Nq

        for i=1:N
            kq[i] = evalkernel(Xq[iq], η.X[i], η.θ)
        end

        Yq[iq] = dot(kq,η.c)

        # # variance.
        # v = L\kq
        # vq[iq] = evalkernel(Xq[iq], Xq[iq], η.θ) - dot(v,v)
    end

    return nothing
end

function query!(    Yq::Vector{T},
                    Xq::Vector{T},
                    η::RKHSProblemType{KT,T}) where {KT, T}

    N = length(η.X)
    @assert !isempty(Xq)

    @assert size(Yq) == size(Xq)

    # Pre-allocate
    kq = Vector{T}(undef,N)
    Nq = length(Xq)

    for iq = 1:Nq

        for i=1:N
            kq[i] = evalkernel(Xq[iq], η.X[i], η.θ)
        end

        Yq[iq] = dot(kq,η.c)

        # # variance.
        # v = L\kq
        # vq[iq] = evalkernel(Xq[iq], Xq[iq], η.θ) - dot(v,v)
    end

    return nothing
end

function query!(    Yq::Vector{T},
                    Xq::Vector{Vector{T}},
                    η::RKHSProblemType{Vector{KT},T}) where {KT, T}

    N = length(η.X)
    @assert !isempty(Xq)

    @assert size(Yq) == size(Xq)

    # Pre-allocate
    kq = Vector{T}(undef,N)
    Nq = length(Xq)

    for iq = 1:Nq

        for i=1:N
            kq[i] = evalkernel(Xq[iq], η.X[i], η.θ[i])
        end

        Yq[iq] = dot(kq,η.c)

        # # variance.
        # v = L\kq
        # vq[iq] = evalkernel(Xq[iq], Xq[iq], η.θ) - dot(v,v)
    end

    return nothing
end

#### posdef checks.
function zerobarrier(x::T, a::T)::T where T <: Real

    #return one(T)/(a*x)^2
    return exp(-(x-a))
    #return abs(x)*a
end

function zerobarrier(x::Complex{T}, a::T)::T where T <: Real

    #return one(T)/abs2(a*x)
    return exp(-(abs(x)-a))
    #return abs(x)*a
end

function zerobarrier(x::Vector{T2}, a::T)::T where {T <: Real, T2}

    #return maximum( zerobarrier(x[d], a) for d = 1:length(x) )
    return sum( zerobarrier(x[d], a) for d = 1:length(x) )/length(x)
end

# function piecewisezerobarrier(x::T, a::T, b::T)::T where T <: Real
#     @assert b > a >= 0.0
#
#     c = (b-1)/a
#
#     if a < x < b
#         return -log(b-c*x)
#     end
#
#     if x > b
#         return Inf
#     end
#
#     return zero(T)
# end

# x must be positive, a is the threshold for 0.0.
function piecewisezerobarrier(x::T, a::T)::T where T <: Real
    @assert x > zero(T)

    if a < x
        return exp(x-a) - one(T)
    end

    return zero(T)
end

function logisticfunc(L, k, x0::T, x)::T where T <: Real

    return L/(1+exp(-k*(x-x0)))
end

"""
y is data here.
"""
function discrepancyscorefromdata(  K::Matrix{T},
                                     y::Vector{T},
                                     σ²::T,
                                     c::Vector{T};
                                     a_cond::T = 1e-10) where T

   discrepancy = norm(( K + σ² .* LinearAlgebra.I )*c - y)

   return exp(discrepancy*a_cond)-one(T)
end

function discrepancyscorefromdata(  U::Matrix{T},
                                     y::Vector{T},
                                     c::Vector{T};
                                     a_cond::T = 1e-10) where T

   discrepancy = norm(U*c - y)

   return exp(discrepancy*a_cond)-one(T)
end

"""
The lower the score, the more well-conditioned K is.
"""
function evalconditioningscore( K::Matrix{T},
                                display_flag::Bool,
                                a_posdef_chk::T,
                                y::Vector{T},
                                σ²::T,
                                c::Vector{T})::T where T

    #s, Q_unused = eigen(K)
    #score = zerobarrier(s, a_posdef_chk)

    # score = piecewisezerobarrier(   cond(K),
    #                                 a_posdef_chk/eps(T))
    score = discrepancyscorefromdata(K, y, σ², c; a_cond = a_posdef_chk)

    if display_flag
        if !isposdef(K)
            println("K not posdef! eigenvalues are:")
            s,Q = eigen(K)
            display(sort(s))
        end
    end

    return score
end

function evalconditioningscore( U::Matrix{T},
                                display_flag::Bool,
                                a_posdef_chk::T,
                                y::Vector{T},
                                c::Vector{T})::T where T

    score = discrepancyscorefromdata(U, y, c; a_cond = a_posdef_chk)

    if display_flag
        if !isposdef(U)
            println("U not posdef! eigenvalues are:")
            s,Q = eigen(U)
            display(sort(s))
        end
    end

    return score
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
"""
    Klobuchar

Struct holding the ionospheric correction terms for the Klobuchar Ionospheric
Model.

Fields:

- `α₀`
- `α₁`
- `α₂`
- `α₃`
- `β₀`
- `β₁`
- `β₂`
- `β₃`
"""
struct Klobuchar
    α₀::Float64
    α₁::Float64
    α₂::Float64
    α₃::Float64
    β₀::Float64
    β₁::Float64
    β₂::Float64
    β₃::Float64
end

"""
    EphemerisEpoch

Struct for the epoch time of the ephemeris.

Fields:

- `timestamp::String`: UTC timestamp string
- `year::Float64`: year
- `month::Float64`: month
- `day::Float64`: day
- `hour::Float64`: hour
- `minute::Float64`: minute
- `second::Float64`: fractional second
"""
struct EphemerisEpoch
    timestamp::String
    year::Float64
    month::Float64
    day::Float64
    hour::Float64
    minute::Float64
    second::Float64
end

"""
    BRDC

Struct holding GPS Broadcast Ephemeris parameters and ionospheric correction
terms.

Fields:

- `prn::Int64`: satellite assigned PRN
- `iono_parms::T`: struct containing ionospheric correction terms
- `epoch::EphemerisEpoch`: ephemeris epoch time structure
- `Af₀::Float64`: clock bias in seconds
- `Af₁::Float64`: clock drift in seconds/second
- `Af₂::Float64`: clock drift rate in seconds/second²
- `IODE::Float64`: issue of data
- `Crs::Float64`:  amplitude of the sine harmonic correction
    * apply to orbital radius
- `Crc::Float64`: amplitude of the Cosine Harmonic Correction
    * use with orbt radius
- `Cus::Float64`: amplitude of the Sine Harmonic Correction
    * apply to Argument of Latitude
- `Cuc::Float64`: amplitude of the Cosine Harmonic
    * apply to Argument of Latitude
- `Cis::Float64`: amplitude of Sine Harmonic Correction
    * apply to th Angle of Inclination (i)
- `Cic::Float64`:  amplitude of Cosine Harmonic Correction
    * apply to th Angle of Inclination (i)
- `ṅ::Float64`: mean motion difference from computed value in rad/s
- `M₀::Float64`: mean anomaly at reference time in rad
- `e::Float64`: eccentricity
- `a::Float64`: semi-major axis in meters
- `Toe::Float64`: referemce time of ephemeris (seconds into GPS week)
- `Toc::Float64`: time of clock
- `Ω::Float64`: Longitude of Ascending Node of orbit plane
- `i::Float64`: inclination angle at reference time in rad
- `di::Float64`: rate of change inclination angle in rad/s
- `dα::Float64`: rate of change of Right Ascension in rad/s
- `ω::Float64`: argument of perigee in rad
- `gps_week::Float64`: GPS week number (use with Toe)
- `health::Float64`: satellite health (0.00 = usable)
"""
struct BRDC{T}
    prn::Int64
    iono_parms::T
    epoch::EphemerisEpoch
    Af₀::Float64
    Af₁::Float64
    Af₂::Float64
    IODE::Float64
    Crs::Float64
    Crc::Float64
    Cus::Float64
    Cuc::Float64
    Cis::Float64
    Cic::Float64
    ṅ::Float64
    M₀::Float64
    e::Float64
    a::Float64
    Toe::Float64
    Toc::Float64
    Ω::Float64
    i::Float64
    di::Float64
    dα::Float64
    ω::Float64
    gps_week::Float64
    health::Float64
end

"""
    replaceD(str)

Replaces a `D` with an `E` from values in RINEX files so they can be parsed
using `parse`.

Required Arguments:

- `str`: string from RINEX file that contains `D` in place of the `E`

Returns:

- `str`: string with all `D` characters converted to `E` which can be parsed
         into `Float64`
"""
function replaceD(str)
    str = replace(str, "D" => "E")
    return str
end

"""
    loadrinex(file_name)

Loads a RINEX broadcast ephemeris file and stores into dictionary where the keys
are the PRN numbers.

Julia implementation based off the MATLAB implementation, writen by Ben K.
Bradley (2009-07-19).

Returns dictionary epoch dictionaries containing `BRDC` structs for each PRN.
The dictionary keys are the PRN values.

Required Arguments:

- `file_name`: name of file

Returns:

- `BRDCs::Dict`
"""
function loadrinex(file_name)
    f = open(file_name, "r")
    # Read header
    endofheader = false
    alphas = missing
    betas = missing
    fline = missing
    while ~endofheader
        fline = readline(f)
        # Check if this is the last line of the header
        endofheader = occursin("END OF HEADER", fline)
        # Check if this is the line with the ionospheric correction terms, αᵢ
        if occursin("ION ALPHA", fline)
            alphas = split(fline)
            alphas = [parse(Float64, replaceD(alphas[1])),
                      parse(Float64, replaceD(alphas[2])),
                      parse(Float64, replaceD(alphas[3])),
                      parse(Float64, replaceD(alphas[4]))]
        end
        # Check if this is the line with the ionospheric correction terms, βᵢ
        if occursin("ION BETA", fline)
            betas = split(fline)
            betas = [parse(Float64, replaceD(betas[1])),
                     parse(Float64, replaceD(betas[2])),
                     parse(Float64, replaceD(betas[3])),
                     parse(Float64, replaceD(betas[4]))]
        end
    end
    iono_parms = Klobuchar(alphas[1], alphas[2], alphas[3], alphas[4],
                           betas[1], betas[2], betas[3], betas[4])
    BRDCs = Dict()
    for i in 1:32
        BRDCs[i] = Dict()
    end
    fline = readline(f)
    while fline != ""
        """
            First Line of PRN Navigation Data
        """
        prn = parse(Int64, fline[1:2], base=10)       # PRN
        # Get emphemeris epoch time
        year = 2000+parse(Int64, fline[4:5], base=10) # Ephemeris year
        month = parse(Int64, fline[7:8], base=10)     # Ephemeris month
        day = parse(Int64, fline[10:11], base=10)     # Ephemeris day
        hour = parse(Int64, fline[13:14], base=10)    # Ephemeris hour
        minute = parse(Int64, fline[16:17], base=10)  # Ephemeris minute
        second = parse(Float64, fline[19:22])         # Ephemeris second (fractional)
        secondint = Int64(floor(second))              # Second truncated
        secondfrac = Int64(floor((second - secondint)*10))  # Fractional part of second
        timestamp = string(year, "-", string(month, pad=2), "-",
                           string(day, pad=2), "T", string(hour, pad=2), ":",
                           string(minute, pad=2), ":", string(secondint, pad=2),
                           ".", string(secondfrac), "Z")
        epoch = EphemerisEpoch(timestamp, year, month, day, hour, minute, second)
        # Read clock information
        Af₀ = parse(Float64, replaceD(fline[23:41]))  # Clock bias (s)
        Af₁ = parse(Float64, replaceD(fline[42:60]))  # Clock drift (s/s)
        Af₂ = parse(Float64, replaceD(fline[61:79]))  # Clock drift rate (s/s²)
        """
            Second Line of PRN Navigation Data
        """
        fline = readline(f)                           # Read line
        IODE = parse(Float64, replaceD(fline[4:22]))  # Issue of data (ephemeris)
        Crs = parse(Float64, replaceD(fline[23:41]))  # Amplitude of the Sine Harmonic Correction
                                                      # Apply to orbital radius
        ṅ = parse(Float64, replaceD(fline[42:60]))    # Mean motion difference from computed value (rad/s)
        M₀ = parse(Float64, replaceD(fline[61:79]))   # Mean anomaly at reference time (rad)
        """
            Third Line of PRN Navigation Data
        """
        fline = readline(f)                           # Read line
        Cuc = parse(Float64, replaceD(fline[4:22]))   # Amplitude of the Cosine Harmonic Correction
                                                      # Apply to Argument of Latitude
        e = parse(Float64, replaceD(fline[23:41]))    # Eccentricity
        Cus = parse(Float64, replaceD(fline[42:60]))  # Amplitude of the Sine Harmonic Correction
                                                      # Apply to Argument of Latitude
        a = parse(Float64, replaceD(fline[61:79]))^2  # Semi-major axis (meters)
        """
            Fourth Line of PRN Navigation Data
        """
        fline = readline(f)                           # Read line
        Toe = parse(Float64, replaceD(fline[4:22]))   # Referemce time of ephemeris (seconds into GPS week)
        Cic = parse(Float64, replaceD(fline[23:41]))  # Amplitude of Cosine Harmonic Correction
                                                      # Apply to th Angle of Inclination (i)
        Ω = parse(Float64, replaceD(fline[42:60]))    # Longitude of Ascending Node of orbit plane
                                                      # at weekly epoch (rad)
        Cis = parse(Float64, replaceD(fline[61:79]))  # Amplitude of Sine Harmonic Correction
                                                      # Apply to th Angle of Inclination (i)
        """
            Fifth Line of PRN Navigation Data
        """
        fline = readline(f)                           # Read line
        i = parse(Float64, replaceD(fline[4:22]))     # Inclination angle at reference time (rad)
        Crc = parse(Float64, replaceD(fline[23:41]))  # Amplitude of the Cosine Harmonic Correction
                                                      # Use with orbt radius
        ω = parse(Float64, replaceD(fline[42:60]))    # Argument of perigee (rad)
        dα = parse(Float64, replaceD(fline[61:79]))   # Rate of change of Right Ascension (rad/s)
        """
            Sixth Line of PRN Navigation Data
        """
        fline = readline(f)                           # Read line
        di = parse(Float64, replaceD(fline[4:22]))    # Rate of change inclination angle (rad/s)
        week = parse(Float64, replaceD(fline[42:60])) # GPS week number (use with Toe)
        if week < 1024
            week += 1024
        end
        """
            Seventh Line of PRN Navigation Data
        """
        fline = readline(f)                           # Read line
        health = parse(Float64, replaceD(fline[23:41])) # Satellite health (0.00 = usable)
        """
            Eigtht Line of PRN Navigation Data
        """
        fline = readline(f)                           # Read line
        Toc = Toe                                     # Time of clock
        fline = readline(f)                           # Read line
        # Save to `BRDC` struct and place into dictionary, `BRDCs`
        BRDCs[prn][timestamp] = BRDC(prn, iono_parms, epoch, Af₀, Af₁, Af₂, IODE, Crs, Crc,
                                     Cus, Cuc, Cis, Cic, ṅ, M₀, e, a, Toe, Toc, Ω, i,
                                     di, dα, ω, week, health)
    end
    return BRDCs
end


text:
# Sample script for plotting fieldlines with handpicked seeds.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot, FieldTracer
using Vlasiator: RE # Earth radius, [m]

file = "../../bulk.0000999.vlsv"
meta = load(file)
(;coordmin, coordmax, ncells) = meta
pcolormesh(meta, "proton/vg_rho", colorscale=Linear)

if ncells[2] == 1
   dim_ = (1,3)
elseif ncells[3] == 1
   dim_ = (1,2)
else
   @error "Not implemented for ncells = $ncells."
end

# regular Cartesian mesh
grid1 = range(coordmin[dim_[1]], coordmax[dim_[1]], length=ncells[dim_[1]])
grid2 = range(coordmin[dim_[2]], coordmax[dim_[2]], length=ncells[dim_[2]])

ns1 = 14 # number of seeding points on the dayside
ns2 = 10 # number of seeding points on the nightside
ns4 = 4  # number of seeding points in the polar region
seeds = Matrix{Float64}(undef, 2, ns1+2*ns2+2*ns4)
for i in 1:ns1
   seeds[1,i] = coordmax[dim_[1]] / ns1 * (i - 1)
   seeds[2,i] = 0.0
end
for i in 1:ns2
   seeds[1,ns1+i] = -10RE
   seeds[2,ns1+i] = coordmin[dim_[2]] + (coordmax[dim_[2]] - coordmin[dim_[2]]) / ns2 * (i - 1)
   seeds[1,ns1+ns2+i] = -30RE
   seeds[2,ns1+ns2+i] = coordmin[dim_[2]] + (coordmax[dim_[2]] - coordmin[dim_[2]]) / ns2 * (i - 1)
end
for i in 1:ns4
   seeds[1,ns1+2*ns2+i] = -20RE + 10RE*(i - 1)
   seeds[2,ns1+2*ns2+i] = 30RE
   seeds[1,ns1+2*ns2+ns4+i] = -20RE + 10RE*(i - 1)
   seeds[2,ns1+2*ns2+ns4+i] = -30RE
end

b = meta["vg_b_vol"]
b1 = reshape(b[dim_[1],:], ncells[dim_[1]], ncells[dim_[2]])
b2 = reshape(b[dim_[2],:], ncells[dim_[1]], ncells[dim_[2]])

for i = axes(seeds,2)
   startx, starty = seeds[:,i]
   x1, y1 = trace(b1, b2, startx, starty, grid1, grid2;
      ds=0.5, maxstep=3000, gridtype="ndgrid")
   x1 ./= RE
   y1 ./= RE
   if length(x1) < 5; continue; end
   line = plot(x1, y1, color="w")
   add_arrow(line[1])
end


text:
using BenchmarkTools
using Plots
using Statistics

# Example 1 - Arbitrary function

f(x) = 2x*sin(10x)*exp(-x)+10

function example_for_loop(arr)
    results = similar(arr)
    @inbounds for i in eachindex(arr)
        results[i] = f(arr[i])
    end
    return results
end
function example_broadcast(arr)
    results = f.(arr)
    return results
end
function example_inline_broadcast(x)
    return 2 .* x .* sin.(10 .* x).*exp.(-x) .+ 10
end

function example_for_loop_simd(arr)
    results = similar(arr)
    @inbounds @simd for i in eachindex(arr)
        results[i] = f(arr[i])
    end
    return results
end

function example_for_loop_threaded(arr)
    results = similar(arr)
    @inbounds Threads.@threads for i in eachindex(arr)
        results[i] = f(arr[i])
    end
    return results
end

using Distributed
@everywhere f(x) = 2x*sin(10x)*exp(-x)+10
function example_mapping_distributed(arr)
    results = pmap(f, arr)
    return results
end

using CUDA

function example_cuda_implementation(arr)
    arr_gpu = cu(arr) # copy array to the GPU
    res_gpu = f.(arr_gpu)
    res_cpu = Array(res_gpu)
    return res_cpu
end

function example_cuda_kernel!(y, x)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(x)
        @inbounds y[i] = f(x[i])
    end
    nothing
end

function run_example_kernel()
    x = cu(rand(2^26))
    y = similar(x)
    kernel = @cuda launch=false example_cuda_kernel!(y, x)
    config = launch_configuration(kernel.fun)
    threads = min(length(x), config.threads)
    blocks = cld(length(x), threads)

    time_taken = @elapsed CUDA.@sync kernel(y, x; threads, blocks)

    x_cpu = Array(x)
    y_cpu = f.(x_cpu)

    correctness = Array(y) .≈ y_cpu
    return time_taken, correctness
end

const non_distributed_funcs = [example_for_loop, example_broadcast, example_inline_broadcast, example_for_loop_simd, example_for_loop_threaded]
const distributed_funcs = [example_mapping_distributed]
const cuda_funcs = [example_broadcast]

function benchmark_funcs(funcs; use_gpu=false)
    ns = 2 .^ (2:2:12)
    result_times = Dict{String, Any}()

    labels = (x->"$(string(Symbol(x)))").(funcs)
    arrays = [rand(n) for n in ns]
    if use_gpu
        arrays .= cu.(arrays) # copy the array to the GPU
    end
    for (label, fn) in zip(labels, funcs)
        r = Dict{Symbol, Any}()
        min_times = []
        std_times = []
        for (n, arr) in zip(ns, arrays)
            if use_gpu
                benchmark_results = @benchmark CUDA.@sync $fn($arr)
            else
                benchmark_results = @benchmark $fn($arr)
            end

            min_time = minimum(benchmark_results.times)
            std_time = std(benchmark_results.times)
            push!(min_times, min_time)
            push!(std_times, std_time)
        end
        result_times[label] = Dict(:min_times => min_times, :std_times=>std_times)
    end
    return ns, result_times
end

function plot_results(ns, result_times; new_plot=true)
    initial_plot_fn = new_plot ? plot : plot!
    plt = initial_plot_fn(size=(500,400), dpi=600)
    for (label, res) in result_times
        plot!(ns, res[:min_times]; label=label, markershape=:auto)
    end
    plot!(;legend=:bottomright,xscale=:log10, yscale=:log10, lw=4)
    xlabel!("n")
    ylabel!("Time (ns)")
    return plt
end

max_stars_repo_name: JamieMair/intro-to-julia-talk

max_issues_repo_name: JamieMair/intro-to-julia-talk

max_forks_repo_name: JamieMair/intro-to-julia-talk


text:
@testset "native_function.jl" begin
    @testset "signature parser" begin
        @testset "$example" for example in (
            "arg0:f32",
            "arg0:f32,arg1:f32",
            "arg0:i64,arg1:i32",
            "arg0:f32[10]",
            "?arg0:i32,arg1:f32[arg0]",
            "arg2:f32[arg0]",
            "?arg0:i32,?arg1:i32,arg2:f32[arg0,arg1]",
            "arg3:f32[arg1,arg0]",
            "arg0:f32,?arg1:i32,arg2:f32[arg1]"
        )
            # This is just a quick check to make sure the parser doesn't error.
            # later integration tests will show it has the right behavour.
            @test DexCall.parse_sig(example) isa Vector{DexCall.Binder}
        end
    end

    @testset "dex_func anon funcs" begin
        @test dex_func"\x:Float. exp x"(0f0) === 1f0
        @test dex_func"\x:Float. (2.0*x, x)"(1.5f0) === (3f0,  1.5f0)
        @test dex_func"\x:Int64 y:Int. I64ToI x + y"(Int64(1), Int32(2)) === Int32(3)
        @test dex_func"\x:((Fin 3)=>Float). sum x"([1f0, 2f0, 3f0]) === 6f0

        @test dex_func"\x:((Fin 3)=>Float). for i. 2.0 * x.i"([1f0, 2f0, 3f0]) isa Vector{Float32}
        @test dex_func"\x:((Fin 3)=>Float). for i. 2.0 * x.i"([1f0, 2f0, 3f0]) == [2f0, 4f0, 6f0]
    end

    @testset "dex_func named funcs" begin
        dex_func"""
        def myTranspose (n: Int) ?-> (m: Int) ?->
                        (x : (Fin n)=>(Fin m)=>Float) : (Fin m)=>(Fin n)=>Float =
            for i j. x.j.i
        """

        myTranspose([1f0 2f0 3f0; 4f0 5f0 6f0]) isa AbstractMatrix{Float32}
        @test myTranspose([1f0 2f0 3f0; 4f0 5f0 6f0]) == [1f0 2f0 3f0; 4f0 5f0 6f0]'

        dex_func"double_it = \x:Float. 2.0 * x"
        @test double_it(4f0) === 8f0
    end

    @testset "dex_func not all implicits at start" begin
        dex_func"def f (a : Float) (n : Int) ?-> (b : (Fin n)=>Float) : Float = a + sum b"
        @test f(100f0, [10f0, 2f0, 0.1f0, 0.1f0, 0.1f0]) === 112.3f0
    end

    @testset "dex_func named const funcs" begin
        @eval dex_func"foo = \x:Int. 1.5"c  # use @eval to run at global scope, so can declare const
        @test isconst(@__MODULE__, :foo)
        @test foo(Int32(4)) === 1.5f0
    end

    @testset "dex_func errors" begin
        @test_throws ArgumentError dex_func"\x:Float. exp x"(0.0)
    end

    @testset "NativeFunction directly" begin
        m = DexModule(raw"def addTwo (n: Int) ?-> (x: (Fin n)=>Float) : (Fin n)=>Float = for i. x.i + 2.0")
        add_two = NativeFunction(m.addTwo)
        @test add_two([1f0, 10f0]) == [3f0, 12f0]
    end
end

max_stars_repo_name: justinchiu/dex-lang

max_issues_repo_name: justinchiu/dex-lang

max_forks_repo_name: justinchiu/dex-lang


text:
function mulliken_charges(scf::SCF)
    mulliken_charges = copy(scf.mole.nuclei.charges)
    orb_charges = sum(2 .* scf.P .* scf.S, dims=2)
    basis_idx = get_basis_idx(scf.mole)
    for (i, q) in enumerate(orb_charges)
        mulliken_charges[basis_idx[i]] -= q
    end
    return mulliken_charges
end


text:
"""
Provides a really simple implementation of Bayesian Optimization
for black box function. The current implementation is loosely based
off of [this](https://github.com/fmfn/BayesianOptimization) python library.

TODO: This could probably be setup to implement the
[MathProgBase](https://github.com/JuliaOpt/MathProgBase.jl)
SolverInterface, but for now I'm going to leave this as is.
"""
module BayesianOptimization

using Logging
using Distributions
using GaussianProcesses
using NLopt

export BayesOpt, optimize!, optimize

# Logging.configure(output=STDOUT, level=DEBUG)

"""
A type for storing our search history, bounds and model.
"""
type BayesOpt
    """
    A function which takes the parameters
    as input and returns a value which is to
    be maximized."""
    f::Function

    """
    The history of parameters tried as an N x M matrix
    where N is the number of parameters and M is the number
    observations.
    """
    X::Array

    """ The history of returned values to function `f` """
    y::Array

    """ The best set of parameters tried """
    xmax::Array

    """ The best result found """
    ymax::Float64

    """
    A collection of bounds for each parameter to restrict the search.
    Each bound can be an abstract array, but is typically a Range{T},
    """
    bounds::Array{AbstractArray}

    """
    The internal model of the dataset which can predict the next value.

    Currently, only a [GaussianProcess](https://github.com/STOR-i/GaussianProcesses.jl)
    is used.
    """
    model::Any
end

"""
The default BayesOpt constructor

Args:
- `f`: the function to optimize
- `bounds`: a collection of bounds (1 per parameter) to direct the search.

Returns:
- `BayesOpt`: a new `BayesOpt` instance.
"""
function BayesOpt{T<:AbstractArray}(f::Function, bounds::T...; noise=-1e8)
    # Initialize our parameters matrix
    X = Array{Float64}(length(bounds), 1)

    # Initialize our results array
    y = Array{Float64}(1)

    # Set our initial observations
    # requiring use to run the function once
    xmax = [rand(b) for b in bounds]
    ymax = f(xmax...)

    X[:,1] = xmax
    y[1] = ymax

    # return our new BayesOpt instance with these reasonable defaults
    # NOTE: see https://github.com/STOR-i/GaussianProcesses.jl for our
    # Gaussian Process arguments.
    BayesOpt(
        f, X, y, xmax, ymax, collect(bounds),
        GP(X[:,1:1], y[1:1], MeanZero(), SE(0.0, 0.0), noise)
    )
end

"""
Performs the optimization loop on a `BayesOpt` instance for the
specified number of iterations with the provided number of restarts.

Args:
- `opt`: the `BayesOpt` type to optimize
- `iter`: the number of iterations to run the optimization loop
- `restarts`: the number of restarts to allow in the acquisition function.
- `noise` : log of the observation noise, default no observation noise

Returns:
- `Array{Real}(nparams, 1)`: the best params found for `f`
- `Real`: the best result found.
"""
function optimize!(opt::BayesOpt, iter=100, restarts=10; noise=-1e8)
    # The last index we looked at is equal to the current number of
    # observation looked at.
    last_index = length(opt.y) + 1

    # Allocate `iter` number of observation to both `X` and `y`
    opt.X = hcat(opt.X, Array{Float64}(length(opt.bounds), iter))
    append!(opt.y, Array{Float64}(iter))

    # Begin our optimization loop
    for i in last_index:length(opt.y)
        Logging.debug("iteration=$i, x=$(opt.xmax), y=$(opt.ymax)")

        # Select the next params
        new_x = acquire_max(opt.model, opt.ymax, restarts, opt.bounds)

        # Run the next params
        new_y = opt.f(new_x...)

        Logging.debug("new x=$new_x, y=$new_y")

        # update X and y with next params and run
        opt.X[:,i] = new_x
        opt.y[i] = new_y

        # update the observations in the internal model (ie: the GP).
        # NOTE: This update process is specific to Gaussian Processes,
        # but could be refactored out.
        # opt.model.x = opt.X[:,1:i+1]
        # opt.model.y = opt.y[1:i+1]
        # GaussianProcesses.update_mll!(opt.model)

        # Just rebuild the model cause apparently the above complains if we update with more data.
        opt.model = GP(opt.X[:,1:i], opt.y[1:i], MeanZero(), SE(0.0, 0.0), noise)

        # Update the max x and y if our new
        # result is better.
        if new_y > opt.ymax
            opt.ymax = new_y
            opt.xmax = new_x
        end
    end

    Logging.info("Optimization completed: x=$(opt.xmax), y=$(opt.ymax)")

    # return our best parameters and result respectively.
    return opt.xmax, opt.ymax
end

"""
A helper function which handles creating the `BayesOpt` instance and
running `optimize!`

Args:
- `f`: the function to optimize
- `bounds`: a collection of bounds (1 per parameter) to direct the search.
- `iter`: the number of iterations to run the optimization loop
- `restarts`: the number of restarts to allow in the acquisition function.
- `noise` : log of the observation noise, default no observation noise

Returns:
- `Array{Real}(nparams, 1)`: the best params found for `f`
- `Real`: the best result found.
- `BayesOpt`: a new `BayesOpt` instance.
"""
function optimize{T<:AbstractArray}(f::Function, bounds::T...; iter=100, restarts=10, noise=-1e8)
    # Create out BayesOpt instance
    opt = BayesOpt(f, bounds...; noise=noise)

    # call `optimize!` for it.
    optimize!(opt, iter, restarts; noise=noise)

    # return the x and y maximums like `optimize!` along with the
    # BayesOpt instance.
    return opt.xmax, opt.ymax, opt
end

"""
Acquire function wraps a tradition [acquisition function](http://arxiv.org/pdf/1012.2599.pdf)
which selects the next x value(s) that will produce the minimum
value of y balancing exploration vs exploitation.

NOTE: While currently `expected_improvement` is the only option, future versions will support
different acquisition function described in the linked paper.

Args:
- `model`: the model which the acquisition function (in this case EI) needs to use.
- `

Returns:

"""
function acquire_max(model, ymax, restarts, bounds)
    # Selective some random parameters to check
    x_max = [rand(b) for b in bounds]

    # Create an NLopt.Opt instance using LBFGS
    # See https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
    # for details.
    opt = Opt(:LD_LBFGS, length(bounds))

    # Set the upper and lower bounds
    lower_bounds!(opt, Float64[minimum(b) for b in bounds])
    upper_bounds!(opt, Float64[maximum(b) for b in bounds])

    # TODO: Set the step size if a bound is a Range{T} with a step

    # set the objective to minimize the acquisition function.
    min_objective!(opt, expected_improvement(model, ymax))

    # initial max found is 0
    ei_max = 0

    # Optimize the acquisition with different
    # random start parameters for some number of restarts
    # to avoid being stuck in a local minima.
    for i in 1:restarts
        # Sample some points at random.
        x_try = Float64[rand(b) for b in bounds]

        #Find the minimum of minus the acquisition function
        (minf,minx,ret) = NLopt.optimize(opt, x_try)

        #Store it if better than previous minimum(maximum).
        if -minf >= ei_max
            x_max = minx
            ei_max = -minf
        end
    end

    # Return the next best expected set of parameters.
    return x_max
end

"""
The expected improvement acquisition function which balances
exporation vs exploitation by selecting the next x that is not only
the most likely to be improve the current best result, but also
to account for the expected magnitude.

NOTE: Since, all function used with NLopt only take some
data and a gradient (which in this case will always be empty
cause we are using a derivative free optimize method) we use a closure
to generate this function with the model and max y value.

Args:
- `model`: any model that returns the mean and variance when predict is
        called on some sample data.
- 'ymax': the current maximum value

Returns:
- `Function`: The `ei` function to pass to `NLopt.optimize`.
"""
function expected_improvement(model, ymax)
    """
    Inner closure fuction that performs the EI.

    Args:
    - `x`: the observation to pass to the models predict.
    - `grad`: the gradient, which should be an empty array since we are using a
            derivative free optimization method, but is required by NLopt.
    """
    function ei(x, grad)
        # Get the mean and variance from the existing observations
        # by calling predict on a gaussion process.
        mean, var = predict(model, reshape(x, length(x), 1))

        if var == 0
            # If the variance is 0 just return 0
            return 0
        else
            # otherwise calculate the ei and return that.
            Z = (mean - ymax)/sqrt(var)
            res = -((mean - ymax) * cdf(Normal(), Z) + sqrt(var) * pdf(Normal(), Z))[1]
            return res
        end
    end
    return ei
end

end # module


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

using RadiationDetectorSignals
using Test

using ArraysOfArrays, FillArrays, StructArrays, Unitful

@testset "detector_waveforms" begin
    nwf = 50
    wfdata = nestedview(rand(128, nwf))
    timedata = Fill(0:0.1:12.7,nwf)

    @test @inferred(ArrayOfRDWaveforms((wfdata, timedata))) isa StructArray
    A = ArrayOfRDWaveforms((timedata, wfdata))

    @test A.time[1] == A[1].time == 0:0.1:12.7
    @test A.value isa ArrayOfSimilarArrays
    @test A.value[1] == A[1].value
end # testset

@testset "detector_waveforms with units" begin
    nwf = 50
    wfdata = nestedview(rand(128, nwf) * u"eV")
    timedata = Fill((0:0.1:12.7) * u"ns",nwf)

    @test @inferred(ArrayOfRDWaveforms((wfdata, timedata))) isa StructArray
    A = ArrayOfRDWaveforms((timedata, wfdata))

    @test A.time[1] == A[1].time == (0:0.1:12.7) * u"ns"
    @test A.value isa ArrayOfSimilarArrays
    @test A.value[1] == A[1].value
end # testset


text:
@testset "robots" begin
    F.def_LOCAL_VARS!()
    F.add_disallow_item()
    @test F.url_curpage() in F.DISALLOW
end


text:
using RandomizedLasso
using Test

@testset "RandomizedLasso.jl" begin
    # Write your tests here.
end


text:
using NaiveBayes
using Test

@testset "NaiveBayes.jl" begin

end


text:
function convertMathMode!(name::String)
    # expects markdown files as input
    @assert getFileEnding(name)=="md" "no markdown file provided"
    lines = Vector{String}
    open(name) do file
        lines = readlines(file)
    end
    rm(name)
    fout = open(name,"w")
    first = true
    for (i,line) in enumerate(lines)
        if line=="\$\$"
            if first==true
                lines[i] = "```math"
                first=false
            else
                lines[i] = "```"
                first=true
            end
        end
        print(fout,lines[i]*"\n")
    end
    close(fout)
    return name
end

function convertJuliaMode!(name::String,envname::String)
    # expects markdown files as input
    @assert getFileEnding(name)=="md" "no markdown file provided"
    lines = Vector{String}
    open(name) do file
        lines = readlines(file)
    end
    rm(name)
    fout = open(name,"w")
    first = true
    for (i,line) in enumerate(lines)
        if line=="```julia"
            lines[i] = "```@example $envname"
        end
        print(fout,lines[i]*"\n")
    end
    close(fout)
    return name
end

function removeComments!(name::String)
    # expects julia file as input, i.e. name = "xxx.jl"
    @assert getFileEnding(name)=="jl" "no julia file provided"
    lines = Vector{String}
    open(name) do file
        lines = readlines(file)
    end
    rm(name)
    fout = open(name,"w")
    remove = zeros(Int64,0)
    keep = zeros(Int64,0)
    for (i,line) in enumerate(lines)
        # display(line)
        if length(line)>0 && Int(line[1])!=35
            # Int()==35 --> character is a hashtag
            print(fout,lines[i]*"\n")
        end
    end
    close(fout)
end

function addSetupAndMakeMarkdown(name::String,envname::String="mysetup")
    # expects julia file as input, i.e. name = "xxx.jl"
    @assert getFileEnding(name)=="jl" "no julia file provided"
    newname = removeFileEnding(name)
    newname = "$(newname)_setup.md"
    # s = split(name,"")
    # nameind = findfirst(x->x==".",s)
    lines = Vector{String}
    open(name) do file
        lines = readlines(file)
    end
    rm(name)
    fout = open(newname,"w")
    print(fout,"```@setup "*envname*"\n")
    [ print(fout,lines[i]*"\n") for i=1:length(lines) ]
    print(fout,"```\n")
    close(fout)
    return newname
end
"""
    notebook2markdown(fname::String)
Workflow to get markdown file from Julia notebook
"""
function notebook2markdown(fname::String)
    @assert getFileEnding(fname)=="ipynb" "no ipynb file provided"
    run(`jupyter-nbconvert --ClearOutputPreprocessor.enabled=True --inplace $fname --to markdown`) # ipynb to markdown
    # remove file ending (if any)
    name = removeFileEnding(fname)
    file_text = convertMathMode!("$name.md")
end

"""
    notebook2markdown_repl(fname::String;envname::String="mysetup")
Workflow to get REPL-based markdown file from Julia notebook
"""
function notebook2markdown_repl(fname::String;envname::String="mysetup")
    @assert getFileEnding(fname)=="ipynb" "no ipynb file provided"
    # remove file ending (if any)
    name = removeFileEnding(fname)
    # create ```setup mysetup ``` entry
    run(`jupyter-nbconvert $fname --ClearOutputPreprocessor.enabled=True --inplace`)
    run(`jupyter-nbconvert $fname --to python`) # ipynb to python
    run(`mv $name.py $name.jl`) # python to julia
    removeComments!("$name.jl")
    file_setup = addSetupAndMakeMarkdown("$name.jl",envname)

    # convert $$...$$ to ```math...```
    run(`jupyter-nbconvert $fname --to markdown`) # ipynb to markdown
    file_text = convertMathMode!("$name.md")
    file_text = convertJuliaMode!(file_text,envname)
    # append the two files
    text = [ readlines(file_setup); readlines(file_text) ]
    rm(file_text)
    rm(file_setup)
    fout = open(name*".md","w")
    [ print(fout,line*"\n") for line in text]
    close(fout)
end

function removeFileEnding(name::String)
    s = split(name,"")
    dotindex = findfirst(x->x==".",s)
    dotindex==nothing ? (return name) : (return name[1:dotindex-1])
end

function getFileEnding(name::String)
    s = split(name,"")
    dotindex = findfirst(x->x==".",s)
    dotindex==nothing ? (return name) : (return lowercase(name[dotindex+1:end]))
end

# Workflow to get REPL-executed markdown file from Julia notebook
    # 1. Create setup
    #     jupyter-nbconvert file --to python
    #     mv file.py file.jl
    #     removeComments!()
    #     add ```@setup```
    # 2. Add math-converted markdown
    #     jupyter-nbconvert file --to markdown
    #     convert2math
    # 3. concatenate results from 1. and 2.


text:
function redGauss(A::Array{Float64,2})
    # First, we make sure to convert all to a Float64 array
    # A = convert(Array{Float64,2},A)
    numRow, numCol = size(A); # Then retrieve some important variable names
    #
    for k in 1:numCol-1
        # Is the pivot zero?
        if A[k,k] == 0.0 # Yes
            nonzero = findall(x->x!==0.0, A[:,k]) # Find all non-zero elements in the current column
            selected = nonzero[1] # We then select the first non-zero element in the column
            A[k,:], A[selected,:] = A[selected,:],A[k,:] # We swap the rows
            pivot = A[k,k]; # And make a new pivot
            #
        else # Nope, everything was fine :D
            pivot = A[k,k]; #The pivot is the k,k since we want an identity matrix
        end
        #
        A[k,:] = A[k,:]./pivot #We scale the whole current row (Pivot row)
        for j in (k+1):numRow # Last, we go from the current+1 row, up until the bottom
            A[j,:] += -A[j,k].*A[k,:] # We multiply the pivot row by the negative of the current row, and add it to the current row
        end
    end
    #
    return A
end

function lgsolve(A::Array{Float64,2})
    # First, we make sure to convert all to a Float64 array
    #A = convert(Array{Float64,2},A)
    numRow, numCol = size(A); # Then retrieve some important variable names
    #x=zeros(numRow,1);
    x = Array{Float64}(undef,numRow)
    #
    for k in 1:numCol-1
        # Is the pivot zero?
        if A[k,k] == 0.0 # Yes
            nonzero = findall(x->x!==0.0, A[:,k]) # Find all non-zero elements in the current column
            selected = nonzero[1] # We then select the first non-zero element in the column
            A[k,:], A[selected,:] = A[selected,:],A[k,:] # We swap the rows
            pivot = A[k,k]; # And make a new pivot
            #
        else # Nope, everything was fine :D
            pivot = A[k,k];
        end
        #
        A[k,:] = A[k,:]./pivot #We scale the whole current row
        for j in (k+1):numRow # Last, we go from the current+1 row, up until the bottom
            A[j,:] += -A[j,k].*A[k,:] # We multiply the pivot row by the negative of the current row, and add it to the current row
        end
    end

# Backwards solution
    x[end] = A[end,end];
    for k in reverse(1:length(x)-1)
        S=0; #initialize
        for m = 1:length(x)
            S += A[k,m]*x[m];
        end
        x[k] = A[k,end] - S; # isolate
    end
    return x
end

function redGauss2(A::Array{Float64,2})
    # First, we make sure to convert all to a Float64 array
    # A = convert(Array{Float64,2},A)
    numRow, numCol = size(A); # Then retrieve some important variable names
    #
    for k in 1:numCol-1
        # Finding the maximum value to take as a pivot
        ab = abs.(A[:,k])
        pivotIndex = findall(x->x == maximum(ab),ab) #Looking for the index of the pivot

        # Taking the first maximum, just in case it repeats
        if length(pivotIndex) > 1
            pivotInd = pivotIndex[1];
        else #No type warn
            pivotInd = pivotIndex;
        end

        pivot = A[pivotInd,k];

        # Rearranging the rows
        A[k,:], A[pivotInd,:] = A[pivotInd,:],A[k,:]

        A[k,:] = A[k,:]./pivot #We scale the whole current row
        for j in (k+1):numRow # Last, we go from the current+1 row, up until the bottom
            A[j,:] += -A[j,k].*A[k,:] # We multiply the pivot row by the negative of the current row, and add it to the current row
        end
    end
    return A
end

function lgsolve2(A::Array{Float64,2})
    # First, we make sure to convert all to a Float64 array
    # A = convert(Array{Float64,2},A);
    numRow, numCol = size(A); # Then retrieve some important variable names
    x=zeros(numRow,1);
    #
    for k in 1:numCol-1
        # Finding the maximum value to take as a pivot
        ab = abs.(A[:,k])
        pivotIndex = findall(x->x == maximum(ab),ab) #Looking for the index of the pivot

        # Taking the first maximum, just in case it repeats
        if length(pivotIndex) > 1
            pivotInd = pivotIndex[1];
        else #No type warn
            pivotInd = pivotIndex;
        end

        pivot = A[pivotInd,k];

        # Rearranging the rows
        A[k,:], A[pivotInd,:] = A[pivotInd,:],A[k,:]

        A[k,:] = A[k,:]./pivot #We scale the whole current row
        for j in (k+1):numRow # Last, we go from the current+1 row, up until the bottom
            A[j,:] += -A[j,k].*A[k,:] # We multiply the pivot row by the negative of the current row, and add it to the current row
        end
    end

    # Backwards solution
    x[end] = A[end,end];
    for k in reverse(1:length(x)-1)
        S=0; #initialize
        for m = 1:length(x)
            S += A[k,m]*x[m];
        end
        x[k] = A[k,end] - S;
    end
    return x
end


text:
function combine_vertices(g1::AbstractGraph, g2::AbstractGraph; extension_alg="cholesky", combine_function=intersect, kwargs...)
    g3 = combine_function(g1, g2)
    if ischordal(g3)
        return g3
    else
        chordal_g, data = chordal_extension(g3, extension_alg)
        return chordal_g
    end
end

function combine_clique(g1, g2; extension_alg="cholesky", combine_function=intersect, kwargs...)
    clique1 = sort(map(sort, maximal_cliques(g1)))
    clique2 = sort(map(sort, maximal_cliques(g2)))
    clique3 = combine_function(clique1, clique2)
    g3 = build_graph_from_clique(clique2)
    if ischordal(g3)
        return g3
    else
        chordal_g, data = chordal_extension(g3, extension_alg)
        return chordal_g
    end
end

const combine_functions = Dict(
    "vertices_intersect" => (g1, g2; extension_alg="cholesky") -> combine_vertices(g1, g2; extension_alg=extension_alg, combine_function=intersect),
    "vertices_union" => (g1, g2; extension_alg="cholesky") -> combine_vertices(g1, g2; extension_alg=extension_alg, combine_function=union),
    "clique_intersect" => (g1, g2; extension_alg="cholesky") -> combine_clique(g1, g2; extension_alg=extension_alg, combine_function=intersect),
    "clique_union" => (g1, g2; extension_alg="cholesky") -> combine_clique(g1, g2; extension_alg=extension_alg, combine_function=union)
)

function combine_graph(g1, g2; how="vertices_intersect", extension_alg="cholesky")
    return combine_functions[how](g1, g2; extension_alg=extension_alg)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using Plots

function render_population(trace::Trace)
  # parse the trace
  parameter = trace[:predator_prey_parameter]
  true_parameter = _trace[:predator_prey_parameter]
  measurements = trace[:measurements]

  # plot the inferred curve
  plot(
    parameter.solution;
    vars=(1, 2),
    label="inferred dynamics",
    alpha=0.5,
    color=:black,
    title="Predator-Prey population",
    xlabel="prey population",
    ylabel="predator population"
  )

  # plot the true curve
  plot!(
    true_parameter.solution;
    vars=(1, 2),
    label="true dynamics",
    alpha=0.5,
    color=:blue,
  )

  # scatter the measurements
  xs = map(p -> p[1], measurements)
  ys = map(p -> p[2], measurements)
  scatter!(xs, ys; label="measurements", color=:red)
end

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
# regularize/z-test.jl

using Test: @testset

list = [
"Aodwt.jl"
"diffl.jl"
"diffs.jl"
]

for file in list
	@testset "$file" begin
		include(file)
	end
end


text:
using ConventionalApp

proj = Project(@__DIR__)
setup_env(proj)
@include_startup(proj)

show(activeproject)

Main.include("src/SampleConventionalApp.jl")
Main.SampleConventionalApp.run_app()

max_stars_repo_name: ma-laforge/SampleConventionalApp

max_issues_repo_name: ma-laforge/SampleConventionalApp

max_forks_repo_name: ma-laforge/SampleConventionalApp


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Path simulation
################################################################################

## Path type definition

immutable Path

    ## Members
    snum::Int64 # state size
    tnum::Int64 # time points
    pathnum::Int64 # number of paths
    sample::Array{Float64, 3} # sample paths

    function Path(start::Array{Float64, 1}, disturb::Array{Float64, 4})

        ## Extract parameters
        snum = length(start);
        tnum = size(disturb)[3] + 1;
        pathnum = size(disturb)[4];
        sample = Array{Float64}(tnum, pathnum, snum);

        ## Simulate paths
        for p = 1:pathnum
            sample[1, p, :] = start;
            for t = 2:tnum
                sample[t, p, :] = disturb[:, :, t - 1, p] * sample[t - 1, p, :];
            end
        end

        return new(snum, tnum, pathnum, sample);

    end

end


text:
"""
    traits(m::Module, T::Assignable)

Returns a set of Can-types that the data type `T` exhibits.  Look through
the composite traits and return the union of all Can-types as such.
See also [`@assign`](@ref).
"""
function traits(m::Module, T::Assignable)
    traits_map = get_traits_map(m)
    base = Set{DataType}()
    for (Tmap, s) in pairs(traits_map)
        if T <: Tmap
            union!(base, s)
        end
    end
    return base
end

"""
    assign(m::Module, T::Assignable, can_type::DataType)

Assign data type `T` with the specified Can-type from a trait.
"""
function assign(m::Module, T::Assignable, can_type::DataType)
    push_traits_map!(m, T, can_type)
    return nothing
end

"""
    @assign <T> with <CanTrait1, CanTrait2, ...>

Assign traits to the data type `T`.  For example:

```julia
@assign Duck with CanFly,CanSwim
```

is translated to something like:

```julia
flytrait(::Duck) = CanFly()
swimtrait(::Duck) = CanSwim()
```

where `x` is the name of the trait `X` in all lowercase, and `T` is the type
being assigned with the trait `X`.
"""
macro assign(T::Union{Expr,Symbol}, with::Symbol, traits::Union{Expr,Symbol})
    usage = "Invalid usage: try something like: @assign Duck with CanFly,CanSwim"
    with === :with || throw(SyntaxError(usage))
    mod = __module__

    assign_impl(mod, T, traits)
end

function assign_impl(mod, T, traits)
    expressions = Expr[]
    trait_syms = traits isa Expr && traits.head == :tuple ? traits.args : [traits]
    for this_can_type in trait_syms
        trait_function = trait_func_name(mod, this_can_type)

        # Add an expression like: <trait>trait(::T) = Can<Trait>()
        # e.g. flytrait(::Duck) = CanFly()
        push!(expressions,
            Expr(:(=),
                Expr(:call, trait_function, Expr(:(::), T)),
                Expr(:call, this_can_type)))

        # e.g. BinaryTraits.assign(MyModule, Duck, CanFly)
        push!(expressions,
                quote
                    BinaryTraits.assign($mod, $T, $this_can_type)
                end)
    end

    expr = quote
        $(expressions...)
    end

    display_expanded_code(expr)
    return esc(expr)
end


text:
"""
    rank_boolean_mat(mat::Vector{BitArray{1}})

Return the rank of `mat`, when this is interpreted as a the list of rows
of a matrix over the field with two elements,
where `true` and `false` represent the identity and the zero element of
the field, respectively.
(Two rows can be added by applying `xor.` to them.)
"""
function rank_boolean_mat(mat::Vector{BitArray{1}})
    length(mat) == 0 && return 0
    ncols = length(mat[1])
    nzheads = []
    vectors = []

    for row in mat
      # Reduce the row with the known basis vectors.
      for j in 1:length(vectors)
        if row[nzheads[j]]
          row = xor.(row, vectors[j])
        end
      end

      j = findnext(row, 1)
      if j != nothing
        # We found a new basis vector.
        push!(vectors, row)
        push!(nzheads, j)
      end
    end

    return length(vectors)
end

"""
    GrayCodeSwitchIndexIterator(n::Int)

Return an iterator for the sequence of those bit positions in an `n`-bit
Gray code where the flip has to take place.
For large enough `n`, the first values are
1, 2, 1, 3, 1, 2, 1, 4, 1, 2, 1, 3, 1, 2, 1, 5, 1, 2, 1, ...
(This is [series A001511 of the OEIS](https://oeis.org/search?q=A001511).)
"""
mutable struct GrayCodeSwitchIndexIterator
    n::Int
    nn::Int
end

GrayCodeSwitchIndexIterator(n::Int) = GrayCodeSwitchIndexIterator(n, div(n, 2))

# Initialize the iterator, return 1, set the state to [1,0,...,0].
function Base.iterate(itr::GrayCodeSwitchIndexIterator)
    itr.n > 0 || return
    itr.nn > 0 || return (1, Int[])
    next = zeros(Int, itr.nn)
    next[1] = 1
    return (1, next)
end

# Define the iteration step.
# Note that `state` is changed in place.
function Base.iterate(itr::GrayCodeSwitchIndexIterator, state::Vector{Int})
    for i in 1:itr.nn
      val = state[i]
      if val == 0
        state[i] = 1
        return (2*i-1, state)
      elseif val == 1
        state[i] = 2
        return (2*i, state)
      elseif val == 2
        state[i] = 3
        return (2*i-1, state)
      else
        state[i] = 0
      end
    end

    # For odd `n > 1`, we start the second half.
    if mod(itr.n, 2) == 1 && itr.n > 1
      itr.n = 0
      return (2*itr.nn+1, state)
    end

    # The iterator is exhausted.
    return
end

"""
    NumberOfProductsInSubspace(data::Dict, J::Vector{Int}, I::Vector{Int}, bound::Int = 15)

Return either `nothing` or a pair `(e, odd)` of nonnegative big integers
that encodes the number ``2^e \\cdot odd``,
which is the cardinality of the set
``X = \\{ (x, y) \\in U \\times U; x \\cdot y \\in V \\}``,
where `V` and `U` are the subspaces of the Singer algebra described by `data`
that are spanned by the subsets of the canonical basis that are given by
`J` and `I`.

`nothing` is returned whenever the cardinality of the set
``\\{ i \\in I; B_i \\cdot B_j \\not\\in V for some j \\in I \\}``
is larger than `bound`.
"""
function NumberOfProductsInSubspace(data::Dict, J::Vector{Int}, I::Vector{Int}, bound::Int = 15)
    q = data[:parameters][1]
    monom = data[:monomials]
    monV = monom[J]
    monU = monom[I]

    # sums:
    # the nonzero monomials $b \cdot b'$ with $b, b' \in U$
    # such that $b b' \not\in V$;
    # facts:
    # the factorizations of the monomials in 'sums' inside $U \times U$
    # factset:
    # the monomials for the set K (the monomials in U that can occur
    # as factors of nonzero products outside V)
    sums = []
    facts = []
    factset = []
    for m in monU
      for mm in monU
        sum = m + mm
        if maximum(sum) < q && ! (sum in monV)
          push!(facts, [m, mm])
          union!(factset, [m, mm])
          union!(sums, [sum])
        end
      end
    end

    l = length(factset)
    if l == 0
      # all products lie in V
      return [2*size(I, 1), 1]
    elseif l > bound
      # give up
      return
    end

    # outer:
    # the monomials for the set D (those products in S that lie outside V,
    # i. e., for which we have to compute whether they vanish or not)
    outer = setdiff(sums, monV)

    # mats:
    # for each monomial m in outer, a matrix M with
    # M[i,j] = 1 if m = factset[i] * factset[j], and M[i,j] = 0 otherwise.
    mats = [[zeros(Bool, l) for j in 1:l] for i in outer]
    for pair in facts
      pos = findfirst(isequal(pair[1]+pair[2]), outer)
      if pos != nothing
        setindex!(mats[pos][findfirst(isequal(pair[1]), factset)],
                  true, findfirst(isequal(pair[2]), factset))
      end
    end

    # Compute the number of pairs in GF(2)^{|I|} \times GF(2)^{|I|}
    # that belong to elements in the space of (x,y) with x*y inside V.
    ll = length(outer)
    mmats = [[BitArray{1}(m[i]) for m in mats] for i in 1:l]
    two = BigInt(2)
    count = two^l
    mat = [BitArray{1}(zero(m[1])) for m in mats]
    for v in GrayCodeSwitchIndexIterator(l)
      mmats_v = mmats[v]
      mat = [xor.(mmats_v[i], mat[i]) for i in 1:ll]
      count = count + two^(l - rank_boolean_mat(mat))
    end

    # Split the number into 2-part and odd part.
    e = 2 * (length(I) - l)
    while iseven(count)
      e = e + 1
      count = div(count, 2)
    end
    return (e, count)
end

function ReducedMultTable(T::Matrix{Int}, subset::Vector{Int}, tozero::Vector{Int})
    @assert tozero ⊆ subset "wrong factor!"

    diff = setdiff(subset, tozero)
    lookup = fill(0, size(T, 1))
    for i in 1:length(diff)
      lookup[diff[i]] = i
    end

    ndiff = length(diff)
    redT = zeros(Int, ndiff, ndiff)
    tozero = union(tozero, 0)
    n = size(T, 1)
    for i in diff
      posi = lookup[i]
      for j in diff
        Tij = T[i,j]
        if Tij != 0
          posj = lookup[j]
          posk = lookup[Tij]
          if posk == 0
            if ! (Tij in subset)
              error("<subset> does not describe a subalgebra")
            end
          else
            redT[posi, posj] = posk
          end
        end
      end
    end

    return redT
end

# Assume the implementation from `bitarray.jl`.
function flipbit!(B::BitArray, i::Int)
    i1, i2 = Base.get_chunks_id(i)
    u = UInt64(1) << i2
    Bc = B.chunks
    Bc[i1] = xor(Bc[i1], u)
end

function RightDerivationsDimensionNonsparse(T::Matrix{Int})
    n = size(T, 1)
    nn = n^2
    A = BitArray{1}[]

    # `heads[i] == 0` means no row in `A` has leading `true` at position `i`.
    # `heads[i] == j > 0` means `A[j]` has leading `true` at position `i`.
    heads = fill(0, nn)
    for i in 1:n
      for j in i:n
        for m in 1:n
          # Create a new row.
          eqn = falses(nn)
          k = T[i,j]
          if k != 0
            flipbit!(eqn, (k-1)*n+m)
          end
          for k in 1:n
            if T[k,j] == m
              flipbit!(eqn, (i-1)*n+k)
            end
            if T[i,k] == m
              flipbit!(eqn, (j-1)*n+k)
            end
          end

          # Reduce the new row with the known ones.
          pos = findnext(eqn, 1)
          while ! isnothing(pos)
            h = heads[pos]
            if h == 0
              # extend the matrix
              push!(A, eqn)
              heads[pos] = length(A)
              break
            else
              eqn = xor.(eqn, A[h])
              pos = findnext(eqn, pos+1)
            end
          end
        end
      end
    end

    return nn - length(A)
end

function LL4QuoDerDimNonsparse(data, maxdim = 50)
    data[:LL] == 4 || return
    q = data[:parameters][1]
    z = data[:parameters][3]
    n = z+1
    T = SingerAlg.MultTable(data)
    V = setdiff(2:n, filter(i -> all(x -> T[i,x] == 0, setdiff(2:n, [n+1-i])), 1:n))
    length(V) <= maxdim || return
    radser = BasesOfRadicalSeries(data)
    u = ReducedMultTable(T, union(V, radser[2]), [n])
    return RightDerivationsDimensionNonsparse(u)
end

#########################################################################
##
##  sparse matrix variant
##
function flipbit_sparse!(S::Set{Int}, i::Int)
    if i in S
      delete!(S, i)
    else
      push!(S, i)
    end
end

function RightDerivationsDimension(T::Matrix{Int})
    n = size(T, 1)
    nn = n^2
    A = Set{Int}[]

    # `heads[i] == 0` means no row in `A` has leading `true` at position `i`.
    # `heads[i] == j > 0` means `A[j]` has leading `true` at position `i`.
    heads = fill(0, nn)
    for i in 1:n
      for j in i:n
        for m in 1:n
          # Create a new row.
          eqn = Set{Int}()
          k = T[i,j]
          if k != 0
            push!(eqn, (k-1)*n+m)
          end
          for k in 1:n
            if T[k,j] == m
              flipbit_sparse!(eqn, (i-1)*n+k)
            end
            if T[i,k] == m
              flipbit_sparse!(eqn, (j-1)*n+k)
            end
          end

          # Reduce the new row with the known ones.
          while length(eqn) != 0
            pos = minimum(eqn)
            h = heads[pos]
            if h == 0
              # extend the matrix
              push!(A, eqn)
              heads[pos] = length(A)
              break
            else
              # replace eqn by the symmetric difference of eqn and A[h]
              symdiff!(eqn, A[h])
            end
          end
        end
      end
    end

    return nn - length(A)
end

function LL4QuoDerDim(data, maxdim = 50)
    data[:LL] == 4 || return
    q = data[:parameters][1]
    z = data[:parameters][3]
    n = z+1
    T = SingerAlg.MultTable(data)
    V = setdiff(2:n, filter(i -> all(x -> T[i,x] == 0, setdiff(2:n, [n+1-i])), 1:n))
    length(V) <= maxdim || return
    radser = BasesOfRadicalSeries(data)
    u = ReducedMultTable(T, union(V, radser[2]), [n])
    return RightDerivationsDimension(u)
end

"""
```
julia> z = 171;  qs = [11, 68];

julia> datalist = [SingerAlg.LoewyStructureInfo(q, SingerAlg.OrderMod(q, z), z) for q in qs];

julia> ids = [setdiff(2:(z+1), [10, 37, 55, 64, 82, 100]),
       [73, 91, 109, 118, 136, 163, 172]];

julia> ideals = [[ids[1], ids[1]], [ids[2], ids[2]]];

julia> labels = ["I", "J"];

julia> SingerAlg.SubquoDerDim(datalist, labels, ideals, 150)

julia> SingerAlg.SubquoDerDim(datalist, labels, ideals, 158)
Dict{Symbol,Array{T,1} where T} with 2 entries:
  :labels => ["I", "J"]
  :derdim => [14179, 14201]

```
"""
function SubquoDerDim(datalist, labels, subspaces, maxdim::T = 50) where T<:Integer
    if length(Set([data[:parameters][3] for data in datalist])) != 1
      error("`datalist` must belong to algebras of the same dimension")
    end

    # Consider only those invariant subspaces that are ideals.
    ideals = Dict[]
    for i in 1:length(subspaces)
      list = subspaces[i]
      good = filter(j -> BasisOfIdeal(datalist[j], list[j]) == list[j], 1:length(list))
      if length(good) == length(list)
        push!(ideals, Dict(:label => labels[i], :subsets => list))
      elseif length(good) != 0
        error("subspace is an ideal for some q but not for all of them?")
      end
    end

    # Sort the info by increasing dimension.
    sort!(ideals, by = r -> length(r[:subsets][1]))

    # Run over pairs '[ I, J ]' of index lists with 'J \subseteq I':
    # Compute the dimension of the derivations of the factor 'A{I} / A{J}'.
    for i in 1:length(ideals)
      l = ideals[i][:subsets]
      for j in 1:(i-1)
        derdims = Int[]
        for k in 1:length(datalist)
          if issubset(ideals[j][:subsets][k], l[k]) &&
             length(l[k]) - length(ideals[j][:subsets][k]) <= maxdim
            t = SingerAlg.MultTable(datalist[k])
            t = SingerAlg.ReducedMultTable(t, l[k], ideals[j][:subsets][k])
            push!(derdims, SingerAlg.RightDerivationsDimension(t))
          else
            push!(derdims, -1)
          end
        end
        if any(x -> x != -1, derdims)
          # 'J' is a subset of 'I' for some 'k' ...
          if -1 in derdims
            # ... but not for all 'k'
            error("earlier decision by containment of inv. ideals?")
          elseif length(datalist) == 1 || length(Set(derdims)) != 1
            # We found a decision,
            # or we are just computing the invariant for one algebra.
            return Dict(:labels => [ideals[i][:label], ideals[j][:label]],
                        :derdim => derdims)
          end
        end
      end
    end

    # We did not find a subquotient that distinguishes some algebras.
    return
end

function IsValidReduction(data, k)
    Z = data[:parameters][3]
    T = SingerAlg.MultTable(data)
    for i in 1:(Z+1)
      if mod(i, k) != 1
        # Omit the antidiagonal.
        for j in 1:(Z+1-i)
          if mod(j, k) != 1
            if T[i,j] != 0
              # The multiplication is *not* zero outside, modulo socle.
              return false
            end
          end
        end
      end
    end
    return true
end

max_stars_repo_name: oscar-system/SingerAlg

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: oscar-system/SingerAlg

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: oscar-system/SingerAlg

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
using GeometricFlux
using Flux
using Flux: onehotbatch, onecold, logitcrossentropy, throttle
using Flux: @epochs
using JLD2
using Statistics
using SparseArrays
using LightGraphs.SimpleGraphs
using LightGraphs: adjacency_matrix
using Random
using DataFrames
using Gadfly

@load joinpath(pkgdir(GeometricFlux), "data/cora_features.jld2") features
@load joinpath(pkgdir(GeometricFlux), "data/cora_labels.jld2") labels
@load joinpath(pkgdir(GeometricFlux), "data/cora_graph.jld2") g

num_nodes = 2708
num_features = 1433
hidden₀ = 1024
hidden₁ = 256
target_cat = 7
epochs = 25

## Preprocessing data
train_X = Matrix{Float32}(features)
train_y = Matrix{Float32}(labels)
adj_mat = Matrix{Float32}(adjacency_matrix(g))
fg = FeaturedGraph(adj_mat)

## Model
model = Chain(GCNConv(fg, num_features => hidden₀, leakyrelu),
              Dropout(0.5),
              GCNConv(fg, hidden₀ => hidden₁, leakyrelu),
              Dropout(0.5),
              Dense(hidden₁, target_cat)
        )

## Loss
loss(x, y) = logitcrossentropy(model(x), y)
accuracy(x, y) = mean(onecold(softmax(model(x))) .== onecold(y))

## Training
ps = Flux.params(model)
train_data = [(train_X, train_y)]
opt = ADAM(0.01)
acc = Float32[]
function evalcb()
    accₜ = accuracy(train_X, train_y)
    push!(acc, accₜ)
    @show accₜ
end

@epochs epochs Flux.train!(loss, ps, train_data, opt, cb=throttle(evalcb, 10))

df = DataFrame(epoch=collect(1:epochs), acc=acc)
plot(df, x=:epoch, y=:acc, Geom.line)


text:
#ДАНО: Робот - в произвольной клетке поля. дополнительно на поле могут находиться внутрение
#перегородки прямоугольной формы, среди которых могут быть и вырожденые прямоугольники (отрезки),
#эти внутренние перегородки изолированы друг от друга и от внешней рамки.
#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки по периметру внешней рамки промакированы

function ugol(r,side1,side2)
    while !isborder(r,side1) || !isborder(r,side2)
        if !isborder(r,side1)
            move!(r,side1)
        else
            move!(r,side2)
        end
    end
end

function mramka(r,side) #обходит по периметру и ставит марки
    while isborder(r,side) == false
        move!(r,side)
        putmarker!(r)
    end
end

function ramka(r)
    ugol(r, Sud, West)
    for side in [Nord, Ost, Sud, West]
        mramka(r,side)
    end
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
using FastGaussQuadrature

"""
ScaleGaussLegendre(N, a, b)

Uses `gausslegendre` from the `FastGaussQuadrature` package, and scales it to
calculate an integral from `a` to `b` instead of -1 to 1.
"""
function ScaleGaussLegendre(N, a, b)
    x,w = gausslegendre(N)
    mid = (b+a)/2
    L = (b-a)/2

    x = L*x .+ mid
    w = L * w

    return x,w
end

"""
ScaleGaussLaguerre(N, a, scale=1.)

Uses `gausslaguerre` from the `FastGaussQuadrature` package, and scales it to
calculate an integral from `a` to Inf instead of 0 to Inf. The `scale` parameter
allows for a typical length scale to be set.
"""
function ScaleGaussLaguerre(N, a, scale=1.)
    x,w = gausslaguerre(N)
    w = scale * exp.(log.(w) .+ x)

    # This is for the generalised Gauss Laguerre - but it should really be
    # done in the original func... but I don't see how that can be done
    # cleanly in matrix form.
    # x,w = gausslaguerre(N,2.0)
    # w = @. scale * exp(log(w) + x - 2*log(x))

    x = x*scale .+ a

    return x,w
end

"""
    Quadrature(en_target, ...)

This is the workhorse of the code. It figures out how to approximate the integral in the Lippmann-Schwinger equations and evaluates the Greens function too.

The standard grid is to put a Legendre region around the divergence at
k=sqrt(2m*en) of size `sym_size`. This includes a special point at the
divergence if on-shell. If this region does not go down to k=0, another Legendre
grid covers that inner range. Finally, a Laguerre grid covers the larger k
values.

The number of points in each region can be configured with `N` and if different
values are required `N_inner`, `N_mid`, and `N_outer`. The outer grid is
configured with `outer_style`, in which a Legendre grid can be chosen instead.

An alternative style is not yet implemented but includes an additional region
just around the divergence.
"""
function Quadrature(en_target ; N=10,
                    N_inner=N,
                    N_mid=N,
                    N_outer=N,
                    sym_size=0.1, #sym_inner_frac=0.01)
                    sym_inner_frac=nothing,
                    outer_style=(:laguerre, 1.0),
                    mass = 1.,
                    ħ = 1.
                    )
    κ = sqrt(2mass * en_target) / ħ

    if real(κ) < sym_size
        sym_size = real(κ)
    end

    Gfunc(x) = 1/(en_target - ħ^2*x^2/2mass)

    regions = []

    if sym_size == 0
        # This should only happen when the energy is negative.
    else
        if sym_inner_frac == nothing
            isreal(en_target) && @assert iseven(N_mid)

            # a) create one region spanning (κ - sym_size) : (κ + sym_szie)
            x,w = ScaleGaussLegendre(N_mid, real(κ) - sym_size, real(κ) + sym_size)

            w = w .* Gfunc.(x)
            push!(regions, (x,w))
        else
            # b) Create two regions from (κ - sym_size) : (κ - sym_inner_size) and
            # similar above. Then apply the special rule above.
            error("Not implemented")
            sym_inner_size = sym_size * sym_inner_frac
        end

        if isreal(κ)
            # Got to have the divergence always
            push!(regions, ([κ], [DivergenceQuad(κ)]))
        end
    end

    # Then create a Gauss-Laguerre region above and a Gauss-Legendre region below (if sym_size != k_target).
    if sym_size != real(κ)
        a = 0
        b = real(κ) - sym_size
        x,w = ScaleGaussLegendre(N_inner, a, b)
        w = w .* Gfunc.(x)
        lower = (x,w)

        pushfirst!(regions, lower)
    end

    if outer_style[1] == :laguerre
        laguerre_scale = outer_style[2]
        x,w = ScaleGaussLaguerre(N_outer, real(κ) + sym_size, laguerre_scale)
        w = w .* Gfunc.(x)
        outer = (x,w)
        push!(regions, outer)
    elseif outer_style[1] == :legendre
        legendre_max = outer_style[2]
        x,w = ScaleGaussLegendre(N_outer, real(κ) + sym_size, real(κ) + sym_size + legendre_max)
        w = w .* Gfunc.(x)
        outer = (x,w)
        push!(regions, outer)
    else
        error("Unknown outer_style $(outer_style[1])")
    end

    # x = vcat(first.(regions)...)
    # u = vcat(last.(regions)...)
    x,u = reduce(regions) do r1,r2
        x1,u1 = r1
        x2,u2 = r2
        x = union(x1,x2)
        u = zeros(ComplexF64, size(x))

        u[indexin(x1,x)] .+= u1
        u[indexin(x2,x)] .+= u2

        (x,u)
    end

    x = Vector{Float64}(x)
    inds = sortperm(x)
    x = x[inds]
    u = u[inds]
    return x,u,regions
end

function DivergenceQuad(qi)
    mass = ħ = 1.
    u = im*pi * 1/(ħ^2/mass * qi)
end

"""
    TaylorApprox(k)

This evaluates the divergent part of the integral, using a Taylor series expansion around that point.
"""
function TaylorApprox(k3)
    @assert length(k3) == 3
    @assert diff(k3) |> x -> (x[1] == x[2])
    qi = k3[2]
    dk = diff(k3)[1]

    @info "Including approximation for the divergent part"
    mass = ħ = 1.
    prefac = 2*mass/ħ^2 * dk/qi
    # f = k_list'.^2 .* Vmat

    u = zeros(3)
    u[:,3] = -prefac/2dk
    u[:,1] = prefac/2dk
    u[:,2] = prefac/2qi

    return u
end

# This was for a test calculating the K matrix instead of the T matrix.
function KMatQuadrature(en_target ; N=10,
                        N_inner=N,
                        N_mid=N,
                        N_outer=N,
                        sym_size=0.1, #sym_inner_frac=0.01)
                        sym_inner_frac=nothing,
                        laguerre_scale=1.0)
    if isreal(en_target)
        # Need to make sure there is a symmetric region in here.

        mass = ħ = 1.
        k_target = sqrt(2mass * en_target) / ħ

        if k_target < sym_size
            sym_size = k_target
        end

        function Gfunc(x)
            denom = en_target - ħ^2*x/2mass
            denom == 0. ? 0. : 1/denom
        end

        regions = []

        # Allow for two options
        if sym_inner_frac == nothing
            @assert isodd(N_mid)
            # a) create one region spanning (k_target - sym_size) : (k_target + sym_szie)
            x,w = ScaleGaussLegendre(N_mid, k_target - sym_size, k_target + sym_size)

            w .*= Gfunc.(x)
            push!(regions, (x,w))
        else
            # b) Create two regions from (k_target - sym_size) : (k_target - sym_inner_size) and
            # similar above. Then apply the special rule above.
            error("Not implemented")
            sym_inner_size = sym_size * sym_inner_frac
        end

        # Then create a Gauss-Laguerre region above and a Gauss-Legendre region below (if sym_size != k_target).
        if sym_size != k_target
            a = 0
            b = k_target - sym_size
            x,w = ScaleGaussLegendre(N_inner, a, b)
            w .*= Gfunc.(x)
            lower = (x,w)

            pushfirst!(regions, lower)
        end

        x,w = ScaleGaussLaguerre(N_outer, k_target + sym_size, laguerre_scale)
        w .*= Gfunc.(x)
        outer = (x,w)
        push!(regions, outer)

        # x = vcat(first.(regions)...)
        # u = vcat(last.(regions)...)
        x,u = reduce(regions) do r1,r2
            x1,u1 = r1
            x2,u2 = r2
            x = union(x1,x2)
            u = zeros(size(x))

            u[indexin(x1,x)] .+= u1
            u[indexin(x2,x)] .+= u2

            (x,u)
        end

        inds = sortperm(x)
        x = x[inds]
        u = u[inds]
        return x,u,regions
    else
        error("Not implemented")

        # Maybe do a GaussLegendre region around the real value of k_target, but
        # down to zero. Then throw in a Laguerre region for the outer.

    end
end


text:
include("Visitor.jl")
include("rules/consultas.jl")
"""
    Validatequery(documentAST::Node)

Válida que el AST de la consulta no tenga errores de sintaxis. Validatequery(Parse(query::String))
"""
function Validatequery(documentAST::Node)
	queryRules=[NotExtensionOnOperation(),NotTypeOnOperation(),NotSchemaOnOperation(),FragmentSubscription(),FragmentNames(),OperationNames(),OperationAnonymous(),SubscriptionFields(),FragmentUnknowNotUsed(),FragmentCycles()]
    vi= Visitante(documentAST)
	vi.visitante.(queryRules)
	return "ok"
end


text:
using Dispersal, Test

@testset "chain growth and allee" begin
    init = [1.0 4.0 7.0;
            2.0 5.0 8.0;
            3.0 6.0 9.0]

    output = ArrayOutput(init; tspan=1:3)

    growth = ExponentialGrowth(; rate=-log(2.0), timestep=1)
    allee = AlleeExtinction(minfounders = 2.0)
    ruleset = Ruleset(Chain(growth, allee))
    sim!(output, ruleset)
    @test output[1] == init
    @test output[2] == [0.0  2.0  3.5;
                        0.0  2.5  4.0;
                        0.0  3.0  4.5]
    @test output[3] == [0.0  0.0  0.0;
                        0.0  0.0  2.0;
                        0.0  0.0  2.25]
end