text:
"""
	histogram(cmd0::String="", arg1=nothing; kwargs...)

Examines the first data column to calculate histogram parameters based on the bin-width provided.
Alternatively, show histograms of GMTimage & GMTgrid objects directly. In the GMTimage case, the
options 'auto=true' or 'thresholds=(0, 0.1)' will find the histogram bounds of UInt16 images
convenient for contrast enhancement (histogram stretch). The values represent the percentage of
countings used to estimate the boundings. The option 'zoom=true' will set 'auto=true' and show
histogram only on the region of interest.

Full option list at [`pshistogram`]($(GMTdoc)histogram.html)

Parameters
----------

- $(GMT.opt_J)
- **A** | **horizontal** :: [Type => Bool]

    Plot the histogram horizontally from x = 0 [Default is vertically from y = 0].
    ($(GMTdoc)histogram.html#a)
- $(GMT.opt_Jz)
- $(GMT.opt_B)
- **C** | **color** | **cmap** :: [Type => Str | GMTcpt]

    Give a CPT. The mid x-value for each bar is used to look-up the bar color.
    ($(GMTdoc)histogram.html#c)
- **D** | **annot** | **annotate** | **counts** :: [Type => Str | Tuple]

    Annotate each bar with the count it represents.
    ($(GMTdoc)histogram.html#d)
- **E** | **width** :: [Type => Bool]			`Arg = width[+ooffset]`

    Use an alternative histogram bar width than the default set via T, and optionally shift all bars by an offset.
    ($(GMTdoc)histogram.html#e)
- **binmethod** | *BinMethod** :: [Type => Str]			`Arg = method`

    Binning algorithm: "scott", "fd", "sturges" or "sqrt" for floating point data. "second", "minute", "hour",
    "day", "week", "month" or "year" for DateTime data.
- **F** | **center** :: [Type => Bool]

    Center bin on each value. [Default is left edge].
    ($(GMTdoc)histogram.html#f)
- **G** | **fill** :: [Type => Number | Str]

    Select filling of bars [if no G, L or C set G=100].
    ($(GMTdoc)histogram.html#g)
- **I** | **inquire** | **bins** :: [Type => Bool | :O | :o | bins=(all=true,) | bins=(no_zero=true,) ]

    Inquire about min/max x and y after binning OR output the binned array.
    ($(GMTdoc)histogram.html#i)
- **L** | **out_range** :: [Type => Str]			`Arg = l|h|b`

    Handling of extreme values that fall outside the range set by **T**.
    ($(GMTdoc)histogram.html#l)
- **N** | **distribution** | **normal** :: [Type => Str]

    Draw the equivalent normal distribution; append desired pen [0.5p,black].
    ($(GMTdoc)histogram.html#n)
- $(GMT.opt_P)
- **Q** | **cumulative** :: [Type => Bool | "r"]

    Draw a cumulative histogram. Append r to instead compute the reverse cumulative histogram.
    ($(GMTdoc)histogram.html#q)
- **R** | **region** :: [Type => Str]

    Specifies the ‘region’ of interest in (r,azimuth) space. r0 is 0, r1 is max length in units.
    ($(GMTdoc)histogram.html#r)
- **S** | **stairs** :: [Type => Str | number]

    Draws a stairs-step diagram which does not include the internal bars of the default histogram.
    ($(GMTdoc)histogram.html#s)
- **T** | **range** | **bin** :: [Type => Str]			`Arg = [min/max/]inc[+n] | file|list]`

    Defines the range of the new CPT by giving the lowest and highest z-value and interval.
    ($(GMTdoc)histogram.html#t)
- **W** | **pen** :: [Type => Str | Tuple]

    Set pen attributes for sector outline or rose plot. [Default is no outline].
    ($(GMTdoc)histogram.html#w)
- **Z** | **kind** :: [Type => Number | Str]

    Choose between 6 types of histograms.
    ($(GMTdoc)histogram.html#z)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function histogram(cmd0::String="", arg1=nothing; first=true, kwargs...)

	arg2 = nothing		# May be needed if GMTcpt type is sent in via C
	N_args = (arg1 === nothing) ? 0 : 1

	gmt_proggy = (IamModern[1]) ? "histogram "  : "pshistogram "
	(length(kwargs) == 0 && cmd0 != "") && return monolitic(gmt_proggy, cmd0, arg1, arg2)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd = ""
	opt_Z = add_opt(d, "", 'Z', [:Z :kind], (counts = "0", count = "0", freq = "1", log_count = "2", log_freq = "3",
	                                         log10_count = "4", log10_freq = "5", weights = "+w"), true, "")
	opt_T = parse_opt_range(d, "", "")[1]		# [:T :range :inc :bin]
	(isa(arg1, GMTimage) || isa(arg1, GMTgrid)) && occursin("/", opt_T) && error("here 'bin' must be a scalar")

	# If inquire, no plotting so do it and return
	opt_I = add_opt(d, "", "I", [:I :inquire :bins], (all = "_O", no_zero = "_o"))
	if (opt_I != "")
		cmd *= opt_I
		if ((r = dbg_print_cmd(d, cmd)) !== nothing)
			return (!isa(arg1, GMTimage) && opt_T != "") ? r * " -T" * opt_T : r
		end
		cmd, arg1, = read_data(d, cmd0, cmd, arg1, " ")
		if (isa(arg1, GMTimage))		# If it's an image with no bin option, default to bin=1
			arg1, cmd = loc_histo(arg1, cmd, opt_T, opt_Z)
		else
			cmd *= opt_Z
			(opt_T != "") && (cmd *= " -T" * opt_T)
		end
		cmd = parse_V(d, cmd)
		return gmt(gmt_proggy * cmd, arg1)
	end

	cmd, _, _, opt_R = parse_BJR(d, cmd, "histogram", O, " -JX14c/14c")
	cmd = parse_JZ(d, cmd)[1]
	cmd = parse_common_opts(d, cmd, [:UVXY :JZ :c :e :f :p :t :params], first)[1]
	cmd = parse_these_opts(cmd, d, [[:A :horizontal], [:F :center], [:Q :cumulative], [:S :stairs]])
	cmd = add_opt_fill(cmd, d, [:G :fill], 'G')
	cmd = add_opt(d, cmd, 'D', [:D :annot :annotate :counts], (beneath = "_+b", font = "+f", offset = "+o", vertical = "_+r"))
	cmd = parse_INW_coast(d, [[:N :distribution :normal]], cmd, "N")
	(show_kwargs[1]) && print_kwarg_opts(symbs, "NamedTuple | Tuple | Dict | String")

	(GMTver > v"6.1.1") && (cmd = add_opt(d, cmd, 'E', [:E :width], (width = "", off = "+o", offset = "+o")))

	# If file name sent in, read it and compute a tight -R if this was not provided
	is_datetime = isa(arg1, Array{<:DateTime})
	(opt_R == "" && !isa(arg1, Vector{DateTime})) && (opt_R = " ")	# So it doesn't try to find the -R in next call
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R)
	cmd, arg1, arg2,  = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_args, arg1, arg2)

	# If we still do not know the bin width, either use the GMT6.2 -E or BinMethod in binmethod()
	if (opt_T == "" && !occursin(" -E", cmd) && (arg1 !== nothing) && !isa(arg1, GMTimage) && !isa(arg1, GMTgrid))
		opt_T = binmethod(d, cmd, arg1, is_datetime)	# This comes without the " -T"
		if (is_datetime)
			t = gmt("pshistogram -I -T" * opt_T, arg1)	# Call with inquire option to know y_min|max
			h = round_wesn(t[1].data)					# Only h[4] is needed
			opt_R *= sprintf("/0/%.12g", h[4])			# Half -R was computed in read_data()
			cmd *= opt_R * " -T" * opt_T
			opt_T = ""		# Clear it because the GMTimage & GMTgrid use a version without "-T" that is added at end
		end
	end

	cmd  = add_opt(d, cmd, 'L', [:L :out_range], (first = "l", last = "h", both = "b"))
	cmd *= add_opt_pen(d, [:W :pen], "W", true)     	# TRUE to also seek (lw|lt,lc,ls)
	if (!occursin("-G", cmd) && !occursin("-C", cmd) && !occursin("-S", cmd))
		cmd *= " -G#0072BD"
		!occursin("-W", cmd) && (cmd *= " -Wfaint")
    elseif (occursin("-S", cmd) && !occursin("-W", cmd))
		cmd *= " -Wfaint"
	end

	limit_L = nothing
	if (isa(arg1, GMTimage))		# If it's an image with no bin option, default to bin=1
		do_clip = (isa(arg1[1], UInt16) && (val = find_in_dict(d, [:full_histo])[1]) === nothing) ? true : false
		do_auto = ((val_auto = find_in_dict(d, [:auto :threshols])[1]) !== nothing) ? true : false	# Automatic bounds detetion
		do_getauto = ((val_getauto = find_in_dict(d, [:getauto :getthreshols])[1]) !== nothing) ? true : false
		do_zoom = ((find_in_dict(d, [:zoom])[1]) !== nothing) ? true : false	# Automatic zoom to interesting region
		(do_zoom && !do_auto) && (val_auto = nothing)		# I.e. 'zoom' sets also the auto mode
		hst, cmd = loc_histo(arg1, cmd, opt_T, opt_Z)
		if (do_clip && (all(hst[3:10,2] .== 0)))  hst[1,2] = 0; hst[2,2] = 0  end
		if (do_auto || do_getauto || do_zoom)
			which_auto = (do_auto) ? val_auto : val_getauto
			limit_L, limit_R = find_histo_limits(arg1, which_auto, 20)
			(do_getauto) && return [limit_L, limit_R]		# If only want the histogram limits, we are done then.
			if (do_zoom)
				mm = extrema(hst, dims=1)		# 1×2 Array{Tuple{UInt16,UInt16},2}
				x_max = min(limit_R * 1.15, hst[end,1])		# 15% to the right but not fall the cliff
				opt_R_ = " -R$(limit_L * 0.85)/$x_max/0/$(mm[2][2] * 1.1) "
				(opt_R != " ") && @warn("'zoom' option overrides the requested region limits and sets its own")
				cmd = replace(cmd, opt_R => opt_R_, count=1)
				opt_R = opt_R_					# It will be needed further down in vlines
			end
		end
		arg1 = hst		# We want to send the histogram, not the GMTimage
    elseif (isa(arg1, GMTgrid))
		if (opt_T != "")
			inc = parse(Float64, opt_T) + eps()		# + EPS to avoid the extra last bin at right with 1 count only
			n_bins = Int(ceil((arg1.range[6] - arg1.range[5]) / inc))
		else
			n_bins = Int(ceil(sqrt(length(arg1.z))))
			inc = (arg1.range[6] - arg1.range[5]) / n_bins + eps()
		end
		(!isa(inc, Real) || inc <= 0) && error("Bin width must be a > 0 number and no min/max")
		hst = zeros(n_bins, 2)
		@inbounds Threads.@threads for k = 1:length(arg1)  hst[Int(floor((arg1[k] - arg1.range[5]) / inc) + 1), 2] += 1  end
		[@inbounds hst[k,1] = arg1.range[5] + inc * (k - 1) for k = 1:n_bins]

		cmd = (opt_Z == "") ? cmd * " -Z0" : cmd * opt_Z
		if (!occursin("+w", cmd))  cmd *= "+w"  end		# Pretending to be weighted is crutial for the trick
		cmd *= " -T$(arg1.range[5])/$(arg1.range[6])/$inc"
		arg1 = hst		# We want to send the histogram, not the GMTgrid
	else
		(opt_T != "") && (opt_T = " -T" * opt_T)		# It lacked the -T so that it could be used in loc_histo()
		cmd *= opt_T * opt_Z
	end

	# The following looks a bit messy but it's needed to auto plotting verical lines with the limits
	show_ = false;		fmt_ = FMT[1];		savefig_ = nothing
	if (limit_L !== nothing)
		haskey(d, :show) && (show_ = (d[:show] != 0))				# Backup the :show val
		d[:show] = false
		haskey(d, :fmt) && (fmt_ = d[:fmt]; delete!(d, :fmt))		# Backup the :fmt val
		((val = find_in_dict(d, [:savefig :figname :name])[1]) !== nothing) && (savefig_ = val)
	end

	out2 = nothing;		Vd_ = 0				# Backup values
	(haskey(d, :Vd)) && (Vd_ = d[:Vd])

	# Plot the histogram
	out1 = finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1, arg2)

	# And if wished, plot the two vertical lines with the limits annotated in them
	if (limit_L !== nothing)
		if (opt_R == " ")					# Set a region for the vlines
			mm = extrema(hst, dims=1)
			opt_R = " -R$(mm[1][1])/$(mm[1][2])/0/$(mm[2][2])"
		end
		vlines!([limit_L], pen="0.5p,blue,dashed", decorated=(quoted = true, n_labels = 1, const_label = "$limit_L", font = 9, pen = (0.5, :red)), R=opt_R[4:end], Vd=Vd_)
		out2 = vlines!([limit_R], pen="0.5p,blue,dashed", decorated=(quoted = true, n_labels = 1, const_label = "$limit_R", font = 9, pen = (0.5, :red)), R=opt_R[4:end], fmt=fmt_, savefig=savefig_, show=show_, Vd=Vd_)
	end
	out = (out1 !== nothing && out2 !== nothing) ? [out1;out2] : ((out1 !== nothing) ? out1 : out2)

end

# ---------------------------------------------------------------------------------------------------
function find_histo_limits(In, thresholds=nothing, width=20)
	# Find the histogram limits of a UInt16 GMTimage that allow to better stretch the histogram
	# THRESHOLDS is an optional Tuple input containing the left and right histo thresholds, in percentage,
	# between which the histogram values will be retained. Defaults are (0.1, 0.4) percent. Note, this
	# assumes the histogram follows some sort of Gaussian distribution. It it's flat, shit occurs.
	# WIDTH is bin width used to obtain a rough histogram that is used to compute the limits.
	if (isa(In, Array{UInt16,3}) || isa(In, Array{UInt8,3}))
		L1 = find_histo_limits(view(In, :, :, 1), thresholds, width)
		L2 = find_histo_limits(view(In, :, :, 2), thresholds, width)
		L3 = find_histo_limits(view(In, :, :, 3), thresholds, width)
		return (L1[1], L1[2], L2[1], L2[2], L3[1], L3[2])
	end
	hst = loc_histo(In, "", string(width), "")[1]
	if (size(hst, 1) >= 5)
		all(hst[2:5,2] .== 0) && (hst[1,2] = 0)	# Here we always check for high counts in zero bin
	end
	max_ = maximum(hst, dims=1)[2]
	(max_ == 0) && error("This histogram had nothing but countings ONLY in first bin. No point to proceed.")
	thresh_l = 0.001;		thresh_r = 0.004
	if (isa(thresholds, Tuple) && length(thresholds) == 2)
		thresh_l, thresh_r = thresholds[:] ./ 100
	end
	thresh_l *= max_
	thresh_r *= max_
	kl = 1;		kr = size(hst, 1)
	while (hst[kl,2] == 0 || hst[kl,2] < thresh_l)  kl += 1  end
	while (hst[kr,2] == 0 || hst[kr,2] < thresh_r)  kr -= 1  end
	return Int(hst[kl,1]), Int(min(hst[kr,1] + width, hst[end,1]))
end

# ---------------------------------------------------------------------------------------------------
function loc_histo(in, cmd::String="", opt_T::String="", opt_Z::String="")
	# Very simple function to compute histograms of images (integers)
	# We put the countings in a Mx2 arrray to trick GMT (pshistogram) to think it's recieving a weighted input.
	(!isa(in[1], UInt16) && !isa(in[1], UInt8)) && error("Only UInt8 or UInt16 image types allowed here")

	inc = (opt_T != "") ? Float64(Meta.parse(opt_T)) : 1.0
	(!isa(inc, Real) || inc <= 0) && error("Bin width must be a > 0 number and no min/max")

	n_bins = (isa(in[1], UInt8)) ? 256 : Int(ceil((maximum(in) + 1) / inc))	# For UInt8 use the full [0 255] range
	hst = zeros(n_bins, 2)
	pshst_wall!(in, hst, inc, n_bins)

	cmd = (opt_Z == "") ? cmd * " -Z0" : cmd * opt_Z
	(!occursin("+w", cmd)) && (cmd *= "+w")			# Pretending to be weighted is crutial for the trick

	return hst, cmd * " -T0/$(n_bins * inc)/$inc"
end

# ---------------------------------------------------------------------------------------------------
function pshst_wall!(in, hst, inc, n_bins::Int)
	# Function barrier for type instability. With the body of this in calling fun the 'inc' var
	# introduces a mysterious type instability and execution times multiply by 3.
	if (inc == 1)
		@inbounds Threads.@threads for k = 1:length(in)  hst[in[k] + 1, 2] += 1  end
    else
		@inbounds Threads.@threads for k = 1:length(in)  hst[Int(floor(in[k] / inc) + 1), 2] += 1  end
	end
	@inbounds Threads.@threads for k = 1:n_bins  hst[k,1] = inc * (k - 1)  end
	return nothing
end

# ---------------------------------------------------------------------------------------------------
function binmethod(d::Dict, cmd::String, X, is_datetime::Bool)
	# Compute bin width for a series of binning alghoritms or intervals when X (DateTime) comes in seconds
	val = ((val_ = find_in_dict(d, [:binmethod :BinMethod])[1]) !== nothing) ? lowercase(string(val_)) : ""
	(!is_datetime) && (min_max = extrema(X))		# X should already be sorted but don't trust that
	if (val == "")
		if (!is_datetime)
			val = "scott"
		else
			min_max = extrema(X)		# X should already be sorted but don't trust that
			rng = (min_max[2] - min_max[1])
			if     (rng < 150)                 val = "second"
			elseif (rng / (60) < 150)          val = "minute"
			elseif (rng / (3600)  < 150)       val = "hour"
			elseif (rng / (86400) < 150)       val = "day"
			elseif (rng / (86400 * 7)  < 150)  val = "week"
			elseif (rng / (86400 * 31) < 150)  val = "month"
			else                               val = "year"
			end
		end
	end

	n_bins = 0.0;	bin = 0
	if     (val == "scott")   n_bins = 3.5 .* std(X) .* length(X)^(-1/3)
	elseif (val == "fd")      n_bins = 2 .* IQR(X) .* length(X)^(-1/3)
	elseif (val == "sturges") n_bins = ceil.(1 .+ log2.(length(X)))
	elseif (val == "sqrt")    n_bins = ceil.(sqrt(length(X)))
	elseif (val == "year")    bin = 86400 * 365.25
	elseif (val == "month")   bin = 86400 * 31
	elseif (val == "week")    bin = 86400 * 7
	elseif (val == "day")     bin = 86400
	elseif (val == "hour")    bin = 3600
	elseif (val == "minute")  bin = 60
	elseif (val == "second")  bin = 1
	elseif (!is_datetime)     error("Unknown BinMethod $val")
	end
	if (bin == 0)
		bin = (min_max[2] - min_max[1]) / n_bins	# Should be made a "pretty" number?
	end
	return sprintf("%.12g", bin)
end

# ---------------------------------------------------------------------------------------------------
histogram!(cmd0::String="", arg1=nothing; kw...) = histogram(cmd0, arg1; first=false, kw...)
histogram(arg1; kw...)  = histogram("", arg1; first=true, kw...)
histogram!(arg1; kw...) = histogram("", arg1; first=false, kw...)

const pshistogram  = histogram			# Alias
const pshistogram! = histogram!			# Alias


text:
function half_adder(A, B)
    S = xor(A, B)
    Cout = A & B
    Cout, S
end

function full_adder(A, B, Cin)
    T = xor(A, B)
    S = xor(Cin, T)
    Cout = (A & B) | (Cin & T)
    Cout, S
end

function nibble_adder(A4, A3, A2, A1, B4, B3, B2, B1)
    carry, S1 = half_adder(A1, B1)
    carry, S2 = full_adder(A2, B2, carry)
    carry, S3 = full_adder(A3, B3, carry)
    carry, S4 = full_adder(A4, B4, carry)

    S4, S3, S2, S1
end

function binary_multiply(input, counter)
    accumulator = 0
    input += 0b0000_0000_0000_0000
    while counter > 0
        if counter & 1 == 1
            accumulator += input
        end
        counter >>= 1
        input <<= 1
    end
    accumulator
end

binary_multiply(0b0000_1100,0b1110_1010)

ext: jl

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
using DifferentialEquations, Test

f_2dlinear = (du,u,p,t) -> (@. du = p*u)
f_2dlinear_analytic = (u0,p,t) -> @. u0*exp(p*t)
prob_ode_2Dlinear = ODEProblem(
                    ODEFunction(f_2dlinear,analytic=f_2dlinear_analytic),
                    rand(4,2),(0.0,1.0),1.01)

alg, kwargs = default_algorithm(prob_ode_2Dlinear;dt=1//2^(4))
sol =solve(prob_ode_2Dlinear;dt=1//2^(4))

@test typeof(sol.alg.algs[1]) <: Tsit5
@test typeof(sol.alg.algs[2]) <: Rosenbrock23

sol =solve(prob_ode_2Dlinear;reltol=1e-1)

@test typeof(sol.alg.algs[1]) <: Tsit5
@test typeof(sol.alg.algs[2]) <: Rosenbrock23

sol =solve(prob_ode_2Dlinear;reltol=1e-7)

@test typeof(sol.alg.algs[1]) <: Vern7
@test typeof(sol.alg.algs[2]) <: Rodas4

sol =solve(prob_ode_2Dlinear;reltol=1e-10)

@test typeof(sol.alg.algs[1]) <: Vern9
@test typeof(sol.alg.algs[2]) <: Rodas5

sol =solve(prob_ode_2Dlinear;alg_hints=[:stiff])

@test typeof(sol.alg) <: Rodas4

sol =solve(prob_ode_2Dlinear;alg_hints=[:stiff],reltol=1e-1)

@test typeof(sol.alg) <: Rosenbrock23

const linear_bigα = parse(BigFloat,"1.01")
f = (du,u,p,t) -> begin
  for i in 1:length(u)
    du[i] = linear_bigα*u[i]
  end
end
(::typeof(f))(::Type{Val{:analytic}},u0,p,t) = u0*exp(linear_bigα*t)
prob_ode_bigfloat2Dlinear = ODEProblem(f,map(BigFloat,rand(4,2)).*ones(4,2)/2,(0.0,1.0))

sol =solve(prob_ode_bigfloat2Dlinear;dt=1//2^(4))
@test typeof(sol.alg.algs[1]) <: Vern9
@test typeof(sol.alg.algs[2]) <: Rodas5

default_algorithm(prob_ode_bigfloat2Dlinear;alg_hints=[:stiff])

sol =solve(prob_ode_bigfloat2Dlinear;alg_hints=[:stiff])

@test typeof(sol.alg) <: Rodas4

sol =solve(prob_ode_bigfloat2Dlinear,nothing;alg_hints=[:stiff])

@test typeof(sol.alg) <: Rodas4

struct FooAlg end

@test_throws ErrorException solve(prob_ode_bigfloat2Dlinear,FooAlg();default_set=true)

struct FooAlg2 <: DiffEqBase.DEAlgorithm end

@test_throws ErrorException solve(prob_ode_bigfloat2Dlinear,FooAlg2();default_set=true)

prob = ODEProblem(f,rand(4,2).*ones(4,2)/2,(0.0,1.0))

sol =solve(prob;alg_hints=[:stiff])

@test typeof(sol.alg) <: Rodas4

sol =solve(prob;alg_hints=[:stiff],reltol=1e-1)

@test typeof(sol.alg) <: Rosenbrock23

sol =solve(prob;alg_hints=[:stiff],callback=CallbackSet())

@test typeof(sol.alg) <: Rodas4

prob = ODEProblem(f,rand(4,2).*ones(4,2)/2,(0.0,1.0))

alg, kwargs = default_algorithm(prob;alg_hints=[:stiff])

@test typeof(alg) <: Rodas4

m = 1.0
ω = 1.0

function mass_system!(du,u,p,t)
    # a(t) = (1/m) w^2 x
    (1/m)*(ω^2)*u[1]
end

v0 = 0.0
u0 = 1.0
tspan = (0.0,10.0)

prob = SecondOrderODEProblem(mass_system!,v0,u0,tspan)
sol = solve(prob)

@test typeof(sol.alg) <: DPRKN6


text:
## This file contains functions to handle the symetries

# A symmetry operation (symop) is a couple (Stilde, τtilde) of a
# unitary (in cartesian coordinates, but not in reduced coordinates)
# matrix Stilde and a translation τtilde such that, for each atom of
# type A at position a, Stilde a + τtilde is also an atom of type A.

# This induces a symmetry in the Brillouin zone that the Hamiltonian
# at S k is unitary equivalent to that at k, which we exploit to
# reduce computations. The relationship is
# S = Stilde'
# τ = -Stilde^-1 τtilde
# (valid both in reduced and cartesian coordinates)

# The full (reducible) Brillouin zone is implicitly represented by
# a set of (irreducible) kpoints (see explanation in docs). Each
# irreducible kpoint k comes with a list of symmetry operations
# (S,τ) (containing at least the trivial operation (I,0)), where S
# is a rotation matrix (/!\ not unitary in reduced coordinates)
# and τ a translation vector. The kpoint is then used to represent
# implicitly the information at all the kpoints Sk. The
# relationship between the Hamiltonians is
# H_{Sk} = U H_k U*, with
# (Uu)(x) = u(S^-1(x-τ))
# or in Fourier space
# (Uu)(G) = e^{-i G τ} u(S^-1 G)
# In particular, we can choose the eigenvectors at Sk as u_{Sk} = U u_k

# We represent then the BZ as a set of irreducible points `kpoints`, a
# list of symmetry operations `ksymops` allowing the reconstruction of
# the full (reducible) BZ, and a set of weights `kweights` (summing to
# 1). The value of an observable (eg energy) per unit cell is given as
# the value of that observable at each irreducible kpoint, weighted by
# kweight

# There is by decreasing cardinality
# - The group of symmetry operations of the lattice
# - The group of symmetry operations of the crystal (model.symmetries)
# - The group of symmetry operations of the crystal that preserves the BZ mesh (basis.symmetries)
# - The set of symmetry operations that we use to reduce the
#   reducible Brillouin zone (RBZ) to the irreducible (IBZ) (basis.ksymops)

# See https://juliamolsim.github.io/DFTK.jl/stable/advanced/symmetries for details.

@doc raw"""
Return the ``k``-point symmetry operations associated to a lattice and atoms.
"""
function symmetry_operations(lattice, atoms, magnetic_moments=[]; tol_symmetry=1e-5)
    symmetries = []
    # Get symmetries from spglib
    Stildes, τtildes = spglib_get_symmetry(lattice, atoms, magnetic_moments;
                                           tol_symmetry=tol_symmetry)

    for isym = 1:length(Stildes)
        S = Stildes[isym]'                  # in fractional reciprocal coordinates
        τ = -Stildes[isym] \ τtildes[isym]  # in fractional real-space coordinates
        τ = τ .- floor.(τ)
        @assert all(0 .≤ τ .< 1)
        push!(symmetries, (S, τ))
    end

    unique(symmetries)
end

"""
Filter out the symmetry operations that respect the symmetries of the discrete BZ grid
"""
function symmetries_preserving_kgrid(symmetries, kcoords)
    function preserves_grid(S)
        all(normalize_kpoint_coordinate(S * k) in kcoords
            for k in normalize_kpoint_coordinate.(kcoords))
    end
    filter(symop -> preserves_grid(symop[1]), symmetries)
end

"""
Implements a primitive search to find an irreducible subset of kpoints
amongst the provided kpoints.
"""
function find_irreducible_kpoints(kcoords, Stildes, τtildes)

    # This function is required because spglib sometimes flags kpoints
    # as reducible, where we cannot find a symmetry operation to
    # generate them from the provided irreducible kpoints. This
    # reimplements that part of spglib, with a possibly very slow
    # algorithm.

    # Flag which kpoints have already been mapped to another irred.
    # kpoint or which have been decided to be irreducible.
    kcoords_mapped = zeros(Bool, length(kcoords))
    kirreds = empty(kcoords)           # Container for irreducible kpoints
    ksymops = Vector{Vector{SymOp}}()  # Corresponding symops

    while !all(kcoords_mapped)
        # Select next not mapped kpoint as irreducible
        ik = findfirst(isequal(false), kcoords_mapped)
        push!(kirreds, kcoords[ik])
        thisk_symops = [identity_symop()]
        kcoords_mapped[ik] = true

        for jk in findall(.!kcoords_mapped)
            isym = findfirst(1:length(Stildes)) do isym
                # If the difference between kred and Stilde' * k == Stilde^{-1} * k
                # is only integer in fractional reciprocal-space coordinates, then
                # kred and S' * k are equivalent k-Points
                all(isinteger, kcoords[jk] - (Stildes[isym]' * kcoords[ik]))
            end

            if !isnothing(isym)  # Found a reducible kpoint
                kcoords_mapped[jk] = true
                S = Stildes[isym]'                  # in fractional reciprocal coordinates
                τ = -Stildes[isym] \ τtildes[isym]  # in fractional real-space coordinates
                τ = τ .- floor.(τ)
                @assert all(0 .≤ τ .< 1)
                push!(thisk_symops, (S, τ))
            end
        end  # jk

        push!(ksymops, thisk_symops)
    end
    kirreds, ksymops
end

"""
Apply a symmetry operation to eigenvectors `ψk` at a given `kpoint` to obtain an
equivalent point in [-0.5, 0.5)^3 and associated eigenvectors (expressed in the
basis of the new kpoint).
"""
function apply_ksymop(ksymop::SymOp, basis, kpoint, ψk::AbstractVecOrMat)
    S, τ = ksymop
    S == I && iszero(τ) && return kpoint, ψk

    # Apply S and reduce coordinates to interval [-0.5, 0.5)
    # Doing this reduction is important because
    # of the finite kinetic energy basis cutoff
    @assert all(-0.5 .≤ kpoint.coordinate .< 0.5)
    Sk_raw = S * kpoint.coordinate
    Sk = normalize_kpoint_coordinate(Sk_raw)
    kshift = convert.(Int, Sk - Sk_raw)
    @assert all(-0.5 .≤ Sk .< 0.5)

    # Check whether the resulting kpoint is in the basis:
    ikfull = findfirst(1:length(basis.kpoints)) do idx
        all(isinteger, basis.kpoints[idx].coordinate - Sk)
    end
    if isnothing(ikfull)
        # Build a new kpoint datastructure:
        Skpoint = build_kpoints(basis, [Sk])[1]
    else
        Skpoint = basis.kpoints[ikfull]
        @assert Skpoint.coordinate ≈ Sk
    end

    # Since the eigenfunctions of the Hamiltonian at k and Sk satisfy
    #      u_{Sk}(x) = u_{k}(S^{-1} (x - τ))
    # their Fourier transform is related as
    #      ̂u_{Sk}(G) = e^{-i G \cdot τ} ̂u_k(S^{-1} G)
    invS = Mat3{Int}(inv(S))
    Gs_full = [G + kshift for G in G_vectors(Skpoint)]
    ψSk = zero(ψk)
    for iband in 1:size(ψk, 2)
        for (ig, G_full) in enumerate(Gs_full)
            igired = index_G_vectors(basis, kpoint, invS * G_full)
            @assert igired !== nothing
            ψSk[ig, iband] = cis(-2π * dot(G_full, τ)) * ψk[igired, iband]
        end
    end

    Skpoint, ψSk
end

"""
Apply a `k`-point symmetry operation (the tuple (S, τ)) to a partial density.
"""
function apply_ksymop(symop::SymOp, basis, ρin)
    S, τ = ksymop
    S == I && iszero(τ) && return ρin
    symmetrize_ρ(basis, ρin, [symop])
end

# Accumulates the symmetrized versions of the density ρin into ρout (in Fourier space).
# No normalization is performed
function accumulate_over_symmetries!(ρaccu, ρin, basis, symmetries)
    T = eltype(basis)
    for (S, τ) in symmetries
        invS = Mat3{Int}(inv(S))
        # Common special case, where ρin does not need to be processed
        if S == I && iszero(τ)
            ρaccu .+= ρin
            continue
        end

        # Transform ρin -> to the partial density at S * k.
        #
        # Since the eigenfunctions of the Hamiltonian at k and Sk satisfy
        #      u_{Sk}(x) = u_{k}(S^{-1} (x - τ))
        # with Fourier transform
        #      ̂u_{Sk}(G) = e^{-i G \cdot τ} ̂u_k(S^{-1} G)
        # equivalently
        #      ̂ρ_{Sk}(G) = e^{-i G \cdot τ} ̂ρ_k(S^{-1} G)
        for (ig, G) in enumerate(G_vectors(basis))
            igired = index_G_vectors(basis, invS * G)
            if igired !== nothing
                @inbounds ρaccu[ig] += cis(-2T(π) * dot(G, τ)) * ρin[igired]
            end
        end
    end  # (S, τ)
    ρaccu
end

# Low-pass filters ρ (in Fourier) so that symmetry operations acting on it stay in the grid
function lowpass_for_symmetry!(ρ, basis; symmetries=basis.model.symmetries)
    for (S, τ) in symmetries
        S == I && iszero(τ) && continue
        for (ig, G) in enumerate(G_vectors(basis))
            if index_G_vectors(basis, S * G) === nothing
                ρ[ig] = 0
            end
        end
    end
    ρ
end

"""
Symmetrize a density by applying all the model symmetries (by default) and forming the average.
"""
@views function symmetrize_ρ(basis, ρin; symmetries=basis.model.symmetries)
    ρin_fourier = r_to_G(basis, ρin)
    ρout_fourier = copy(ρin_fourier)
    for σ = 1:size(ρin, 4)
        lowpass_for_symmetry!(ρin_fourier[:, :, :, σ], basis; symmetries=symmetries)
        ρout_fourier[:, :, :, σ] = accumulate_over_symmetries!(zero(ρin_fourier[:, :, :, σ]),
                                                               ρin_fourier[:, :, :, σ], basis, symmetries)
    end
    G_to_r(basis, ρout_fourier ./ length(symmetries))
end
# symmetrize the stress tensor, which is a rank-2 contravariant tensor in reduced coordinates
function symmetrize_stresses(lattice, symmetries, stresses)
    stresses_symmetrized = zero(stresses)
    for (S, τ) in symmetries
        S_reduced = inv(lattice) * S * lattice
        stresses_symmetrized += S_reduced' * stresses * S_reduced
    end
    stresses_symmetrized /= length(symmetries)
    stresses_symmetrized
end

function check_symmetric(basis, ρin; tol=1e-10, symmetries=ρin.basis.model.symmetries)
    for symop in symmetries
        @assert norm(symmetrize_ρ(ρin, [symop]) - ρin) < tol
    end
end

""""
Convert a `basis` into one that doesn't use BZ symmetry.
This is mainly useful for debug purposes (e.g. in cases we don't want to
bother thinking about symmetries).
"""
function unfold_bz(basis::PlaneWaveBasis)
    if all(length.(basis.ksymops_global) .== 1)
        return basis
    else
        kcoords = []
        for (ik, kpt) in enumerate(basis.kcoords_global)
            for (S, τ) in basis.ksymops_global[ik]
                push!(kcoords, normalize_kpoint_coordinate(S * kpt))
            end
        end
        new_basis = PlaneWaveBasis(basis.model,
                                   basis.Ecut, basis.fft_size, basis.variational,
                                   kcoords, [[identity_symop()] for _ in 1:length(kcoords)],
                                   basis.kgrid, basis.kshift, basis.symmetries, basis.comm_kpts)
    end
end

# find where in the irreducible basis `basis_irred` the kpoint `kpt_unfolded` is handled
function unfold_mapping(basis_irred, kpt_unfolded)
    for ik_irred = 1:length(basis_irred.kpoints)
        kpt_irred = basis_irred.kpoints[ik_irred]
        for symop in basis_irred.ksymops[ik_irred]
            Sk_irred = normalize_kpoint_coordinate(symop[1] * kpt_irred.coordinate)
            k_unfolded = normalize_kpoint_coordinate(kpt_unfolded.coordinate)
            if (Sk_irred ≈ k_unfolded) && (kpt_unfolded.spin == kpt_irred.spin)
                return ik_irred, symop
            end
        end
    end
    error("Invalid unfolding of BZ")
end

function unfold_array_(basis_irred, basis_unfolded, data, is_ψ)
    if basis_irred == basis_unfolded
        return data
    end
    if !(basis_irred.comm_kpts == basis_irred.comm_kpts == MPI.COMM_WORLD)
        error("Brillouin zone symmetry unfolding not supported with MPI yet")
    end
    data_unfolded = similar(data, length(basis_unfolded.kpoints))
    for ik_unfolded in 1:length(basis_unfolded.kpoints)
        kpt_unfolded = basis_unfolded.kpoints[ik_unfolded]
        ik_irred, symop = unfold_mapping(basis_irred, kpt_unfolded)
        if is_ψ
            # transform ψ_k from data into ψ_Sk in data_unfolded
            kunfold_coord = kpt_unfolded.coordinate
            @assert normalize_kpoint_coordinate(kunfold_coord) ≈ kunfold_coord
            _, ψSk = apply_ksymop(symop, basis_irred,
                                  basis_irred.kpoints[ik_irred], data[ik_irred])
            data_unfolded[ik_unfolded] = ψSk
        else
            # simple copy
            data_unfolded[ik_unfolded] = data[ik_irred]
        end
    end
    data_unfolded
end

function unfold_bz(scfres)
    basis_unfolded = unfold_bz(scfres.basis)
    ψ = unfold_array_(scfres.basis, basis_unfolded, scfres.ψ, true)
    eigenvalues = unfold_array_(scfres.basis, basis_unfolded, scfres.eigenvalues, false)
    occupation = unfold_array_(scfres.basis, basis_unfolded, scfres.occupation, false)
    E, ham = energy_hamiltonian(basis_unfolded, ψ, occupation;
                                scfres.ρ, eigenvalues, scfres.εF)
    @assert E.total ≈ scfres.energies.total
    new_scfres = (; basis=basis_unfolded, ψ, ham, eigenvalues, occupation)
    merge(scfres, new_scfres)
end


text:
function validate_serialization(sys::IS.SystemData; time_series_read_only = false)
    #path, io = mktemp()
    # For some reason files aren't getting deleted when written to /tmp. Using current dir.
    filename = "test_system_serialization.json"
    @info "Serializing to $filename"

    try
        if isfile(filename)
            rm(filename)
        end
        IS.prepare_for_serialization!(sys, filename; force = true)
        data = IS.serialize(sys)
        open(filename, "w") do io
            JSON3.write(io, data)
        end
    catch
        rm(filename)
        rethrow()
    end

    # Make sure the code supports the files changing directories.
    test_dir = mktempdir()
    path = mv(filename, joinpath(test_dir, filename))

    t_file = splitext(basename(path))[1] * "_" * IS.TIME_SERIES_STORAGE_FILE
    mv(t_file, joinpath(test_dir, t_file))
    v_file = splitext(basename(path))[1] * "_" * IS.VALIDATION_DESCRIPTOR_FILE
    mv(v_file, joinpath(test_dir, v_file))

    data = open(path) do io
        JSON3.read(io, Dict)
    end

    orig = pwd()
    try
        cd(dirname(path))
        sys2 = IS.deserialize(
            IS.SystemData,
            data;
            time_series_read_only = time_series_read_only,
        )
        # Deserialization of components should be directed by the parent of SystemData.
        # There isn't one in IS, so perform the deserialization in the test code.
        for component in data["components"]
            type = IS.get_type_from_serialization_data(component)
            comp = IS.deserialize(type, component)
            IS.add_component!(sys2, comp, allow_existing_time_series = true)
        end
        return sys2, IS.compare_values(sys, sys2, compare_uuids = true)
    finally
        cd(orig)
    end
end

@testset "Test JSON serialization of system data" begin
    for in_memory in (true, false)
        sys = create_system_data_shared_time_series(; time_series_in_memory = in_memory)
        _, result = validate_serialization(sys)
        @test result
    end
end

@testset "Test prepare_for_serialization" begin
    sys = create_system_data_shared_time_series()
    directory = joinpath(mktempdir(), "dir2")
    IS.prepare_for_serialization!(sys, joinpath(directory, "sys.json"))
    @test IS.get_ext(sys.internal)["serialization_directory"] == directory
end

@testset "Test JSON serialization of with read-only time series" begin
    sys = create_system_data_shared_time_series(; time_series_in_memory = false)
    sys2, result = validate_serialization(sys; time_series_read_only = true)
    @test result
end

@testset "Test JSON serialization of with mutable time series" begin
    sys = create_system_data_shared_time_series(; time_series_in_memory = false)
    sys2, result = validate_serialization(sys; time_series_read_only = false)
    @test result
end

@testset "Test verion info" begin
    data = IS.serialize_julia_info()
    @test haskey(data, "julia_version")
    @test haskey(data, "package_info")
end


text:
struct InputData
    n    ::Int
    m    ::Int
    pieces ::Matrix{Int}
end

function parseInput(raw_input; order=[1,2,3,4])
    lines = split(raw_input, '\n')

    firstLine = split(lines[1], ' ')

    n = parse(Int, firstLine[1])
    m = parse(Int, firstLine[2])

    pieces = Matrix{Int}(undef, n*m, 4)

    for i = 2:n*m+1
        line = split(lines[i],' ')
        for j =1:4
            pieces[i-1,j] = parse(Int, line[order[j]])
        end
    end

    return InputData(n, m, pieces)
end

function getInputData(filename;order=[1,2,3,4])

    inputData = nothing

    open(filename, "r") do openedFile
        inputData = parseInput(read(openedFile, String);order=order)
    end

    return inputData
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
#=

  Running all 1..50 Euler problems in Julia.

  This is a file where I've collected the best variants of the 1..50
  Euler problems. Mostly it's a test to see how fast it will be
  by not handling 50 separate Julia programs.

  Running this in Julia REPL:
  * first run:
  """
  Total time: 1.39524527
  6.209198 seconds (27.90 M allocations: 1.515 GiB, 5.93% gc time)
  """

  * second run (everything JIT:ed):
  """
  ...
  Total time: 1.0628387979999998
  3.564508 seconds (20.37 M allocations: 1.121 GiB, 5.59% gc time)
  """

  This Julia model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/Julia_progs/

=#

using Memoization
include("Euler.jl")

using BenchmarkTools
const run_benchmark = false

const res = []
const total_time = [0.0]

function run_euler(p)
    if run_benchmark
        println("\nrun_euler($p)")
        # @run_euler_bench p
        # @benchmark $p
        # b = @benchmarkable $p
        # tune!(b)
        # display(run(b))
        # b = @benchmark $p
        #dump(b)
        # Note the '$p()' construct!

        b = @benchmark $p()
        display(b)
        m = mean(b)
        println("$p mean:$m\n\n")
        push!(res, [p,m])
    else
        # GC.gc(false)
        t = @timed p()
        # GC.gc(true)
        # @printf "%s: %d: %2.8fs\n" p t.value t.time
        push!(res, [p,t.time])
        total_time[1] += t.time
    end
end

# 0.00000276s
function euler1c()
   s = 0
   for i in 1:999
       if i % 3 == 0 || i % 5 == 0
           s+=i
       end
   end
   s
end
run_euler(euler1c)

const fibmem2b = Dict{Int,Int}()
function fib2b(n::Int)
    get!(fibmem2b, n) do
        n < 3 ? 1 : fib2b(n-1) + fib2b(n-2)
    end
end

# 8.108e-6s
function euler2b()
    n = 1
    f = 0
    s = 0
    while f < 4_000_000
        f = fib2b(n)
        if f % 2 == 0
            s += f
        end
        n += 1
    end
    return s
end
run_euler(euler2b)

# 0.00169180s
function euler3c()
    return maximum(factors(600851475143))
end
run_euler(euler3c)

# 0.00958293s
function euler4a()
    max = 0
    from = 100
    to = 999
    for i in from:to, j in i:to
        ij = i*j
        if ij > max && palindromic_number(ij)
            max = ij
        end
    end
    return max;
end
run_euler(euler4a)

# 0.00000113s
function euler5a()
    a = 1
    for i in 2:20
        a = lcm(a,i)
    end
    return a
end
run_euler(euler5a)

# 0.00000085s
function euler6b()
    sum(1:100)^2 - sum((1:100).|>i->i^2)
end
run_euler(euler6b)

# 0.00981122s
function euler7b()
    p = primes(200000) # slightly cheating...
    return p[10001];
end
run_euler(euler7b)

# digits works with this
n = 7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450

# 0.04791558s
function euler8a()
    # a = split(n,"").|>i->parse(Int8,i)
    a = digits(n)
    len = length(a)
    max = 0
    for i in 1:(len-5+1)
        p = 1
        for j = 0:4
            p *= a[i+j]
            if p==0
                continue;
            end
        end
        if p > max
            max = p;
        end
    end
    return max;
end
run_euler(euler8a)

# 0.00683682s
function euler9a()
    for c in 1:500
        for b in 1:c
            for a in 1:b
                if a+b+c === 1000 && a^2 + b^2 - c^2 == 0
                    return a*b*c;
                end
            end
        end
    end
end

run_euler(euler9a)

# 0.04506652s
function euler10c()
    return sieve(2_000_000)|>sum
end
run_euler(euler10c)

m = [
     8  2 22 97 38 15  0 40  0 75  4  5  7 78 52 12 50 77 91  8;
    49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48  4 56 62  0;
    81 49 31 73 55 79 14 29 93 71 40 67 53 88 30  3 49 13 36 65;
    52 70 95 23  4 60 11 42 69 24 68 56  1 32 56 71 37  2 36 91;
    22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80;
    24 47 32 60 99  3 45  2 44 75 33 53 78 36 84 20 35 17 12 50;
    32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70;
    67 26 20 68  2 62 12 20 95 63 94 39 63  8 40 91 66 49 94 21;
    24 55 58  5 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72;
    21 36 23  9 75  0 76 44 20 45 35 14  0 61 33 97 34 31 33 95;
    78 17 53 28 22 75 31 67 15 94  3 80  4 62 16 14  9 53 56 92;
    16 39  5 42 96 35 31 47 55 58 88 24  0 17 54 24 36 29 85 57;
    86 56  0 48 35 71 89  7  5 44 44 37 44 60 21 58 51 54 17 58;
    19 80 81 68  5 94 47 69 28 73 92 13 86 52 17 77  4 89 55 40;
     4 52  8 83 97 35 99 16  7 97 57 32 16 26 26 79 33 27 98 66;
    88 36 68 87 57 62 20 72  3 46 33 67 46 55 12 32 63 93 53 69;
     4 42 16 73 38 25 39 11 24 94 72 18  8 46 29 32 40 62 76 36;
    20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74  4 36 16;
    20 73 35 29 78 31 90  1 74 31 49 71 48 86 81 16 23 57  5 54;
     1 70 54 71 83 51 54 69 16 92 33 48 61 43 52  1 89 19 67 48
]

# 0.07859598s
function euler11a()
    len = size(m)[1]
    mt = m'
    n = 4
    maxp = 0

    # Rows and columns
    for row = 1:len-1
        for col in 1:len-n
            p1 = prod(m[row,col:col+n-1])
            if p1 > maxp
                maxp = p1
            end
            p2 = prod(mt[row,col:col+n-1])
            if p2 > maxp
                maxp = p2
            end
        end
    end

    # diag_down
    for j in 1:17
        for i in 1:17
            pa = []
            for a in 0:n-1
                push!(pa,m[a+i,a+j])
            end
            p = prod(pa)
            if p > maxp
                maxp = p
            end
        end
    end

    # diag_up
    for j in 2:17
        for i = 5:20
            pa = []
            for a = 0:3
                push!(pa,m[i-a,j+a])
            end
            p = prod(pa)
            if p > maxp
                maxp = p
            end
        end
    end

    return maxp;
end

run_euler(euler11a)

# 0.13744243s
function euler12b()
    len = 0
    i = 0
    tnum = 0
    while len <= 500
        i+=1
        tnum += i
        f = factors_map(tnum)
        len = prod(values(f).|>(i->i+1))
    end
    return tnum
end

#
# This is basically factors + collect together so it's not general
# enough for placing in Euler.jl
#
function factors_map(n)
    if n == 1
        return Dict(1=>1)
    end
    m = Dict();
    while n % 2 === 0
        if !haskey(m,2)
            m[2] = 0
        end
        m[2] += 1;
        n = round(Int, n/2);
    end
    t = 3
    while n > 1 && t < ceil(Int,sqrt(n))
        while n % t == 0
            if !haskey(m,t)
                m[t] = 0
            end
            m[t] += 1
            n /= t
        end
        t += 2
    end
    if n > 1
        if !haskey(m,n)
            m[n] = 0
        end
        m[n] += 1
    end
    return m;
end

# 0.13744243s
function euler12b()
    len = 0
    i = 0
    tnum = 0
    while len <= 500
        i+=1
        tnum += i
        f = factors_map(tnum)
        len = prod(values(f).|>(i->i+1))
    end
    return tnum
end
run_euler(euler12b);

ns =[37107287533902102798797998220837590246510135740250,
     46376937677490009712648124896970078050417018260538,
     74324986199524741059474233309513058123726617309629,
     91942213363574161572522430563301811072406154908250,
     23067588207539346171171980310421047513778063246676,
     89261670696623633820136378418383684178734361726757,
     28112879812849979408065481931592621691275889832738,
     44274228917432520321923589422876796487670272189318,
     47451445736001306439091167216856844588711603153276,
     70386486105843025439939619828917593665686757934951,
     62176457141856560629502157223196586755079324193331,
     64906352462741904929101432445813822663347944758178,
     92575867718337217661963751590579239728245598838407,
     58203565325359399008402633568948830189458628227828,
     80181199384826282014278194139940567587151170094390,
     35398664372827112653829987240784473053190104293586,
     86515506006295864861532075273371959191420517255829,
     71693888707715466499115593487603532921714970056938,
     54370070576826684624621495650076471787294438377604,
     53282654108756828443191190634694037855217779295145,
     36123272525000296071075082563815656710885258350721,
     45876576172410976447339110607218265236877223636045,
     17423706905851860660448207621209813287860733969412,
     81142660418086830619328460811191061556940512689692,
     51934325451728388641918047049293215058642563049483,
     62467221648435076201727918039944693004732956340691,
     15732444386908125794514089057706229429197107928209,
     55037687525678773091862540744969844508330393682126,
     18336384825330154686196124348767681297534375946515,
     80386287592878490201521685554828717201219257766954,
     78182833757993103614740356856449095527097864797581,
     16726320100436897842553539920931837441497806860984,
     48403098129077791799088218795327364475675590848030,
     87086987551392711854517078544161852424320693150332,
     59959406895756536782107074926966537676326235447210,
     69793950679652694742597709739166693763042633987085,
     41052684708299085211399427365734116182760315001271,
     65378607361501080857009149939512557028198746004375,
     35829035317434717326932123578154982629742552737307,
     94953759765105305946966067683156574377167401875275,
     88902802571733229619176668713819931811048770190271,
     25267680276078003013678680992525463401061632866526,
     36270218540497705585629946580636237993140746255962,
     24074486908231174977792365466257246923322810917141,
     91430288197103288597806669760892938638285025333403,
     34413065578016127815921815005561868836468420090470,
     23053081172816430487623791969842487255036638784583,
     11487696932154902810424020138335124462181441773470,
     63783299490636259666498587618221225225512486764533,
     67720186971698544312419572409913959008952310058822,
     95548255300263520781532296796249481641953868218774,
     76085327132285723110424803456124867697064507995236,
     37774242535411291684276865538926205024910326572967,
     23701913275725675285653248258265463092207058596522,
     29798860272258331913126375147341994889534765745501,
     18495701454879288984856827726077713721403798879715,
     38298203783031473527721580348144513491373226651381,
     34829543829199918180278916522431027392251122869539,
     40957953066405232632538044100059654939159879593635,
     29746152185502371307642255121183693803580388584903,
     41698116222072977186158236678424689157993532961922,
     62467957194401269043877107275048102390895523597457,
     23189706772547915061505504953922979530901129967519,
     86188088225875314529584099251203829009407770775672,
     11306739708304724483816533873502340845647058077308,
     82959174767140363198008187129011875491310547126581,
     97623331044818386269515456334926366572897563400500,
     42846280183517070527831839425882145521227251250327,
     55121603546981200581762165212827652751691296897789,
     32238195734329339946437501907836945765883352399886,
     75506164965184775180738168837861091527357929701337,
     62177842752192623401942399639168044983993173312731,
     32924185707147349566916674687634660915035914677504,
     99518671430235219628894890102423325116913619626622,
     73267460800591547471830798392868535206946944540724,
     76841822524674417161514036427982273348055556214818,
     97142617910342598647204516893989422179826088076852,
     87783646182799346313767754307809363333018982642090,
     10848802521674670883215120185883543223812876952786,
     71329612474782464538636993009049310363619763878039,
     62184073572399794223406235393808339651327408011116,
     66627891981488087797941876876144230030984490851411,
     60661826293682836764744779239180335110989069790714,
     85786944089552990653640447425576083659976645795096,
     66024396409905389607120198219976047599490197230297,
     64913982680032973156037120041377903785566085089252,
     16730939319872750275468906903707539413042652315011,
     94809377245048795150954100921645863754710598436791,
     78639167021187492431995700641917969777599028300699,
     15368713711936614952811305876380278410754449733078,
     40789923115535562561142322423255033685442488917353,
     44889911501440648020369068063960672322193204149535,
     41503128880339536053299340368006977710650566631954,
     81234880673210146739058568557934581403627822703280,
     82616570773948327592232845941706525094512325230608,
     22918802058777319719839450180888072429661980811197,
     77158542502016545090413245809786882778948721859617,
     72107838435069186155435662884062257473692284509516,
     20849603980134001723930671666823555245252804609722,
     53503534226472524250874054075591789781264330331690
]

# 0.01485372s
function euler13a()
    return parse(Int,string(sum(ns))[1:10])
end

run_euler(euler13a)

function hailstone(n)
    if n % 2 == 0
        return round(Int,n / 2)
    else
        return 3*n+1;
    end
end

# Using an array instead:
# 0.04143852s
function euler14c()
    limit = 1_000_000
    hash = zeros(limit)
    maxN = 0
    maxLen = 0
    for n = 2:1_000_000
        m = n
        mlen = 1
        while m > 1
            if m <= limit && hash[m] > 0
                mlen = hash[m]+mlen-1
                m = 1
            else
                m = hailstone(m)
                mlen += 1
            end
        end
        if hash[n] == 0
            hash[n] = mlen
        end
        if mlen > maxLen
            maxN = n
            maxLen = mlen
        end
    end
    return maxN;
end

run_euler(euler14c)

#=
function prodlist(from,to)
   return prod(from:to)
end
# 0.01229874s
function euler15a()
     return prodlist(BigInt(21),40) / prodlist(2,20);
end
run_euler(euler15a)
=#

function euler15b()
    return prod(BigInt(21):40) / prod(2:20);
  end
run_euler(euler15b)

# 0.00027433s
function euler16a()
    # return split(string(BigInt(2)^1000),"").|>(i->parse(Int,i))|>sum
    return digits(BigInt(2)^1000)|>sum
end

run_euler(euler16a);

function english(n)
    divs  =      [1_000_000_000, 1_000_000,  1_000,       100]
    divnames  =  ["billion", "million", "thousand", "hundred"]
    prefixes  =  ["0", "twen", "thir", "for", "fif", "six", "seven", "eigh", "nine"]
    _ordinals  = ["first", "second", "third", "fourth", "fifth", "sixth", "seventh",
                      "eighth", "ninth", "tenth", "eleventh", "twelfth", "thirteenth",
                      "fourteenth","fifteenth", "sixteenth", "seventeenth",
                      "eighteenth", "nineteenth"]
    cardinals =  ["one", "two", "three", "four", "five", "six", "seven",
                      "eight", "nine", "ten", "eleven", "twelve", "thirteen", "fourteen",
                      "fifteen", "sixteen", "seventeen", "eighteen", "nineteen"];

    s = ""
    printed = 0
    if n < 0
        s = "minus" * s;
        n = -n
    end
    d = 0
    for i in 1:length(divs)
        d = floor(Int,n / divs[i])
        n %= divs[i]
        if d != 0
            s *= english(d) * divnames[i];
            printed = 1
        end
    end

    if n > 0 && printed == 1
        s *= "and";
    end
    if n == 0
        # dummy
    elseif n > 19
        d = floor(Int,n / 10)
        n %= 10
        s *= prefixes[d] * "ty" * english(n)
    else
        s *= cardinals[n]
    end

    return s
end

# 0.00138457s
function euler17a()
    return sum((1:1000).|>i->length(english(i)))
end

run_euler(euler17a)

# Recursive function for euler18()
function pp(row, col, sum, tri, max_val)
    if sum > max_val
        max_val = sum
    end
    row += 1
    if row <= length(tri)
        for i in 0:1
            max_val = pp(row, col + i, sum + tri[row][col + i], tri, max_val)
        end
    end
    return max_val
end

# 0.00015198s
function euler18a()
    triangle  =
         [[75],
         [95,64],
         [17,47,82],
         [18,35,87,10],
         [20, 4,82,47,65],
         [19, 1,23,75, 3,34],
         [88, 2,77,73, 7,63,67],
         [99,65, 4,28, 6,16,70,92],
         [41,41,26,56,83,40,80,70,33],
         [41,48,72,33,47,32,37,16,94,29],
         [53,71,44,65,25,43,91,52,97,51,14],
         [70,11,33,28,77,73,17,78,39,68,17,57],
         [91,71,52,38,17,14,91,43,58,50,27,29,48],
         [63,66, 4,68,89,53,67,30,73,16,69,87,40,31],
         [ 4,62,98,27,23, 9,70,98,73,93,38,53,60, 4,23]];
    max_val = 0
    max_val = pp(1, 1, triangle[1][1], triangle, max_val)

    return max_val
end

run_euler(euler18a)

using Dates

# 0.00001321s
function euler19a()
    count = 0;
    for year in 1901:2000
        for month in 1:12
            if Dates.dayofweek(Date(year, month, 1)) == 7
                count += 1
            end
        end
    end

    return count
end

run_euler(euler19a)

# 0.00029684s
function euler20a()
    return sum(split(factorial(BigInt(100))|>string,"").|>j->parse(Int,j))
end

run_euler(euler20a)

# 0.07385759s
function euler21c()
    n = 9999
    s = zeros(Int,n)
    for i in 2:n
        s[i] = sum(all_divisors3(i))
    end
    a = []
    for i in 2:n
        # ignore perfect numbers...
        # if s[i] <= n && i==s[s[i]] && i!= s[i]
        if s[i] <= n && i!= s[i] && i==s[s[i]]
            push!(a,i)
        end
    end
    return sum(a)
end

run_euler(euler21c);

function to_code(s)
    return [Int(s[i])-64 for i in 1:length(s)]
end

# 0.00231866s
function euler22a()
    words = sort(split(readline("euler22_names.txt"),","))
    a = 1;
    s = 0;
    for word in words
        word = replace(word,"\""=>"")
        s += a*sum(to_code(word));
        a += 1;
    end
    return s;

end

# 0.02438064s
function euler22b()
    a = 1
    return sum(map(word->(word=replace(word,"\""=>"");a=a+1; (a-1)*sum(to_code(word))
                   ),
                  sort(split(readline("euler22_names.txt"),","))))

end

run_euler(euler22a)

# 0.00790154s
function euler23a()
    limit = 20161;
    arr = ones(limit)
    for i in 2:limit
        for j in i*2:i:limit
            arr[j] = arr[j] + i
        end
    end

    abundant = [i for i in 12:limit if arr[i] > i]
    for a in abundant
        for b in abundant
            if b > a || a + b >= limit
                break;
            else
                arr[a + b] = 0
            end
        end
    end

    return sum([i for i in 1:limit if arr[i] != 0])
end

run_euler(euler23a)

#
function euler24a()
    a = Array(0:9)
    c = 1
    while c < 1000000
        a = next_permutation(a)
        c += 1
    end

    return parse(Int, join(a,""))
end
run_euler(euler24a)

function fib_len(n)
    return fibmemBigInt(n)|>string|>length
end

#
# Using some heuristics to find the upper limit
# (from my Picat solution).
# 0.08765257s
function euler25b()
    target = 1000;
    foundUpper = 0;
    i = 1;
    fibLen = 0;
    step = 43;
    # Get the upper limit
    while fibLen < target && foundUpper == 0
        fibLen = fib_len(step*i);
        if fibLen > target
            foundUpper = i
            break
        end
        i += 1
    end

    # Now check all numbers from Step*(FoundUpper-1) .. Step*FoundUpper
    # The target must be in that interval.
    f = step*(foundUpper-1)
    fibLen = fib_len(f);
    while fibLen < target && f <= step*foundUpper
        f += 1
        fibLen = fib_len(f)
    end

    return f
end

run_euler(euler25b)

function get_rep_len(n)
    foundRemainders = zeros(n)
    value = 1;
    position = 1;
    while foundRemainders[value] == 0 && value != 0
        foundRemainders[value] = position
        value = round(Int,1+((value*10) % n)) # added 1
        position += 1
    end
    return position-foundRemainders[value]
end

# Checks only primes
# 0.00048207s
function euler26b()
    maxLen = 0
    maxD = 0
    for d in 2:1000
        if isPrime(d)
            len = get_rep_len(d)
            if len > maxLen
                maxLen = len
                maxD = d
            end
        end
    end
    return maxD
end

run_euler(euler26b)

# 0.04097529s
# with isPrimeCached: 0.07010393s (slower)
function euler27a()
    t = 999
    bestLen = 0
    bestA = 0
    bestB = 0
    for a in -t:t
        for b in -t:t
            len = 0;
            pp = len^2 + a*len + b
            # while pp > 1 && isPrimeCached(pp) # slower
            while pp > 1 && isPrime(pp)
                len += 1
                pp = len^2 + a*len + b;
            end
            if len > bestLen
                bestLen = len
                bestA = a
                bestB = b
            end
        end
    end

    return(bestA * bestB)
end

run_euler(euler27a)

#
function euler28a()
    s = 1
    n = 3
    while n <= 1001
        s += 4*n^2 - 6*n+6
        n+=2
    end
    return s
end

run_euler(euler28a)

# 0.00789801s
function euler29a()
    min = 2
    max = 100
    h = Dict()
    for a = min:max
        for b = min:max
            h[BigInt(a)^BigInt(b)] = 1
        end
    end

    return keys(h)|>length
end
run_euler(euler29a)

# 0.07169755ss
function euler30a()
    t = 0
    m = 5
    for n in 10:6 * (9^5)
        # nn = (split(string(n),"").|>i->parse(Int,i)^m)|>sum
        nn = (digits(n).|>i->i^m)|>sum
        if n === nn
            t += n
        end
    end
    return t
end

run_euler(euler30a)

# DP approach
function coins(c,money,m)
    len = length(c)
    if m == len
        return 1
    end

    sum1 = 0
    for i in m:len
        if money - c[i] == 0
            sum1 += 1
        end
        if  money - c[i] > 0
            sum1 += coins(c, money-c[i], i)
        end
    end
    return sum1;
end

# 0.00058022s
function euler31a()
    c = [200,100,50,20,10,5,2,1]
    return coins(c, 200, 1)
end

run_euler(euler31a)

# 0.17556319s
function euler32a()
    prod_hash = Dict{Int,Int}()
    for a in 2:98
        as = string(a)
        if occursin("0",as) continue end
        for b in a+1:9876
            p = a*b
            if occursin("0",string(p)) continue end
            l = "$as$b$p"
            if length(l) == 9 && length(Set(l)) == 9 && !occursin("0",l)
                prod_hash[p] = 1
            end
        end
    end
    return keys(prod_hash)|>sum
end

GC.gc(true)
GC.gc(false)
run_euler(euler32a)
GC.gc(true)

# 0.00000079s
function euler33a()
    s = 1;
    for y in 1:9
        for z in 1:9
            x = 9.0*y*z/(10.0*y-z)
            if floor(x)==x && y/z < 1.0 && x < 10.0
                s = (s*y)/z;
            end
        end
    end

    return 1/s;
end

run_euler(euler33a)

function factorial_sum(n)
    # return (split(string(n),"").|>i->factorial(parse(Int,i)))|>sum
    return (digits(n).|>i->factorial(i))|>sum

end

# 0.01983716s
function euler34c()
    return [n for n in 10:100000 if n == factorial_sum(n)]|>sum
end

run_euler(euler34c)

# Rotate an array/string
function rotate(n,i)
    s = string(n)
    ii = i+1
    return parse(Int,join(vcat(s[ii:end],s[1:ii-1])))
end

#
# Note: It's a little faster when using prime_set as
#       a parameter for is_circular_prime than
#       relying of the global prime_set in the function.
#       That's surprising...
#
function is_circular_prime(n, prime_set)
    # s = split(string(n),"") # num_to_list(n)
    # len = length(s)
    len = nlen(n)
    v = 0
    for i = 1:len
        v = rotate(n,i)
        if !(v in prime_set)
            return false;
        end
    end
    return v in prime_set
end

# slightly slower than is_circular_prime
function is_circular_prime2(n, prime_set)
    s = split(string(n),"") # num_to_list(n)
    len = length(s)
    v = n
    for i = 2:len
        v = parse(Int,join(vcat([s[j] for j in i:len], [s[j] for j in 1:i-1]),""))
        if !(v in prime_set)
            return false;
        end
    end
    return v in prime_set
end

# Faster if global
prime_set = Set(primes(1_000_000))

# 0.08751774s
function euler35a()
    # prime_set = Set(primes(1_000_000))
    numCircularPrimes = 0
    for n in prime_set
        if is_circular_prime(n,prime_set)
            numCircularPrimes += 1
        end
    end

    return numCircularPrimes
end

run_euler(euler35a)

# 0.06355465s
function euler36a()
    s = 0
    for n in 1:999999
        if palindromic_number(n) &&
            palindromic_list(digits(n,base=2))
            # palindromic_list(string(n,base=2))
            s += n;
        end
    end
    return s
end

run_euler(euler36a)

# 0.00016527s
function euler38a()
    for n in 9876:-1:9
        println
        s = string(n)
        i = 2
        while length(s) < 9
            s *= string(n*i)
            i += 1
        end
        if length(s) == 9 && is_pandigital(s)
            return parse(Int,s)
        end
    end
    return nothing
end

run_euler(euler38a)

# using argmax for finding the max value
# 0.05935469s
function euler39b()
    n = 1000;
    # squares = Set((1:n).|>i->i*i) # slightly slower
    squares = Dict((1:n).|>i->i*i=>1)
    valid = []
    for x in keys(squares)
        for y in keys(squares)
            if x < y && haskey(squares,x+y) && (sqrt(x) + sqrt(x) + sqrt(x+y)) < 1000
                push!(valid,[x, y])
            end
        end
    end

    counts = Dict()
    for x in valid
        c = floor(Int,sqrt(x[1]) + sqrt(x[2]) + sqrt(x[1]+x[2]))
        get!(counts,c,0)
        counts[c] += 1
   end

    # find max count
    return argmax(counts)
end

run_euler(euler39b)

# 0.00941826s
function euler40a()
    i = 1
    dlen = 1
    p = 1
    index = 10; # Index = 10, 100, 1000, ..., 1000000
    while dlen <= 1000000
        i += 1
        istr = string(i)
        istrlen = length(istr)
        if dlen+istrlen>=index
            p *= parse(Int,istr[index-dlen])
            index *= 10
        end
        dlen += istrlen
    end

    return p
end

run_euler(euler40a)

# 0.03632750s
#=
function euler41a()
    # Simplification:
    # n=9 is not possible since 1+2+3+4+5+6+7+8+9=45 is divisible by 3
    #  n=8 is not possible since 1+2+3+4+5+6+7+8=36 is divisible by 3
    n = 7
    while n >= 4
        p = Array(1:n)
        for pp in sort(all_permutations(p),rev=true)
            v = parse(Int,join(pp,""))
            if isPrime(v)
                return v
            end
        end
        n += 1
    end

    return nothing
end

run_euler(euler41a)
=#
using Combinatorics
function euler41b()
    # Simplification:
    # n=9 is not possible since 1+2+3+4+5+6+7+8+9=45 is divisible by 3
    #  n=8 is not possible since 1+2+3+4+5+6+7+8=36 is divisible by 3
    n = 7
    m = 0
    while m == 0 && n >= 4
        p = Array(1:n)
        # for pp in sort(all_permutations(p),rev=true)
        for pp in sort(collect(permutations(p)),rev=true)
            v = parse(Int,join(pp,""))
            if isPrime(v)
                return v
            end
        end
        n += 1
    end

    return nothing
end
run_euler(euler41b)

# n'th triangle number
function triangle_number(n)
    return floor(Int,(n*(n-1)) / 2)
end

# get the score for a name
function get_score(name)
    return ([name[i] for i in 1:length(name)]
           .|>c->(Int(c)-64))|>sum
end

# 0.09908583s
function euler42a()
    t20 = Set((1:20).|>i->triangle_number(i))
    words = split( (readline("euler42_words.txt").|>word->replace(word,"\""=>"")),",")
    return (words
           .|>word-> get_score(word) in t20
           )|>sum
end

run_euler(euler42a)

# Using next_permutation
# 0.05306371s
function euler43b()
    primes = [2,3,5,7,11,13,17]
    p = [1,0,2,3,4,5,6,7,8,9]
    sum = 0
    rev = Array(9:-1:0)
     while p != rev
        i = 0
        found = true
        while i < 7 && found == true
            if (100*p[i+2] + 10*p[i+3] + p[i+4]) % primes[i+1] != 0
                found = false
                break
            end
            i += 1
        end

        if found
            sum += parse(Int,join(p,""))
        end
        p = next_permutation(p)
    end

    return sum
end

run_euler(euler43b)

# 0.04486643s
function euler44b()
    n = 2500
    s = (1:2500).|>n->floor(Int,n*(3*n-1) / 2)
    t = Set(s)
    d = 10_000_000
    for j in 1:n
        for k = 1:n
            if s[j] >= s[k]
                continue;
            end
            a = s[j]+s[k]
            if a >= d
                continue
            end
            b = abs(s[j]-s[k])
            if b >= d
                continue
            end
            if a in t && b in t
                d = b
                break

            end
        end
    end
    return d
end

run_euler(euler44b)

pent(n) = floor(Int, n*(3*n-1) / 2)
tri(n) = floor(Int, n*(n+1) / 2)
hex(n) = n*(2*n-1)

# 0.00013708s
function euler45a()
    t = 285+1
    tt = tri(t)
    p = 165
    pp = pent(p)
    h = 143
    hh = hex(h)
    while tt != pp || pp != hh
        t += 1
        tt = tri(t)
        if tt > pp p+=1; pp = pent(p) end
        if pp > hh h+=1; hh = hex(h) end
        if tt > hh h+=1; hh = hex(h) end
    end

    return tt
end

run_euler(euler45a)

# 0.00057412s
function euler46a()
    res = 0
    gotit = false
    for i in 3:2:10000
        if !isPrime(i) && !gotit
            s = round(Int,sqrt(i/2))
            found = 0
            for j in 1:s
                if found === 0
                    ts = j*j*2
                    if isPrime(abs(i-ts))
                        found = 1
                    end
                end
            end
            if found == 0
                res = i
                gotit = true
                break
            end
        end
    end
    return res
end

run_euler(euler46a)

# 0.01653001s
function euler47a()
    maxn = 1_000_000
    f = zeros(maxn)
    for i in 2:maxn-1
        if f[i] == 0
            for j in 2*i:i:maxn-1
                f[j] += 1
            end
        end
    end

    goal = [4,4,4,4]
    found = 0
    for i in 1:maxn-3
        if [f[i],f[i+1],f[i+2],f[i+3]] == goal
            found = i
            break
        end
    end
    return found
end

run_euler(euler47a)

# 0.00162826s
function euler48a()
    sum = 0
    t = 10_000_000_000
    for i = 1:1_000
        n = i
        for j in 2:i
            n = (n * i) % t
        end
        sum = (sum + n) % t
    end

    return sum
end

run_euler(euler48a)

function check_perms(n, diff)
    allperms = all_permutations2(digits(n))
    if length(allperms) > 0
        p1 = get_element(n, allperms, diff)
        if p1 !== nothing
            p2 = get_element(p1, allperms, diff)
            if p2 !== nothing
                return [n, p1, p2]
            end
        end
    end

    return nothing
end

function get_element(n, ll, diff)
    for p in ll
        pp = parse(Int,join(p,""))
        # if isPrime(pp) && pp > n && pp-n == diff
        if pp > n && pp-n == diff && isPrime(pp)
            return pp
        end
    end
    return nothing
end

# 0.10001170s
function euler49a()
    diff = 3330
    res = 0
    for n in 1001:2:9999
        if n != 1487 && isPrime(n)
            c = check_perms(n, diff)
            if c !== nothing
                res = c
                break
            end
        end
    end

    return parse(Int,join(res,""))

end

run_euler(euler49a)

# 0.00700822s
function euler50a()
    n = 10000
    p = primes(n)
    for len in 550:-1:21
        for offset in 1:549
            pp = ((offset+1:offset+len).|>i->p[i])|>sum
            if pp < 1000000 && isPrime(pp)
                return pp
            end
        end
    end

    return nothing
end

run_euler(euler50a)

println("\nResult in reverse time order:")
for r in sort(res,by=x->x[2],rev=true)
  println("$(r[1]): $(r[2])")
end

println("Total time: $(total_time[1])")

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using Documenter, DocStringExtensions, DynamicEnergyBudgets, Microclimate, Photosynthesis

makedocs(
    modules = [DynamicEnergyBudgets],
    sitename = "DynamicEnergyBudgets.jl",
    pages = Any[
        "Home" => "index.md",
    ],
    clean = false,
)

deploydocs(
    repo = "github.com/rafaqz/DynamicEnergyBudgets.jl.git",
)


text:
nSamples = 5 # The number of samples you have data for
smpl = NewChronAgeData(nSamples)
smpl.Name      =  ("KJ08-157", "KJ04-75", "KJ09-66",  "KJ04-72", "KJ04-70",)
smpl.Height   .=  [     -52.0,      44.0,      54.0,      82.0,      93.0,]
smpl.Height_sigma .= [    3.0,       1.0,       3.0,       3.0,       3.0,]
smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
smpl.Path = "../examples/DenverUPbExampleData/" # Where are the data files?
smpl.inputSigmaLevel = 2 # i.e., are the data files 1-sigma or 2-sigma. Integer.
smpl.Age_Unit = "Ma" # Unit of measurement for ages and errors in the data files
smpl.Height_Unit = "cm" # Unit of measurement for Height and Height_sigma

# Remove outliers (if any)
smpl = screen_outliers(smpl, maxgap=50, make_plots=false)

# Distribution boostrapping from chron strat object
BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl)
@test isa(BootstrappedDistribution, Array{Float64,1})

# Estimate the eruption age distributions for each sample  - - - - - - - -

# Configure distribution model here
distSteps = 10^5 # Number of steps to run in distribution MCMC
distBurnin = floor(Int,distSteps/2) # Number to discard

# Run MCMC to estimate saturation and eruption/deposition age distributions
@time tMinDistMetropolis(smpl,distSteps,distBurnin,BootstrappedDistribution; make_plots=false)

# Run stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - - -
# Configure the stratigraphic Monte Carlo model
config = NewStratAgeModelConfiguration()
# If you in doubt, you can probably leave these parameters as-is
config.resolution = 10.0 # Same units as sample height. Smaller is slower!
config.bounding = 0.5 # how far away do we place runaway bounds, as a fraction of total section height
(bottom, top) = extrema(smpl.Height)
npoints_approx = round(Int,length(bottom:config.resolution:top) * (1 + 2*config.bounding))
config.nsteps = 100000 # Number of steps to run in distribution MCMC
config.burnin = 100000*npoints_approx # Number to discard
config.sieve = round(Int,npoints_approx) # Record one out of every nsieve steps
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

# Run the stratigraphic MCMC model
@time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config)

# Test that results match expectation, within some tolerance
@test isa(mdl.Age, Array{Float64,1})
@test isapprox(mdl.Age, [66.06, 66.05, 66.03, 66.02, 66.01, 66.0, 65.98, 65.97, 65.96, 65.94, 65.94, 65.93, 65.93, 65.9], atol=0.1)
@test isapprox(mdl.Age_025CI, [66.0, 65.98, 65.96, 65.95, 65.94, 65.93, 65.93, 65.92, 65.92, 65.91, 65.91, 65.89, 65.88, 65.82], atol=0.15)
@test isapprox(mdl.Age_975CI, [66.09, 66.09, 66.08, 66.08, 66.07, 66.07, 66.05, 66.04, 66.02, 65.97, 65.97, 65.96, 65.96, 65.95], atol=0.15)

# Data about hiatuses
nHiatuses = 2 # The number of hiatuses you have data for
hiatus = NewHiatusData(nHiatuses) # Struct to hold data
hiatus.Height         = [-7.0, 35.0 ]
hiatus.Height_sigma   = [ 0.0,  0.0 ]
hiatus.Duration       = [ 0.3,  0.3 ]
hiatus.Duration_sigma = [ 0.05, 0.05]

# Run the model. Note the additional `hiatus` arguments
@time (mdl, agedist, hiatusdist, lldist) = StratMetropolisDist(smpl, hiatus, config)

# Test that results match expectation, within some tolerance
@test isa(mdl.Age, Array{Float64,1})
@test isapprox(mdl.Age, [66.08, 66.07, 66.07, 66.07, 66.02, 66.01, 66.01, 66.01, 65.94, 65.94, 65.93, 65.93, 65.92, 65.9], atol=0.1)
@test isapprox(mdl.Age_025CI, [66.05, 66.04, 66.03, 66.02, 65.94, 65.94, 65.93, 65.93, 65.91, 65.9, 65.9, 65.89, 65.88, 65.82], atol=0.15)
@test isapprox(mdl.Age_975CI, [66.1, 66.1, 66.1, 66.1, 66.08, 66.08, 66.08, 66.08, 65.98, 65.96, 65.96, 65.96, 65.96, 65.95], atol=0.15)

## --- As above, but treat everything as a gaussian/weighted mean

# Tel tMinDistMetropolis to treat these as gaussians, using the first row of data file
smpl.Age_DistType.=1

# Run MCMC to estimate saturation and eruption/deposition age distributions
@time tMinDistMetropolis(smpl,distSteps,distBurnin,BootstrappedDistribution; make_plots=false)

# Run the stratigraphic MCMC model
@time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config)
@test isa(mdl.Age, Array{Float64,1})
@test isapprox(mdl.Age, [65.97, 65.97, 65.96, 65.95, 65.95, 65.94, 65.93, 65.92, 65.92, 65.91, 65.9, 65.89, 65.87, 65.85], atol=0.1)
@test isapprox(mdl.Age_025CI, [65.86, 65.86, 65.85, 65.85, 65.84, 65.84, 65.84, 65.83, 65.83, 65.83, 65.82, 65.77, 65.74, 65.72], atol=0.15)
@test isapprox(mdl.Age_975CI, [66.07, 66.06, 66.06, 66.06, 66.06, 66.05, 66.04, 66.03, 66.02, 66.0, 65.99, 65.98, 65.97, 65.96], atol=0.15)


text:
# create_link
struct Create_LinkParms
    clientId::Int32
    lockDevice::Bool
    lock_timeout::UInt32
    device::String
end

const Device_Link = UInt32
const Device_ErrorCode = UInt32
const Device_Flags = UInt32

struct Create_LinkResp
    error::Device_ErrorCode
    link::Device_Link
    abortPort::UInt16
    maxRecvSize::UInt32
end

const create_link = make_rpc(395183, 1, 10, Create_LinkResp, Create_LinkParms)

# destroy_link
const destroy_link = make_rpc(395183, 1, 23, Device_ErrorCode, Device_Link)

# device_write
struct Device_WriteParms
    lid::Device_Link
    io_timeout::UInt32
    lock_timeout::UInt32
    flags::Device_Flags
    data::Vector{UInt8}
end

struct Device_WriteResp
    error::Device_ErrorCode
    size::UInt32
end

const device_write = make_rpc(395183, 1, 11, Device_WriteResp, Device_WriteParms)

# device_read
struct Device_ReadParms
    lid::Device_Link
    request_size::UInt32
    io_timeout::UInt32
    lock_timeout::UInt32
    flags::Device_Flags
    termChar::UInt8
end

struct Device_ReadResp
    error::Device_ErrorCode
    reason::UInt32
    data::Vector{UInt8}
end

const device_read = make_rpc(395183, 1, 12, Device_ReadResp, Device_ReadParms)


text:
"""
maximum_weight_matching(g::Graph, w::Dict{Edge,Real} -> Dict{Edge,Int64}

Given a graph `g` and an edgemap `w` containing weights associated to edges,
returns a matching with the maximum total weight.
`w` is a dictionary that maps edges i => j to weights.
If no weight parameter is given, all edges will be considered to have weight 1
(results in max cardinality matching)

The efficiency of the algorithm depends on the input graph:
  - If the graph is bipartite, then the LP relaxation is integral.
  - If the graph is not bipartite, then it requires a MIP solver and
  the computation time may grow exponentially.

The package JuMP.jl and one of its supported solvers is required.

Returns MatchingResult containing:
  - a solve status (indicating whether the problem was solved to optimality)
  - the optimal cost
  - a list of each vertex's match (or -1 for unmatched vertices)
"""
function maximum_weight_matching end

function maximum_weight_matching(g::Graph,
          solver::JuMP.OptimizerFactory,
          w::AbstractMatrix{U} = default_weights(g)) where {U <:Real}

    model = Model(with_optimizer(solver))
    n = nv(g)
    edge_list = collect(edges(g))

    # put the edge weights in w in the right order to be compatible with edge_list
    for j in 1:n
      for i in 1:n
        if i > j  && w[i,j] > zero(U) && w[j,i] < w[i,j]
          w[j,i] = w[i,j]
        end
        if Edge(i,j) ∉ edge_list
          w[i,j] = zero(U)
        end
      end
    end

    if is_bipartite(g)
      @variable(model, x[edge_list] >= 0) # no need to enforce integrality
    else
      @variable(model, x[edge_list] >= 0, Int) # requires MIP solver
    end
    @objective(model, Max, sum(x[e]*w[src(e),dst(e)] for e in edge_list))

    @constraint(model, c1[i=1:n], sum(x[Edge(minmax(i,j))] for j in neighbors(g,i)) <= 1)
    optimize!(model)
    status = JuMP.termination_status(model)
    status != MOI.OPTIMAL && error("JuMP solver failed to find optimal solution.")
    solution = value.(x)
    cost = objective_value(model)
    return MatchingResult(cost, dict_to_arr(n, solution, edge_list))
end

""" Returns an array of mates from a dictionary that maps edges to {0,1} """
function dict_to_arr(n::Int64, solution::JuMP.Containers.DenseAxisArray{U,1,Tuple{Array{E,1}}}, edge_list::AbstractVector{E}) where {U<: Real, E<: Edge}
  mate = fill(-1,n)
  for e in edge_list
    if solution[e] >= 1 - 1e-5 # Some tolerance to numerical approximations by the solver.
        mate[src(e)] = dst(e)
        mate[dst(e)] = src(e)
    end
  end
  return mate
end

function default_weights(g::G) where {G<:AbstractGraph}
  m = spzeros(nv(g),nv(g))
  for e in edges(g)
    m[src(e),dst(e)] = 1
  end
  return m
end


text:
global const check_syms = Tuple{Symbol,Symbol}[]

isexpr(x::Expr) = true
isexpr(x::Any) = false
isexpr(x::Expr, s::Symbol) = x.head == s
isexpr(x::Any, s::Symbol) = false

istype(x::Expr) = x.head == :curly
istype(::Symbol) = true

macro fun(expr)
  @assert isexpr(expr)
  if expr.head == :(=)
    lhs,rhs = expr.args
    if isexpr(lhs, :(:))
      error("found : where :: was expected")
    end
    @assert isexpr(lhs, :(::))
    call,rettype = lhs.args
    @assert istype(rettype)
    @assert isexpr(rhs, :tuple)
    rhargs = rhs.args
  elseif expr.head == :(::)
    call,rettype = expr.args
    @assert istype(rettype)
    rhargs = call.args[2:end]
  else
    error("Not a @fun expression ($(expr.head))")
  end

  @assert isexpr(call, :call)
  fname = call.args[1]
  jargs = call.args[2:end]
  jdict = Dict{Symbol,Expr}()
  for i in jargs
    ta = isexpr(i, :kw) ? i.args[1] : i
    if isexpr(ta, :(:))
      error("found : where :: was expected")
    end
    @assert isexpr(ta, :(::))
    @assert isa(ta.args[1], Symbol)
    jdict[ta.args[1]] = ta
  end

  a = quote
    #export $(esc(fname))
    $(esc(call)) = ccall(Libdl.dlsym($lib, $(string(fun_prefix, fname))), $rettype, ())
  end

  b = a.args[2]
  @assert isexpr(b, :(=))
  cc = b.args[2]
  @assert isexpr(cc, :ccall)

  ccall_params = cc.args
  @assert isexpr(ccall_params[3], :tuple)
  ccall_types = ccall_params[3].args
  for i in rhargs
    ta2 = isa(i, Symbol) ? jdict[i] : i
    if isexpr(ta2, :(:))
      error("found : where :: was expected")
    end
    @assert isexpr(ta2, :(::))
    push!(ccall_types, esc(ta2.args[2]))
    push!(ccall_params, esc(ta2.args[1]))
  end

  push!(check_syms, (lib, symbol(fun_prefix, fname)))

  return a
end


text:
@testset "Quadratic and hermitian forms" begin
  include("QuadForm/Basic.jl")
  include("QuadForm/Lattices.jl")
  include("QuadForm/Genus.jl")
  include("QuadForm/Morphism.jl")
  include("QuadForm/Enumeration.jl")
  include("QuadForm/MassQuad.jl")
  include("QuadForm/MassHerm.jl")
  include("QuadForm/Quad.jl")
  include("QuadForm/QuadBin.jl")
  include("QuadForm/Herm.jl")
  include("QuadForm/Torsion.jl")
  include("QuadForm/PadicLift.jl")
end


text:
using Printf

function _menu(items)
    for (ind, item) in enumerate(items)
        @printf "  %2i) %s\n" ind item
    end
end

_ok(::Any,::Any) = false
function _ok(reply::AbstractString, itemcount)
    n = tryparse(Int, reply)
    return isnull(n) || 0 ≤ get(n) ≤ itemcount
end

"Prompt to select an item from the items"
function _selector(items, prompt::AbstractString)
    isempty(items) && return ""
    reply = -1
    itemcount = length(items)
    while !_ok(reply, itemcount)
        _menu(items)
        print(prompt)
        reply = strip(readline(STDIN))
    end
    return items[parse(Int, reply)]
end

items = ["fee fie", "huff and puff", "mirror mirror", "tick tock"]
item = _selector(items, "Which is from the three pigs: ")
println("You chose: ", item)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using AdversarialDriving
using AutomotiveSimulator
using AutomotiveVisualization
using Random
using Test
using Distributions

# Check out roadways
# a = render([Tint_roadway])
# a = render([ped_roadway, crosswalk])

## Test the struct types and construction
pos_noise = Noise((0.3, 0.1), 2.)
neg_noise = Noise((-0.3, -0.1), -2.)
@test pos_noise.pos[1] == 0.3
@test pos_noise.pos[2] == 0.1
@test pos_noise.vel == 2.

ped = NoisyPedState(VehicleState(VecSE2(25, -10, π/2), ped_roadway, 3.), pos_noise)
ped_entity = Entity(ped, PEDESTRIAN_DEF, 1)
@test ped isa NoisyState
@test posg(ped) == VecSE2(25, -10, π/2)
@test vel(ped) == 3.
@test noise(ped) == pos_noise
@test laneid(ped_entity) == 2

veh = BlinkerState(VehicleState(VecSE2(50., 0., 0.), ped_roadway, 3.), true, [3,4], neg_noise)
bv_entity = Entity(veh, VehicleDef(), 1)
@test veh isa NoisyState
@test posg(veh) == VecSE2(50., 0, 0.)
@test vel(veh) == 3.
@test blinker(veh)
@test goals(veh) == [3,4]
@test noise(veh) == neg_noise
@test laneid(bv_entity) == 1

## Test getters for blinker noise and goals
@test blinker(bv_entity)
@test blinker(bv_entity.state)
@test !blinker(ped_entity)
@test !blinker(ped_entity.state)

@test goals(bv_entity) == [3,4]
@test goals(bv_entity.state) == [3,4]
@test goals(ped_entity) == [2]
@test goals(ped_entity.state) == [2]

@test noise(bv_entity) == neg_noise
@test noise(bv_entity.state) ==  neg_noise
@test noise(ped_entity) == pos_noise
@test noise(ped_entity.state) == pos_noise

## Test update_veh_state
new_state = VehicleState(VecSE2(50., 50., 0.), ped_roadway, 30.)

new_ped = update_veh_state(ped, new_state)
@test new_ped isa NoisyPedState
@test new_ped.veh_state == new_state
@test noise(new_ped) == pos_noise

new_veh = update_veh_state(veh, new_state)
@test new_veh isa BlinkerState
@test new_veh.veh_state == new_state
@test noise(new_veh) == neg_noise

## Test the noise

# Check the noisy_entity function
noisy_ped = noisy_entity(ped_entity, ped_roadway)
@test laneid(noisy_ped) == 2
@test posf(noisy_ped).s == posf(ped_entity).s + noise(ped_entity).pos[1]
@test posf(noisy_ped).t == posf(ped_entity).t + noise(ped_entity).pos[2]
@test posf(noisy_ped).ϕ == posf(ped_entity).ϕ
@test vel(noisy_ped) == vel(ped_entity) + noise(ped_entity).vel
@test posg(noisy_ped).x == 24.9
@test posg(noisy_ped).y ≈ -9.7

noisy_veh = noisy_entity(bv_entity, ped_roadway)
@test laneid(noisy_veh) == 1
@test posf(noisy_veh).s == posf(bv_entity).s + noise(bv_entity).pos[1]
@test posf(noisy_veh).t == posf(bv_entity).t + noise(bv_entity).pos[2]
@test posf(noisy_veh).ϕ == posf(bv_entity).ϕ
@test vel(noisy_veh) == vel(bv_entity) + noise(bv_entity).vel
@test posg(noisy_veh).x ≈ 49.7
@test posg(noisy_veh).y ≈ -0.1

# Check the noisy_scene function
scene = Scene([ped_entity, bv_entity])
n_scene = noisy_scene(scene, ped_roadway)
@test n_scene[1] == noisy_ped
@test n_scene[2] == noisy_veh

## Check constructors

# Test the BlinkerVehicle Constructor
veh = BlinkerVehicle(roadway = Tint_roadway, lane = 5, s = 30., v = 20., id = 1, goals = Tint_goals[5], blinker = true)
state = VehicleState(VecSE2(polar(20,-π/2) + VecE2(1.5,0), π/2), 20.)
@test all(isapprox.(posg(veh), posg(state)))
@test vel(veh) == vel(state)
@test posf(veh) == posf(veh.state.veh_state)
@test blinker(veh)
@test goals(veh) == [5,6]
@test laneid(veh) == 5
@test can_have_goal(veh, 6, Tint_roadway)
@test can_have_goal(veh, 5, Tint_roadway)
@test !can_have_goal(veh, 1, Tint_roadway)
@test !can_have_goal(veh, 2, Tint_roadway)

# Test the pedestrian Constructor
ped = NoisyPedestrian(roadway = ped_roadway, lane = 2, s = 2., v = 2., id = 1)
@test posf(ped).s == 2.
@test posg(ped).y ≈ -8.
@test posg(ped).x == 25.
@test vel(ped) == 2.
@test laneid(ped) == 2

## Check the behavior of the actions for the Pedestrian Control
simple_acc = PedestrianControl(a = (0.1, 0.1))
acc_state = propagate(ped_entity, simple_acc, ped_roadway, 0.1)
@test vel(acc_state) > vel(ped_entity)
@test posf(acc_state).s > posf(ped_entity).s
@test posf(acc_state).t > posf(ped_entity).t
@test posf(acc_state).ϕ > posf(ped_entity).ϕ

test_action = PedestrianControl(a = (0.2, 0.2), da = (-0.1,-0.1), noise = neg_noise)
new_state = propagate(ped_entity, test_action, ped_roadway, 0.1)
new_ped = Entity(new_state, PEDESTRIAN_DEF, 1)
@test laneid(new_ped) == 2
@test vel(acc_state) == vel(new_state)
@test all(posg(acc_state) .== posg(new_state))
@test noise(new_state) == test_action.noise

## Check the behavior of actions for BlinkerVehicle
simple_acc = BlinkerVehicleControl(a = 0.1)
acc_state = propagate(veh, simple_acc, Tint_roadway, 0.1)
@test vel(acc_state) > vel(veh)
@test posf(acc_state).s > posf(veh).s
@test posf(acc_state).t == posf(veh).t
@test posf(acc_state).ϕ == posf(veh).ϕ

test_action = BlinkerVehicleControl(a = 0.2, da = -0.1, toggle_goal = true, toggle_blinker = true, noise = neg_noise)
new_state = propagate(veh, test_action, Tint_roadway, 0.1)
new_veh = Entity(new_state, VehicleDef(), 1)
@test laneid(new_veh) == 6
@test blinker(new_state) == false
@test vel(acc_state) == vel(new_state)
@test all(posg(acc_state) .== posg(new_state))
@test noise(new_state) == test_action.noise

## Test the models

# Test the contruction of an IDM Model
Tmodel = TIDM(  yields_way = Tint_yields_way,
                intersection_enter_loc = Tint_intersection_enter_loc,
                intersection_exit_loc = Tint_intersection_exit_loc,
                goals = Tint_goals,
                should_blink = Tint_should_blink
             )

@test Tmodel.idm isa IntelligentDriverModel
@test Tmodel.noisy_observations == false
test_action = BlinkerVehicleControl(a = 0.2, da = -0.1, toggle_goal = false, toggle_blinker = true, noise = neg_noise)
Tmodel.next_action = test_action
track_longitudinal!(Tmodel.idm, 15., 0., 10.)
actual_action = rand(Random.GLOBAL_RNG, Tmodel)
@test actual_action.a == Tmodel.idm.a
@test actual_action.da == test_action.da
@test actual_action.toggle_goal == test_action.toggle_goal
@test actual_action.toggle_blinker == test_action.toggle_blinker
@test actual_action.noise == test_action.noise

# Test the construction of an AdversarialPedestrian model
ped_model = AdversarialPedestrian()
@test ped_model.idm isa IntelligentDriverModel
ped_model.next_action = PedestrianControl(a = (0.2, 0.2), da = (-0.1,-0.1), noise = neg_noise)
track_longitudinal!(ped_model.idm, 15., 0., 10.)
actual_action = rand(Random.GLOBAL_RNG, ped_model)
@test actual_action.a[1] == ped_model.idm.a
@test actual_action.a[2] == 0.
@test all(isapprox.(actual_action.da, (-0.1,-0.1)))
@test actual_action.noise == neg_noise

## Test lane prediction
veh1 = BlinkerVehicle(roadway = Tint_roadway, lane = 2, s = 40., v = 20., id = 1, goals = Tint_goals[2], blinker = false)
veh2 = BlinkerVehicle(roadway = Tint_roadway, lane = 2, s = 40., v = 20., id = 1, goals = Tint_goals[2], blinker = true)
veh3 = BlinkerVehicle(roadway = Tint_roadway, lane = 2, s = 50., v = 20., id = 1, goals = Tint_goals[2], blinker = false)
veh4 = BlinkerVehicle(roadway = Tint_roadway, lane = 2, s = 50., v = 20., id = 1, goals = Tint_goals[2], blinker = true)

@test lane_belief(veh1, Tmodel, Tint_roadway) == 2
@test lane_belief(veh2, Tmodel, Tint_roadway) == 1
@test lane_belief(veh3, Tmodel, Tint_roadway) == 2
@test lane_belief(veh4, Tmodel, Tint_roadway) == 2

## Check collisions
egovehicle = BlinkerVehicle(roadway = Tint_roadway, lane = 5, s = 50., v = 9., id = 5, goals = Tint_goals[5], blinker = true)
egovehicle2 = BlinkerVehicle(roadway = Tint_roadway, lane = 5, s = 35., v = 9., id = 5, goals = Tint_goals[5], blinker = true)
vehicles = [BlinkerVehicle(roadway = Tint_roadway, lane = 4, s = 50., v = 20., id = 1, goals = Tint_goals[5], blinker = false),
            BlinkerVehicle(roadway = Tint_roadway, lane = 2, s = 50., v = 19., id = 2, goals = Tint_goals[5], blinker = false),
            BlinkerVehicle(roadway = Tint_roadway, lane = 2, s = 10., v = 19., id = 3, goals = Tint_goals[5], blinker = false),
            BlinkerVehicle(roadway = Tint_roadway, lane = 4, s = 40., v = 14., id = 4, goals = Tint_goals[5], blinker = true)]

scene_ego = Scene([egovehicle, vehicles[2:end]...])
scene_any = Scene([egovehicle2, vehicles...])

@test any_collides(scene_any)
@test any_collides(scene_ego)
@test ego_collides(5, scene_ego)
@test !ego_collides(3, scene_any)

egovehicle = BlinkerVehicle(roadway = ped_roadway, lane = 1, s = 25., v=20., id = 2, blinker = false, goals = ped_goals[1])
pedestrian = NoisyPedestrian(roadway = ped_roadway, lane = 2, s = 9.5, v = 2., id = 1)
scene = Scene([pedestrian, egovehicle])

@test any_collides(scene)
@test ego_collides(2, scene)

## Construct a whole scene and let it run (T-Intersection)
egovehicle = BlinkerVehicle(roadway = Tint_roadway, lane = 5, s = 35., v = 9., id = 5, goals = Tint_goals[5], blinker = true)
vehicles = [BlinkerVehicle(roadway = Tint_roadway, lane = 3, s = 30., v = 20., id = 1, goals = Tint_goals[5], blinker = false),
            BlinkerVehicle(roadway = Tint_roadway, lane = 2, s = 15., v = 19., id = 2, goals = Tint_goals[5], blinker = false),
            BlinkerVehicle(roadway = Tint_roadway, lane = 2, s = 10., v = 19., id = 3, goals = Tint_goals[5], blinker = false),
            BlinkerVehicle(roadway = Tint_roadway, lane = 4, s = 40., v = 14., id = 4, goals = Tint_goals[5], blinker = true)]

scene = Scene([egovehicle, vehicles...])

@test !any_collides(scene)
@test !ego_collides(5, scene)

models = Dict(i => TIDM(Tint_TIDM_template) for i=1:5)

nticks, timestep = 100, 0.1
scenes = AutomotiveSimulator.simulate(scene, Tint_roadway, models, nticks, timestep)
@test length(scenes) == nticks + 1

# test end_of_road
@test end_of_road(scenes[end][1], Tint_roadway, Inf)
@test end_of_road(scenes[end][2], Tint_roadway, Inf)
@test end_of_road(scenes[end][3], Tint_roadway, Inf)
@test end_of_road(scenes[end][4], Tint_roadway, Inf)
@test end_of_road(scenes[end][5], Tint_roadway, Inf)

egovehicle = BlinkerVehicle(roadway = Tint_roadway, lane = 5, s = 35., v = 9., id = 5, goals = Tint_goals[5], blinker = true)
@test end_of_road(egovehicle, Tint_roadway, 35)

# Note: use the following code to make a gif
# The cars are not removed at the end of the roadway, the mpd takes care of that
# using Reel
# animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
#     i = Int(floor(t/dt)) + 1
#     render([Tint_roadway, scenes[i]], canvas_width=1200, canvas_height=800)
# end
# write("Tint_roadway_animated.gif", animation)

## Construct a scene and let it run with the pedestrian and the car
pedestrian = NoisyPedestrian(roadway = ped_roadway, lane = 2, s = 5., v = 2., id = 1)
egovehicle = BlinkerVehicle(roadway = ped_roadway, lane = 1, s = 0., v=10., id = 2, blinker = false, goals = ped_goals[1])
scene = Scene([pedestrian, egovehicle])

models = Dict(1=>AdversarialPedestrian(idm = IntelligentDriverModel(v_des = 2)), 2 => ped_TIDM_template)

nticks = 50
timestep = 0.1
scenes = AutomotiveSimulator.simulate(scene, ped_roadway, models, nticks, timestep)
@test length(scenes) == nticks + 1

# animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
#     i = Int(floor(t/dt)) + 1
#     render([ped_roadway, crosswalk, scenes[i]], canvas_width=1200, canvas_height=800)
# end
# write("ped_roadway_animated.gif", animation)

## Run where noise causes the system to fail
pedestrian = NoisyPedestrian(roadway = ped_roadway, lane = 2, s = 5., v = 2., id = 1, noise = Noise(pos=(-100, 0), vel=-3.))
egovehicle = BlinkerVehicle(roadway = ped_roadway, lane = 1, s = 0., v=10., id = 2, blinker = false, goals = ped_goals[1])
scene = Scene([pedestrian, egovehicle])

vehicle_model = TIDM(ped_TIDM_template, noisy_observations = true)
ped_model = AdversarialPedestrian(
            idm = IntelligentDriverModel(v_des = 2),
            next_action = PedestrianControl(noise =  Noise(pos=(-100, 0)))
            )
models = Dict(1 => ped_model,  2 => vehicle_model)

nticks = 50
timestep = 0.1
scenes = AutomotiveSimulator.simulate(scene, ped_roadway, models, nticks, timestep)
@test length(scenes) == nticks + 1

# animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
#     i = Int(floor(t/dt)) + 1
#     render([ped_roadway, crosswalk, scenes[i]], canvas_width=1200, canvas_height=800)
# end
# write("ped_roadway_animated.gif", animation)

## Test that car tracks the object if it is outside the intersection
pedestrian = NoisyPedestrian(roadway = ped_roadway, lane = 2, s = 5., t=10., v = 2., id = 1, )
egovehicle = BlinkerVehicle(roadway = ped_roadway, lane = 1, s = 0., v=10., id = 2, blinker = false, goals = ped_goals[1])

scene = Scene([pedestrian, egovehicle])

vehicle_model = TIDM(ped_TIDM_template)
ped_model = AdversarialPedestrian(idm = IntelligentDriverModel(v_des = 2))
models = Dict(1 => ped_model,  2 => vehicle_model)

nticks = 50
timestep = 0.1
scenes = AutomotiveSimulator.simulate(scene, ped_roadway, models, nticks, timestep)
@test length(scenes) == nticks + 1

# animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
#     i = Int(floor(t/dt)) + 1
#     render([ped_roadway, crosswalk, scenes[i]], canvas_width=1200, canvas_height=800)
# end
# write("ped_roadway_animated.gif", animation)

## Test of the blindspot
@test !in_blindspot(VecSE2(5, 0, 0.), Blindspot(π/12., π/6), VecSE2(25, -3, π/2))
@test in_blindspot(VecSE2(5, 0, -π/6), Blindspot(π/12., π/6), VecSE2(25, -3, π/2))
@test !in_blindspot(VecSE2(5, 0, π/6), Blindspot(π/12., π/6), VecSE2(25, -3, π/2))
@test in_blindspot(VecSE2(5, 0, 0.), Blindspot(π/12., π/3), VecSE2(25, -3, π/2))

## Construct a disturbance model for the adversaries
disturbances = Sampleable[Normal(0,1), Bernoulli(0.5), Bernoulli(0.5), Bernoulli(0.), Bernoulli(0.)]
sut_agent = BlinkerVehicleAgent(rand_up_left(id=1, s_dist=Uniform(15.,25.), v_dist=Uniform(10., 29)), TIDM(Tint_TIDM_template))
right_adv = BlinkerVehicleAgent(rand_right(id=3, s_dist=Uniform(15.,35.), v_dist=Uniform(15., 29.)), TIDM(Tint_TIDM_template), disturbance_model = disturbances)
mdp = DrivingMDP(sut_agent, [right_adv], Tint_roadway, 0.2, γ = 0.95)

@test length(agents(mdp)) == 2
@test mdp.per_timestep_penalty == 0
@test mdp.v_des == 25


text:
# y: data
# yhat_generator: function that samples data given the parameters
# algo parameters: Namedtuple of algorithm parameters
# max_time: max time allowed for the algo
# N_samples: number of samples desired
function ABC(y,yhat_generator,algo_parameters,N_samples)
    CPUtic()

    # TODO: implement max_time
    eta_y = algo_parameters[:eta](y)
    q = length(algo_parameters[:prior])
    result = zeros(N_samples,q)
    size_y = size(y)
    i = 1

    while i <= N_samples
        # sample parameters
        p = rand.(algo_parameters[:prior])

        # generate sample
        yhat = yhat_generator(p)

        # in case there is no soln given the parameters
        if size_y != size(yhat)
            continue
        end

        # check acceptance
        if algo_parameters[:d](eta_y,algo_parameters[:eta](yhat)) < algo_parameters[:epsilon]
            result[i,:] .= p
            i += 1
        end

    end
    t_end = CPUtoc_modified(false)
    return result, t_end
end