text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ 4dd67540-6696-11eb-313d-49aaec40c97b
using ScikitLearn

# ╔═╡ 07836df0-669b-11eb-36e1-c7a130d62f2a
using Statistics

# ╔═╡ ba206ea0-6695-11eb-0046-77b8e1fcd083
begin
    using PlutoUI
    using CSV # This is a pacakge we use for loading CSV Files.
    using DataFrames
    @sk_import model_selection:train_test_split
    @sk_import linear_model:LogisticRegression
    @sk_import metrics:confusion_matrix
    @sk_import metrics:classification_report
    @sk_import model_selection:cross_val_score
    @sk_import metrics:plot_confusion_matrix
    pathtodata = joinpath("covid_cleaned.csv")
    data = CSV.read(pathtodata, DataFrame)
end

# ╔═╡ bce85980-669b-11eb-0ce0-051cccead17d
using PyPlot

# ╔═╡ 84be3e80-6696-11eb-1c2d-695ff4573f91
md"""
The dataset looks like this 👇. By clicking on the `more` button in the table you can see more rows and columns of the data"""

# ╔═╡ 751bd5c0-6699-11eb-1ea5-51bcfdab90a1
md""" ### Our Dataset"""

# ╔═╡ 74483970-6696-11eb-0da7-2fa268ef4d25
begin
    X = convert(Array, data[!, Not(:covid_res)])
    y = convert(Array, data[!, :covid_res])
end;

# ╔═╡ 1da3a1d0-6697-11eb-32e7-d7b544ffbb28
md"""Our **target variable** is a column of zeros and ones:"""

# ╔═╡ fcdf6b00-6696-11eb-3b70-a94d2316c8d7
DataFrame(; covid_res=y)

# ╔═╡ 41745550-6697-11eb-062a-35fad903d41c
md""" #### Our Input variables
"""

# ╔═╡ 3b2027b0-6697-11eb-3a5e-7fa81a886201
DataFrame(X, names(data[!, Not(:covid_res)]))

# ╔═╡ 8d7322fe-6697-11eb-0280-cb7d8ce27fa2
md"""The below line of code is splitting the data into a training set and test set"""

# ╔═╡ 7d973b10-6697-11eb-36af-b94eb1278ddd
X_train, X_test, y_train, y_test = train_test_split(X, y; test_size=0.5, random_state=42);

# ╔═╡ 10d8d500-6698-11eb-285f-dd02bc6eb853
md"""##### Taking a look at the trainning X"""

# ╔═╡ 0336b070-6698-11eb-3cbe-a92e5ae467be
DataFrame(X_train, names(data[!, Not(:covid_res)]))

# ╔═╡ 4f1f4190-6699-11eb-285f-5d9a5dd0f7bc
md"""#### Looking at Test X"""

# ╔═╡ 071a0520-6698-11eb-353f-3b96d46773e1
DataFrame(X_test, names(data[!, Not(:covid_res)]))

# ╔═╡ 4b93d0e0-6699-11eb-21b3-9df9ce7bbe44
md"""#### Y Train"""

# ╔═╡ 63ea6ff0-6699-11eb-0384-49a7e5cce81d
DataFrame(; covid_res=y_train)

# ╔═╡ a8832ad0-6699-11eb-348f-8ddd45591076
md""" ### Y Test"""

# ╔═╡ b3440160-6699-11eb-30b1-c18f20f4ae6f
DataFrame(; covid_res=y_test)

# ╔═╡ f27fd1ce-6697-11eb-0bf7-293f5fe447e2
md""" #### Fitting the Model to Data"""

# ╔═╡ e259c180-6697-11eb-3ce6-8506c5dc4e3f
begin
    simplelogistic = LogisticRegression(; max_iter=300)
    fit!(simplelogistic, X_train, y_train)
end;

# ╔═╡ da6bf630-6699-11eb-0668-33d393e9b2ac
begin
    y_pred = predict(simplelogistic, X_train)

    with_terminal() do
        print(classification_report(y_train, y_pred))
    end
end

# ╔═╡ a79700a0-669a-11eb-0aa8-4f8ed5b18798
c = cross_val_score(LogisticRegression(), X_train, y_train; cv=5)

# ╔═╡ de055420-669a-11eb-2dfd-a5fdbe772ffb
md"""**Average acurracy of the model with cross validation on training dataset is $(mean(c)) and SD is $(std(c))**"""

# ╔═╡ 97e42f10-669b-11eb-28fc-4bf98f97d749
ctest = cross_val_score(LogisticRegression(), X_test, y_test; cv=5)

# ╔═╡ a74c6da0-669b-11eb-1426-618da4ecb341
md"""**Average acurracy of the model with cross validation on test dataset is $(mean(ctest)) and SD is $(std(ctest))**"""

# ╔═╡ d53e0e30-669b-11eb-0928-637f482b2a51
md"""### Confusion Matrix"""

# ╔═╡ c61d5ff0-669b-11eb-2584-1fb0fd92f02a
begin
    plot_confusion_matrix(simplelogistic, X_train, y_train)
    gcf()
end

# ╔═╡ Cell order:
# ╠═4dd67540-6696-11eb-313d-49aaec40c97b
# ╠═07836df0-669b-11eb-36e1-c7a130d62f2a
# ╟─84be3e80-6696-11eb-1c2d-695ff4573f91
# ╟─751bd5c0-6699-11eb-1ea5-51bcfdab90a1
# ╟─ba206ea0-6695-11eb-0046-77b8e1fcd083
# ╠═74483970-6696-11eb-0da7-2fa268ef4d25
# ╟─1da3a1d0-6697-11eb-32e7-d7b544ffbb28
# ╟─fcdf6b00-6696-11eb-3b70-a94d2316c8d7
# ╟─41745550-6697-11eb-062a-35fad903d41c
# ╟─3b2027b0-6697-11eb-3a5e-7fa81a886201
# ╟─8d7322fe-6697-11eb-0280-cb7d8ce27fa2
# ╠═7d973b10-6697-11eb-36af-b94eb1278ddd
# ╟─10d8d500-6698-11eb-285f-dd02bc6eb853
# ╟─0336b070-6698-11eb-3cbe-a92e5ae467be
# ╟─4f1f4190-6699-11eb-285f-5d9a5dd0f7bc
# ╟─071a0520-6698-11eb-353f-3b96d46773e1
# ╟─4b93d0e0-6699-11eb-21b3-9df9ce7bbe44
# ╟─63ea6ff0-6699-11eb-0384-49a7e5cce81d
# ╟─a8832ad0-6699-11eb-348f-8ddd45591076
# ╟─b3440160-6699-11eb-30b1-c18f20f4ae6f
# ╟─f27fd1ce-6697-11eb-0bf7-293f5fe447e2
# ╠═e259c180-6697-11eb-3ce6-8506c5dc4e3f
# ╟─da6bf630-6699-11eb-0668-33d393e9b2ac
# ╠═a79700a0-669a-11eb-0aa8-4f8ed5b18798
# ╟─de055420-669a-11eb-2dfd-a5fdbe772ffb
# ╠═97e42f10-669b-11eb-28fc-4bf98f97d749
# ╟─a74c6da0-669b-11eb-1426-618da4ecb341
# ╟─bce85980-669b-11eb-0ce0-051cccead17d
# ╟─d53e0e30-669b-11eb-0928-637f482b2a51
# ╟─c61d5ff0-669b-11eb-2584-1fb0fd92f02a

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
"""
    grdview(cmd0::String="", arg1=nothing, arg2=nothing, arg3=nothing; kwargs...)

Reads a 2-D grid and produces a 3-D perspective plot by drawing a mesh, painting a
colored/grayshaded surface made up of polygons, or by scanline conversion of these polygons
to a raster image.

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

- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_R)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **G** | **drape** | **drapefile** :: [Type => Str | GMTgrid | a Tuple with 3 GMTgrid types]

    Drape the image in drapefile on top of the relief provided by relief_file.
    ($(GMTdoc)grdview.html#g)
- **I** | **shade** | **shading** | **intensity** :: [Type => Str | GMTgrid]		``Arg = GMTgrid | filename``

    Gives the name of a grid file or GMTgrid with intensities in the (-1,+1) range,
    or a grdgradient shading flags.
    ($(GMTdoc)grdview.html#i)
- **N** | **plane** :: [Type => Str | Int]		``Arg = (level [,fill])``

    Draws a plane at this z-level.
    ($(GMTdoc)grdview.html#n)
- $(GMT.opt_P)
- **Q** | **surftype** | **surf** :: [Type => Str | Int] ``Arg = mesh=Bool, surface=Bool, image=Bool, wterfall=(:rows|cols,[fill])``

    Specify **m** for mesh plot, **s** for surface, **i** for image.
    ($(GMTdoc)grdview.html#q)
- **S** | **smoothfactor** :: [Type => Number]

    Used to resample the contour lines at roughly every (gridbox_size/smoothfactor) interval..
    ($(GMTdoc)grdview.html#s)
- **T** | **tiles** | **no_interp** :: [Type => Str | NT]	``Arg = (skip|skip_nan=Bool, outlines=Bool|pen)``

    Plot image without any interpolation.
    ($(GMTdoc)grdview.html#t)
- **W** | **pens** | **pen** :: [Type => Str]	``Arg = (contour=Bool|pen, mesh=Bool|pen, facade=Bool|pen)``

    Draw contour, mesh or facade. Append pen attributes.
    ($(GMTdoc)grdview.html#w)

- **isgeog** :: [Type => Any]

    When drapping an image that has projection info over a grid that is in geographics but does not carry any
    information about this fact we may need to use this option to help the program finding the common BoundingBox.
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_f)
- $(GMT.opt_n)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function grdview(cmd0::String="", arg1=nothing; first=true, kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdview", cmd0, arg1)
	arg2 = nothing;	arg3 = nothing;	arg4 = nothing;	arg5 = nothing;

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode
	common_insert_R!(d, O, cmd0, arg1)			# Set -R in 'd' out of grid/images (with coords) if limits was not used

	has_opt_B = (find_in_dict(d, [:B :frame :axis :axes], false)[1] !== nothing)
	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "grdview", O, " -JX" * split(def_fig_size, '/')[1] * "/0")
	(startswith(opt_J, " -JX") && !contains(opt_J, "/")) && (cmd = replace(cmd, opt_J => opt_J * "/0")) # When sub-regions
	(!has_opt_B && isa(arg1, GMTimage) && (isimgsize(arg1) || CTRL.limits[1:4] == zeros(4)) && opt_B == def_fig_axes) &&
		(cmd = replace(cmd, opt_B => ""))	# Dont plot axes for plain images if that was not required

	cmd, = parse_common_opts(d, cmd, [:UVXY :c :f :n :p :t :params], first)
	cmd  = add_opt(d, cmd, 'S', [:S :smooth])
	if ((val = find_in_dict(d, [:N :plane])[1]) !== nothing)
		cmd *= " -N" * parse_arg_and_pen(val, "+g", false)
	end
	cmd = add_opt(d, cmd, 'Q', [:Q :surf :surftype],
				  (mesh=("m", add_opt_fill), surface="_s", surf="_s", img=("i",arg2str), image="i", nan_alpha="_c", monochrome="_+m", waterfall=(rows="my", cols="mx", fill=add_opt_fill)))
	cmd = add_opt(d, cmd, 'W', [:W :pens :pen], (contour=("c", add_opt_pen),
	              mesh=("m", add_opt_pen), facade=("f", add_opt_pen)) )
	cmd = add_opt(d, cmd, 'T', [:T :no_interp :tiles], (skip="_+s", skip_nan="_+s", outlines=("+o", add_opt_pen)) )
	(!occursin(" -T", cmd)) ? cmd = parse_JZ(d, cmd)[1] : del_from_dict(d, [:JZ])	# Means, even if we had one, ignore silently
	cmd = add_opt(d, cmd, "%", [:layout :mem_layout], nothing)

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)		# Find how data was transmitted

	(isa(arg1, Array{<:Real})) && (arg1 = mat2grid(arg1))

	cmd, N_used, arg1, arg2, arg3 = common_get_R_cpt(d, cmd0, cmd, opt_R, got_fname, arg1, arg2, arg3, "grdview")
	cmd, arg1, arg2, arg3, arg4 = common_shade(d, cmd, arg1, arg2, arg3, arg4, "grdview")
	cmd, arg1, arg2, arg3, arg4, arg5 = parse_G_grdview(d, [:G :drape :drapefile], cmd0, cmd, arg1, arg2, arg3, arg4, arg5)

	_cmd, K = finish_PS_nested(d, ["grdview " * cmd], K)
	finish_PS_module(d, _cmd, "", K, O, true, arg1, arg2, arg3, arg4, arg5)
end

# ---------------------------------------------------------------------------------------------------
function parse_G_grdview(d::Dict, symbs::Array{<:Symbol}, cmd0::String, cmd::String, arg1, arg2, arg3, arg4, arg5)
	(show_kwargs[1]) && return print_kwarg_opts(symbs, "GMTgrid | Tuple | String"), arg1, arg2, arg3, arg4, arg5
	if ((val = find_in_dict(d, symbs)[1]) !== nothing)
		function range_it(val)
			cmd, N_used = put_in_slot(cmd, val, 'G', [arg1, arg2, arg3, arg4])
			if     (N_used == 1)  arg1 = val
			elseif (N_used == 2)  arg2 = val
			elseif (N_used == 3)  arg3 = val
			elseif (N_used == 4)  arg4 = val
			end
			return cmd, arg1, arg2, arg3, arg4
		end

		if (isa(val, String) || isa(val, GMTimage))
			if (isa(val, String) && guess_T_from_ext(val) != " -Ti")
				cmd *= " -G" * val
			else
				if (cmd0 != "")
					prj = (startswith(cmd0, "@earth_r")) ? prj4WGS84 : getproj(cmd0, proj4=true)
				else
					prj = (isa(arg1, GMTgrid)) ? getproj(arg1) : ""
				end
				t = split(scan_opt(cmd, "-R"), '/')
				Iname = drape_prepare(d, val, ["-projwin", t[1], t[4], t[2], t[3]], prj)
				cmd *= " -G" * Iname
				(!contains(cmd, " -Qi")) && (cmd *= " -Qi300")	# Otherwise GMT crashes because grdview goes through the "MESH" branch
			end
		elseif (isa(val, GMTgrid))			# A single drape grid (arg1-3 may be used already)
			cmd, arg1, arg2, arg3, arg4 = range_it(val)
		elseif (isa(val, Tuple) && length(val) == 3)
			cmd, N_used = put_in_slot(cmd, val[1], 'G', [arg1, arg2, arg3, arg4, arg5])
			cmd *= " -G -G"					# Because the above only set one -G and we need 3
			if     (N_used == 1)  arg1 = val[1];	arg2 = val[2];		arg3 = val[3]
			elseif (N_used == 2)  arg2 = val[1];	arg3 = val[2];		arg4 = val[3]
			elseif (N_used == 3)  arg3 = val[1];	arg4 = val[2];		arg5 = val[3]
			end
		else
			error("Wrong way of setting the drape (G) option.")
		end
	end
	return cmd, arg1, arg2, arg3, arg4, arg5
end

# ---------------------------------------------------------------------------------------------------
function drape_prepare(d::Dict, fname, opts::Vector{AbstractString}, prj::String)
	# Deal with the option of drapping an image, which can be smaller, larger or with fifferent projection.
	prj_img = getproj(fname, proj4=true)
	(prj_img == "") && return fname			# If drape image has no RefSys just return its name and let it all be used

	(prj == "" && find_in_dict(d, [:isgeog])[1] !== nothing) && (prj = prj4WGS84)

	# Layout is "TRB" so all matrices are contrary to Julia order. opts=[-projwin xmin ymax xmax ymin]
	G_I = (prj == "" || prj == prj_img) ? gdaltranslate(fname, opts) :
	                                      gdalwarp(fname, ["-t_srs", prj, "-te", opts[2], opts[5], opts[4], opts[3]])

	W = parse(Float64, opts[2]);	E = parse(Float64, opts[4])
	S = parse(Float64, opts[5]);	N = parse(Float64, opts[3])
	Wbk, Ebk, Sbk, Nbk = W, E, S, N
	dx_W = G_I.range[1] - W;		dx_E = G_I.range[2] - E
	dy_S = G_I.range[3] - S;		dy_N = G_I.range[4] - N
	pad_W = ceil(Int, abs(dx_W) / G_I.inc[1]);		pad_E = ceil(Int, abs(dx_E) / G_I.inc[1])
	pad_S = ceil(Int, abs(dy_S) / G_I.inc[2]);		pad_N = ceil(Int, abs(dy_N) / G_I.inc[2])
	# Recompute the WESN such that the increments don't change (original -R was in GRID increment multiples)
	W = G_I.range[1] - pad_W * G_I.inc[1];			E = G_I.range[2] + pad_E * G_I.inc[1]
	S = G_I.range[3] - pad_S * G_I.inc[2];			N = G_I.range[4] + pad_N * G_I.inc[2]
	if (abs(W - Wbk) > 2*G_I.inc[1] || abs(E - Ebk) > 2*G_I.inc[1] || abs(S - Sbk) > 2*G_I.inc[2] || abs(N - Nbk) > 2*G_I.inc[2])
		img_is_geo = contains(prj_img, "longlat") || contains(prj_img, "latlong")
		msg = "Grid and draping image do not carry enough information about their Referencing Systems"
		(!img_is_geo) && (msg *= "\n\tDraping image is not in geogs but grid probably is. Please use option `isgeog=true`")
		error(msg)
	end
	if (pad_W > 0 || pad_E > 0 || pad_S > 0 || pad_N > 0)
		img_new = (size(G_I, 3) == 1) ? fill(UInt8(255), size(G_I,1)+pad_W+pad_E, size(G_I,2)+pad_S+pad_N) :
										fill(UInt8(255), size(G_I,1)+pad_W+pad_E, size(G_I,2)+pad_S+pad_N, size(G_I,3))
		n = 0
		for l = 1:size(img_new,3)
			@simd for row = pad_N+1:(size(G_I,2)+pad_N)
				@simd for col = pad_W+1:(size(G_I,1)+pad_W)
					@inbounds img_new[col,row,l] = G_I.image[n += 1]
				end
			end
		end
		G_I = mat2img(img_new, G_I)
		G_I.x = linspace(W, E, size(img_new,1)+G_I.registration)
		G_I.y = linspace(S, N, size(img_new,2)+G_I.registration)
		G_I.inc = [G_I.x[2]-G_I.x[1], G_I.y[2]-G_I.y[1]]
		G_I.range[1:4] = [W, E, S, N]
	end

	#def_name = joinpath(tempdir(), "GMTjl_2grdview.tiff")
	def_name = "/vsimem/GMTjl_2grdview.tiff"	# I'm amazed that this works
	gdalwrite(def_name, G_I)
	ressurectGDAL()
	return def_name
end

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

ext: jl


text:
# The number 3797 has an interesting property. Being prime itself, it is
# possible to continuously remove digits from left to right, and remain prime
# at each stage: 3797, 797, 97, and 7. Similarly we can work from right to
# left: 3797, 379, 37, and 3.
#
# Find the sum of the only eleven primes that are both truncatable from left to
# right and right to left.
#
# NOTE: 2, 3, 5, and 7 are not considered to be truncatable primes.

using ProjectEulerSolutions

# Helper function to generate an array of truncated integers given an input
# integer.
function truncatednumbers(n::Integer)::Array{Integer}
    numdigits = ndigits(n) - 1

    n_left = n
    n_right = n
    trunc_nums = Integer[]

    while numdigits > 0
        n_left = fld(n_left, 10)
        n_right = n_right % (10^numdigits)
        push!(trunc_nums, n_left)
        push!(trunc_nums, n_right)
        numdigits -= 1
    end

    return trunc_nums
end

# Use Sieve of Eratosthenes to get primes, store in set, then cycle through
# them creating truncated numbers and checking whether all truncated numbers
# are prime.  If so, move the prime into the verified set.
function p037solution(n::Integer=10)::Integer

    primes = Set{Integer}(sieve_eratosthenes(1_000_000))
    verified = Set{Integer}()

    num_truncs = 0
    for prime in primes
        if num_truncs == n
            break
        end
        if ndigits(prime) > 1
            digs = digits(prime)
            if digs[end] == 9 || mapreduce(x-> x % 2, *, digs[1:end-1]) == 0
                # If any digit (other than the first) is even or first digit is
                # a 9, ignore number
                continue
            end

            truncnums = truncatednumbers(prime)
            if mapreduce(prime -> in(prime, primes), &, truncnums)
                push!(verified, prime)
                num_truncs += 1
            end
        end
    end

    return sum(verified)
end

p037 = Problems.Problem(p037solution)

Problems.benchmark(p037, 11)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# test atomic units equiality based in valence
let equality_tests = [
        :((C().valence == H().valence) == false)
        :((H().valence == H().valence) == true)
        :((C() == H()) == false)
        :((C() == C()) == true)
        :((C(10) == H(10)) == false)
        :((C(10) == C(15)) == true)
    ]
    for test in equality_tests
        @eval @test $test
    end
end

# test atomic units constructors against every native numeric type
let n = 1
    nums_data = [
        (int8(n), Int8),
        (uint8(n), Uint8),
        (int16(n), Int16),
        (uint16(n), Uint16),
        (int32(n), Int32),
        (uint32(n), Uint32),
        (int64(n), Int64),
        (uint64(n), Uint64),
        (int128(n), Int128),
        (uint128(n), Uint128),
        (true, Bool),
        (complex(true, true), Bool),
        (float16(n), Float16),
        (float32(n), Float32),
        (float64(n), Float64),
        (BigInt(n), BigInt),
        (BigFloat(n), BigFloat),
        (n//n, Rational{Int}),
        (n//n + n//big(n)*im, Rational{BigInt}),
        (n+n*im, Int),
        (Inf*im, Float64)
    ]
    for (value, T) in nums_data
        @test Carbon{T} == typeof(C(value))
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CUTEst"
version = v"2.0.3"

# Collection of sources required to build ThinASLBuilder
sources = [
    ArchiveSource("https://github.com/ralna/ARCHDefs/archive/v2.0.3x.tar.gz",
                  "6583e27f84338447767bbdf4335514c8836ae4ad54f5e66280307e8b57189cff"),
    ArchiveSource("https://github.com/ralna/SIFDecode/archive/v2.0.3.tar.gz",
                  "3a4aa817e1bf4e3595d0e4378da6172b65f02861f3a7c39f9da632a5cc31b1b2"),
    ArchiveSource("https://github.com/ralna/CUTEst/archive/v2.0.3.tar.gz",
                  "d21a65c975302296f9856c09034cf46edc5da34b6efd96eed6cc94af6d2c8a55"),
]

# Bash recipe for building across all platforms
script = raw"""
echo "building for ${target}"

# setup
mkdir -p ${prefix}/libexec
mkdir -p ${bindir}
mkdir -p ${libdir}
cp -r ARCHDefs-2.0.3x ${prefix}/libexec/
cp -r SIFDecode-2.0.3 ${prefix}/libexec/
cp -r CUTEst-2.0.3 ${prefix}/libexec/
export ARCHDEFS=${prefix}/libexec/ARCHDefs-2.0.3x
export SIFDECODE=${prefix}/libexec/SIFDecode-2.0.3
export CUTEST=${prefix}/libexec/CUTEst-2.0.3

# build SIFDecode
cd $SIFDECODE

if [[ "${target}" == *-linux* || "${target}" == *-freebsd* ]]; then
  echo "6" > sifdecode.opts   # PC64
  echo "2" >> sifdecode.opts  # Linux
  echo "5" >> sifdecode.opts  # gfortran
elif [[ "${target}" == *-apple* ]]; then
  echo "13" > sifdecode.opts  # macOS
  echo "2" >> sifdecode.opts  # gfortran
elif [[ "${target}" == *-mingw* ]]; then
  cd $ARCHDEFS
  cd $SIFDECODE
  echo "6" > sifdecode.opts   # PC64
  echo "1" >> sifdecode.opts  # Windows
  echo "3" >> sifdecode.opts  # gfortran
fi
echo "nny" >> sifdecode.opts
./install_sifdecode < sifdecode.opts

# build CUTEst
cd $CUTEST
if [[ "${target}" == *-linux* || "${target}" == *-freebsd* ]]; then
  echo "6" > cutest.opts   # PC64
  echo "2" >> cutest.opts  # Linux
  echo "5" >> cutest.opts  # gfortran
  echo "2" >> cutest.opts  # build all tools except Matlab
  echo "7" >> cutest.opts  # gcc
  export MYARCH=pc64.lnx.gfo
elif [[ "${target}" == *-apple* ]]; then
  echo "13" > cutest.opts  # macOS
  echo "2" >> cutest.opts  # gfortran
  echo "2" >> cutest.opts  # build all tools except Matlab
  echo "5" >> cutest.opts  # gcc
  export MYARCH=mac64.osx.gfo
elif [[ "${target}" == *-mingw* ]]; then
  echo "5" > cutest.opts   # PC64
  echo "1" >> cutest.opts  # Windows
  echo "3" >> cutest.opts  # gfortran
  echo "2" >> cutest.opts  # build all tools except Matlab
  echo "5" >> cutest.opts  # gcc
  export MYARCH=pc64.mgw.gfo
fi
echo "nnydy" >> cutest.opts
./install_cutest < cutest.opts

# build shared libs
all_load="--whole-archive"
noall_load="--no-whole-archive"
extra=""
if [[ "${target}" == *-apple-* ]]; then
    all_load="-all_load"
    noall_load="-noall_load"
    extra="-Wl,-undefined -Wl,dynamic_lookup -headerpad_max_install_names"
fi
cd $CUTEST/objects/$MYARCH/double
gfortran -fPIC -shared ${extra} -Wl,${all_load} libcutest.a -Wl,${noall_load} -o libcutest_double.${dlext}
cd $CUTEST/objects/$MYARCH/single
gfortran -fPIC -shared ${extra} -Wl,${all_load} libcutest.a -Wl,${noall_load} -o libcutest_single.${dlext}

ln -s $SIFDECODE/bin/sifdecoder ${bindir}/
ln -s $SIFDECODE/objects/$MYARCH/double/slct ${bindir}/
ln -s $SIFDECODE/objects/$MYARCH/double/clsf ${bindir}/
install_license $CUTEST/lgpl-3.0.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# can't build shared libs on Windows, which imposes all symbols to be defined
platforms = expand_gfortran_versions(filter!(p -> !isa(p, Windows), supported_platforms()))

# The products that we will ensure are always built
products = [
    ExecutableProduct("sifdecoder", :sifdecoder),
    ExecutableProduct("slct", :slct),
    ExecutableProduct("clsf", :clsf),
]

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
home = homedir()
artifacts_dir = joinpath(homedir(), ".julia/artifacts")
mkpath(artifacts_dir)
datapath = joinpath(pwd(), "data")
write(
    joinpath(artifacts_dir, "Overrides.toml"),
    "113db057ea8fd9dd2ed6814748e832a38653cb43 = $(repr(datapath))"
)
import Pkg
Pkg.activate(".")
Pkg.develop(url="https://github.com/Wild-Area/VsRecorderBase.jl.git")
Pkg.develop(url="https://github.com/sunoru/SimpleI18n.jl.git")
Pkg.develop(url="https://github.com/sunoru/Tesseract.jl.git")
Pkg.activate("VsRecorderBase")
Pkg.develop("Tesseract")


text:
#RSDeltaSigma: simulateHBF() algorithm
#-------------------------------------------------------------------------------

"""`y = simulateHBF(x, f1, f2, mode=0)`

Simulate a Saramaki half-band filter.

The `f1` and `f2` vectors contain coefficients for the structure.
(`f1` and `f2` can also be struct arrays like those returned from designHBF.m
i.e. struct arrays whose .val fields contain the coeffiecients.)

The `mode` flag determines whether the input is filtered, interpolated,
or decimated according to the following table:
 - `mode = 0`: Plain filtering, no interpolation or decimation.
 - `mode = 1`: The input is interpolated.
 - `mode = 2`: The output is decimated, even samples are taken.
 - `mode = 3`: The output is decimated, odd samples are taken.
"""
function simulateHBF(x, f1, f2, mode::Int=0)
	if isa(f1, Number) && isnan(f1)
		(f1, f2) = exampleHBF(4)
	end
	if isa(f1, Vector{CoeffCSD})
		f1 = values(f1)
	end
	if isa(f2, Vector{CoeffCSD})
		f2 = values(f2)
	end
	x = x[:]
	f1 = f1[:]; f2 = f2[:] #Redundant unless `val`s are not always simple numbers
	n1 = length(f1); n2 = length(f2)
	f2imp = vcat(f2[end:-1:1], f2)
	if 0 == mode #Insert zeros
		f2imp = vcat(f2imp, zeros(2*n2))
		f2imp = f2imp[1:4*n2-1]
	end
	F2 = _tf(f2imp, vcat(1, zeros(length(f2imp)-1)), 1)

	if 0 == mode #Plain
		(up,) = _lsim(F2,x)
		y = 0.5*delay(x,2*n2-1) + f1[1]*up
		for i in 2:n1
			(up,) = _lsim(F2,up)
			(up,) = _lsim(F2,up)
			y = f1[i]*up + delay(y, 4*n2-2)
		end
	elseif 1 == mode #Interpolating
		up = zeros(size(x))
		nz = 2*n2-1
		for i in n1:-1:1
			if i==1
				(up,) = _lsim(F2,up+f1[i]*x)
				nz = n2-1
			else
				(up,) = _lsim(F2,up+f1[i]*x)
				(up,) = _lsim(F2,up)
			end
			x = delay(x,nz)
		end
		y = [2*up'; x']; y = y[:] #Interleave the upper and lower streams

	elseif mode in [2, 3] #Decimating
		x = x[1:2*div(end,2)]
		if 3 == mode
			y = 0.5*x[1:2:end]
			up = x[2:2:end]
			nz = n2-1
		else
			y = 0.5*x[2:2:end]
			up = x[1:2:end]
			nz = n2
		end
		for i in 1:n1
			if i==1
				(up,) = _lsim(F2,up)
			else
				(up,) = _lsim(F2,up)
				(up,) = _lsim(F2,up)
			end
			y = f1[i]*up + delay(y,nz)
			nz = 2*n2-1
		end
	else
		throw("$mode is not a valid value for the mode variable.")
	end
	return y
end

#Last line


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 3636e5d4-0276-11eb-2346-b77b042e1240
using Yao, YaoPlots

# ╔═╡ 1aba2418-0273-11eb-2fb6-a12d0afc33c6
md"# Applications
We'll discuss a few applications of qubits and quantum circuits. Namely the Bell States, Superdense Coding and Quantum Teleportation."

# ╔═╡ 84abe9cc-0273-11eb-0e38-c3a2b42c0848
md"## Bell States
The states,

$ \frac{|00> \;+\; |11>}{\sqrt2} ,\;\;$ $ \frac{|01> \;+\; |10>}{\sqrt2} \;\;$

$ \frac{|00> \;-\; |11>}{\sqrt2} \;$ and $ \frac{|01> \;-\; |10>}{\sqrt2} \;\;$

are known as the bell states. They are made by the bell circuit. The bell circuit looks like this."

# ╔═╡ 24727e9e-0276-11eb-344c-25a2ba5f138c
begin
	bellcircuit = chain(2, put(1=>H), control(1, 2=>X))
	plot(bellcircuit)
end

# ╔═╡ 6a117428-0276-11eb-156f-cf8597070a34
md"As you can see, the circuit takes in two qubits as input, and operates on them to give the bell states."

# ╔═╡ 932449a0-0277-11eb-0805-c533ef1ceebb
md"#### Feeding qubits to a circuit in Yao
There are many ways to create qubits to feed to quantum circuits in Yao."

# ╔═╡ 1a64822e-0278-11eb-1ed7-8f7f03623cee
q1 = ArrayReg(bit"00") #creating the system of two qubits with state |00>.

# ╔═╡ 6af54232-027d-11eb-087e-3b2516a14eb8
state(q1) #state of a qubit in vector form

# ╔═╡ 97a2930a-0278-11eb-1e5b-c1dc5b203022
md" **Note:** _**Normalizing** basically means making the summation of squares of probability amplitudes, equal to 1._"

# ╔═╡ 1e60af8c-0279-11eb-21fd-77cbe8640448
md"Other ways of doing it are"

# ╔═╡ 23f8e17c-0287-11eb-34ab-25db7b60da63
ArrayReg(bit"00") + ArrayReg(bit"11") |> normalize! #Equivalent to (1/√2)*(|00>+|11>)

# ╔═╡ 3951f92c-0279-11eb-1097-0b9935025360
zero_state(2) #2 qubits, both with the state |0> and |0>

# ╔═╡ 033aa410-0279-11eb-3742-c7b3084f847b
md" There! We have a system of two qubits! Let's try feeding the qubits to the Bell circuit we made!"

# ╔═╡ 4a6f0f16-027b-11eb-0d9e-77a3d5c6f1ee
a = (q1 |> bellcircuit) #Passing the qubit q1 through the bell circuit

# ╔═╡ f6d9ad26-0285-11eb-115c-db22262790db
state(a)

# ╔═╡ 7f471eb2-0286-11eb-3fa4-45952f77d4eb
md"### Reverse Bell Circuit
A circuit which reverses the effects of the bell circuit on a qubit. It's represented in a circuit as follows"

# ╔═╡ a6ab9b2e-0286-11eb-382e-07d449a188c6
begin
	reversebellcircuit = chain(2, control(1,2=>X), put(1=>H))
	plot(reversebellcircuit)
end

# ╔═╡ cfa3d672-0286-11eb-3ce5-47a7afe1d336
md"You can input the output state you got from the bell circuit, into the reverse bell circuit, and you'll get back your original state. Why not give it a try?"

# ╔═╡ 43241144-031c-11eb-01f4-d91cefa66359
md"Let's pass the qubits we got by passing two qubits in the Bell circuit, into the Reverse Bell Circuit."

# ╔═╡ 0f1b1ea0-031d-11eb-326c-65cd5aa41630
let
		res = (a |> reversebellcircuit)
		state(res)
end

# ╔═╡ f7b24140-0320-11eb-010d-db85da69ed7c
md"What do you think is the effect of single qubit gates on a qubit, which is entangled with another qubit? Lets check it out!"

# ╔═╡ 30a26136-0321-11eb-3375-af7b708a196d
begin
	singlequbitcircuit = chain(2, put(1=>X))
	plot(singlequbitcircuit)
end

# ╔═╡ 57be66c2-0326-11eb-0090-39e6b4f26b19
begin
	bellstate = ArrayReg(bit"00") + ArrayReg(bit"11") |> normalize!
	re = bellstate |> singlequbitcircuit
	state(re)
end

# ╔═╡ 68a743ea-0322-11eb-0039-d9a5ab8ad3f9
md"Can you notice how the circuit behaves as if its operating on one qubit? If the input is |{0}[1]> + |{1}[1]>, it behaves like the bits in the curly braces are the state of one qubit, while the bits in sqare brackets are the state of the other qubit, and it operates on them accordingly?"

# ╔═╡ 4d4432a2-031d-11eb-28f1-f39a26b8b789
md"### Super-dense Coding and Quantum Teleportation"

# ╔═╡ 8012368e-031d-11eb-1aa6-f53e9ca5f88f
md"Suppose we have two entangled qubits in a bell state, represented by
$ \frac{|00> \;+\; |11>}{\sqrt2} .$ Alice and Bob are two friends. Alice gets one qubit and Bob gets the other one. Both of them travel far apart with their qubits, and not measuring their qubits so as to preserve the entangled state. The above information is the premise for both:
1. Super-dense coding
2. Quantum Teleportation"

# ╔═╡ 669e0b32-031e-11eb-231f-956bf10c9850
md"##### 1. Super-dense coding
Alice wants to send Bob two classical bits of information. That means, one of the four states 00, 01, 10, 11. How does she do this?
Lets say she achieves this by passing her qubits through one of the following gates, corresponding to the information she wants to send.

1. If she wants to send Bob 00, then she'll send her qubit to Bob, as it is
2. If she wants to send Bob 01, then she'll send her qubit to Bob, after passing it through X gate.
3. If she wants to send Bob 10, then she'll send her qubit to Bob, after passing it through Z gate.
4. If she wants to send Bob 11, then she'll send her qubit to Bob, after passing it through Y gate.

Bob will then run both the qubits through the reverse bell circuit and measure them and he'll get what Alice wanted to send him."

# ╔═╡ 65c3b24c-0324-11eb-11bb-8dfba3e35770
md"Measuring qubits in Yao can be done by using the measure function.

Assume you've the a qubit q and you want to measure it, the syntax is q |> r->measure(r, nshots=number_of_runs)
We can't determine the probability amplitudes of a qubit. But, we can run the measurement many times, and count the frequency of each bit combination it gives.

If you don't get that, just know the syntax and keep nshots=1000.

Now, lets try measuring the result of the previous assignment!"

# ╔═╡ 0e2e6f62-0325-11eb-1637-116783c77ea8
measuredqubits = re |> r->measure(r, nshots=1000)

# ╔═╡ 89e76d2e-04b4-11eb-0704-17b457103932
let
	using StatsBase
	hist = fit(Histogram, Int.(measuredqubits), 1:3)
	h = hist.weights
	"The number of times the measurement result was 01 is $(h[1]/10)% and the number of times the measurement result was 10 is $(h[2]/10)"
end

# ╔═╡ d623229e-04b3-11eb-37e4-65be46f8a1c2
md"Sometimes the measurement gives 01 and sometimes 10. The probability of getting a 01 is same as getting the probability of 10, which is $ (1/\sqrt2)^2 $ that is 0.5 ."

# ╔═╡ 42eda908-04b4-11eb-1d71-2b2df3ab0b82
md"Wanna see how we can test that?"

# ╔═╡ ef6e1276-03d8-11eb-08b6-1fc082d158f6
md"Implementing the superdense coding."

# ╔═╡ 124eca24-03d9-11eb-3bca-c188b5d09404
begin
	Alice_and_Bobs_entangled_qubits = ArrayReg(bit"00") + ArrayReg(bit"11") |> normalize!
	input = ["00" "01" "10" "11"][rand(1:4)] #Assuming Alice wants to send one of these inputs to Bob.
end

# ╔═╡ 3e4d90dc-03df-11eb-371e-130d4339b472
state(Alice_and_Bobs_entangled_qubits)

# ╔═╡ ae574c4c-03da-11eb-2eb0-b1d06f9d61de
begin
	if(input == "00")
		Alices_circuit = chain(2)
	elseif(input == "01")
		Alices_circuit = chain(2, put(1=>X))
	elseif(input == "10")
		Alices_circuit = chain(2, put(1=>Z))
	elseif(input == "11")
		Alices_circuit = chain(2, put(1=>Y))
	end
	plot(Alices_circuit)
end

# ╔═╡ adf60706-03db-11eb-2206-c7f5ac69ab41
Bobs_part = ((Alice_and_Bobs_entangled_qubits |> Alices_circuit) |> reversebellcircuit) |> r->measure(r, nshots=1000)
#The content in the first roung bracket, outputs qubits, which are then fed to reversebellcircuit you made before, which is then fed to the measure function.

# ╔═╡ 8570a5d0-03b7-11eb-0f75-09aad5a4b811
md"##### 2. Quantum Teleportation
Alice now has one more qubit, in the state $ \;a|\;0> + \;b|\;1>.$ She wants to send her qubit to Bob, by changing the state of Bob's entangled qubit to Alice's extra qubit.

Confusing? Let's try naming the qubits. Alice and Bob have the qubits A and B respectively, which are entangled in the bell state. Alice has another qubit C.

Alice wants to send her qubit to Bob by changing the state of B, to that of C. She wants the state of the qubit B to be $ \;a|\;0> + \;b|\;1>. $ Alice doesn't know the values of a and b. Also, Bob can't make a measurement in any case as doing so will destroy B's state.

How does Alice do it? Well, she first passes both her qubits through the Reverse Bell circuit. She then measures both of her qubits. She gets one of the following outputs: $ |\;00>,\; |\;01>,\; |\;10>\;or\; |\;11>, $ each with probability 1/4. Alice sends these two classical bits to Bob, via any classical means, example, she texts him or calls and tells him. Bob knows that corresponding to her message, there are 4 options, and he takes action accordingly for each.
1. He gets the message, $ |\;00>, $ from Alice and knows that his qubit\(qubit B\) changed its state to that of qubit C.
2. He gets the message, $ |\;01>, $ and knows that his qubit changed its state to $ \;a|\;1> +\; b|\;0>, $ and he passes his qubit through X gate to change its state to that of qubit C.
3. He gets the message, $ |\;10>, $ and knows that his qubit changed its state to $ \;a|\;0> -\; b|\;1>, $ and he passes his qubit through Z gate to change its state to that of qubit C.
4. He gets the message, $ |\;11>, $ and knows that his qubit changed its state to $ \;a|\;1> -\; b|\;0>, $ and he passes his qubit through Y gate to change its state to that of qubit C."

# ╔═╡ 37d4357e-03be-11eb-3e1f-dda5a699af33
begin
	Alices_and_Bobs_entangled_qubits = ArrayReg(bit"00") + ArrayReg(bit"11") |> normalize!
	Alicequbit = rand_state(1) |> normalize!; #This function creates a qubit with a random state.
	state(Alicequbit)
end

# ╔═╡ 6ca5a9cc-03be-11eb-2afa-295c81c21e35
state(Alices_and_Bobs_entangled_qubits)

# ╔═╡ af31a570-03be-11eb-1af1-03950ee1c51f
begin
	teleportationcircuit = chain(3, control(1,2=>X), put(1=>H))
	plot(teleportationcircuit)
end

# ╔═╡ 6498549a-03bf-11eb-0348-a9072154b85e
begin
	feeding = join(Alices_and_Bobs_entangled_qubits, Alicequbit) |> teleportationcircuit
	state(feeding)
end

# ╔═╡ 14abac76-03c1-11eb-0634-afc86a739215
md"The join( ) function is used to join multiple qubits. It takes input from left to right."

# ╔═╡ b2e2fea0-03e4-11eb-1a20-157447a12ce5
Alices_measuredqubits = measure_remove!(feeding, 1:2)

# ╔═╡ 5331ad20-03e5-11eb-1244-1b494c4afa02
if(Alices_measuredqubits == bit"00")
	Bobs_qubit = feeding
elseif(Alices_measuredqubits == bit"01")
	Bobs_qubit = feeding |> chain(1, put(1=>Z))
elseif(Alices_measuredqubits == bit"10")
	Bobs_qubit = feeding |> chain(1, put(1=>X))
else
	Bobs_qubit = feeding |> chain(1, put(1=>Y))
end

# ╔═╡ 8e222eee-03e2-11eb-3419-edfe12b7c61a
md"The difference between **_r->measure(r, location of qubits(optional, measures all qubits if absent), nshots=number of runs)_** and *measure_remove!(qubit, the location of qubits(optional, measures all qubits if absent))* is that the latter collapses the state of the qubit, and removes the 0 + 0im rows of the state vector, while the former only collapses the state of the qubit."

# ╔═╡ 49604160-03c4-11eb-2480-ff7abbe90508
state(Bobs_qubit)

# ╔═╡ bec22f54-03e7-11eb-34fd-0fce9ea7da8f
md"Is Alice's qubit same as Bob's qubit now?
You can see that your for your self!"

# ╔═╡ 14a7320c-03e8-11eb-0909-1354fd574659
[state(Alicequbit) state(Bobs_qubit)]

# ╔═╡ 340a69fe-068c-11eb-3b37-ed2694843113
md"Left side : State of Alice's qubit. Right side : State of Bob's qubit. Almost equivalent!"

# ╔═╡ Cell order:
# ╠═3636e5d4-0276-11eb-2346-b77b042e1240
# ╟─1aba2418-0273-11eb-2fb6-a12d0afc33c6
# ╟─84abe9cc-0273-11eb-0e38-c3a2b42c0848
# ╠═24727e9e-0276-11eb-344c-25a2ba5f138c
# ╟─6a117428-0276-11eb-156f-cf8597070a34
# ╟─932449a0-0277-11eb-0805-c533ef1ceebb
# ╠═1a64822e-0278-11eb-1ed7-8f7f03623cee
# ╠═6af54232-027d-11eb-087e-3b2516a14eb8
# ╟─97a2930a-0278-11eb-1e5b-c1dc5b203022
# ╟─1e60af8c-0279-11eb-21fd-77cbe8640448
# ╠═23f8e17c-0287-11eb-34ab-25db7b60da63
# ╠═3951f92c-0279-11eb-1097-0b9935025360
# ╟─033aa410-0279-11eb-3742-c7b3084f847b
# ╠═4a6f0f16-027b-11eb-0d9e-77a3d5c6f1ee
# ╠═f6d9ad26-0285-11eb-115c-db22262790db
# ╟─7f471eb2-0286-11eb-3fa4-45952f77d4eb
# ╠═a6ab9b2e-0286-11eb-382e-07d449a188c6
# ╟─cfa3d672-0286-11eb-3ce5-47a7afe1d336
# ╟─43241144-031c-11eb-01f4-d91cefa66359
# ╠═0f1b1ea0-031d-11eb-326c-65cd5aa41630
# ╟─f7b24140-0320-11eb-010d-db85da69ed7c
# ╠═30a26136-0321-11eb-3375-af7b708a196d
# ╠═57be66c2-0326-11eb-0090-39e6b4f26b19
# ╟─68a743ea-0322-11eb-0039-d9a5ab8ad3f9
# ╟─4d4432a2-031d-11eb-28f1-f39a26b8b789
# ╟─8012368e-031d-11eb-1aa6-f53e9ca5f88f
# ╟─669e0b32-031e-11eb-231f-956bf10c9850
# ╟─65c3b24c-0324-11eb-11bb-8dfba3e35770
# ╠═0e2e6f62-0325-11eb-1637-116783c77ea8
# ╟─d623229e-04b3-11eb-37e4-65be46f8a1c2
# ╟─42eda908-04b4-11eb-1d71-2b2df3ab0b82
# ╠═89e76d2e-04b4-11eb-0704-17b457103932
# ╟─ef6e1276-03d8-11eb-08b6-1fc082d158f6
# ╠═124eca24-03d9-11eb-3bca-c188b5d09404
# ╠═3e4d90dc-03df-11eb-371e-130d4339b472
# ╠═ae574c4c-03da-11eb-2eb0-b1d06f9d61de
# ╠═adf60706-03db-11eb-2206-c7f5ac69ab41
# ╟─8570a5d0-03b7-11eb-0f75-09aad5a4b811
# ╠═37d4357e-03be-11eb-3e1f-dda5a699af33
# ╠═6ca5a9cc-03be-11eb-2afa-295c81c21e35
# ╠═af31a570-03be-11eb-1af1-03950ee1c51f
# ╠═6498549a-03bf-11eb-0348-a9072154b85e
# ╟─14abac76-03c1-11eb-0634-afc86a739215
# ╠═b2e2fea0-03e4-11eb-1a20-157447a12ce5
# ╠═5331ad20-03e5-11eb-1244-1b494c4afa02
# ╟─8e222eee-03e2-11eb-3419-edfe12b7c61a
# ╠═49604160-03c4-11eb-2480-ff7abbe90508
# ╟─bec22f54-03e7-11eb-34fd-0fce9ea7da8f
# ╠═14a7320c-03e8-11eb-0909-1354fd574659
# ╟─340a69fe-068c-11eb-3b37-ed2694843113

max_stars_repo_name: logankilpatrick/tutorials

max_issues_repo_name: logankilpatrick/tutorials

max_forks_repo_name: logankilpatrick/tutorials


text:
abstract type AbstractAmplitude end

"""
$(SIGNATURES)

Simulates the amplitude of interference signal which is
static over time.
"""
function propagate(jnr::Unitful.Gain, n0, sample_frequency, Δt, rng)
    jnr
end

@inline function get_amplitude(jnr::Unitful.Gain, n0, sample_frequency)
    sqrt(uconvertp(NoUnits, jnr) * n0 * sample_frequency)
end

@inline function get_jammer_to_noise_ratio(jnr::Unitful.Gain)
    jnr
end

struct DynamicJNR{T <: Real, VC <: AbstractVector{<: Unitful.Gain}, F <: Function} <: AbstractAmplitude
    jnrs::VC
    ampl::T
    sample_freq::typeof(1.0Hz)
    time::typeof(1.0s)
    ampl_filter::F
end

function DynamicJNR(jnrs, sample_freq; ampl_filter = x -> x, time = 0.0s)
    DynamicJNR(jnrs, 0.0, float(sample_freq), time, ampl_filter)
end

function propagate(dynamic_jnr::DynamicJNR, n0, sample_frequency, Δt, rng)
    next_time = dynamic_jnr.time + Δt
    jnr = get_sampled_value(next_time, dynamic_jnr.sample_freq, dynamic_jnr.jnrs)
    ampl = get_amplitude(jnr, n0, sample_frequency)
    ampl_filtered = dynamic_jnr.ampl_filter(ampl)
    DynamicJNR(dynamic_jnr.jnrs, ampl_filtered, dynamic_jnr.sample_freq, next_time, dynamic_jnr.ampl_filter)
end

@inline function get_amplitude(dynamic_jnr::DynamicJNR, n0, sample_frequency)
    dynamic_jnr.ampl
end

@inline function get_amplitude(dynamic_jnr::DynamicJNR)
    dynamic_jnr.ampl
end

@inline function get_jammer_to_noise_ratio(dynamic_jnr::DynamicJNR)
    dynamic_jnr
end

@inline get_amplitude(cn0::Unitful.Level, n0) = calc_amplitude_from_cn0(cn0, n0)

struct DynamicCN0{T <: Real, VC <: AbstractVector{<: Unitful.Level}, F <: Function} <: AbstractAmplitude
    cn0s::VC
    ampl::T
    sample_freq::typeof(1.0Hz)
    time::typeof(1.0s)
    ampl_filter::F
end

function DynamicCN0(cn0s, sample_freq; ampl_filter = x -> x, time = 0.0s)
    DynamicCN0(cn0s, 0.0, float(sample_freq), time, ampl_filter)
end

function propagate(dynamic_cn0::DynamicCN0, n0, Δt)
    next_time = dynamic_cn0.time + Δt
    cn0 = get_sampled_value(next_time, dynamic_cn0.sample_freq, dynamic_cn0.cn0s)
    ampl = get_amplitude(cn0, n0)
    ampl_filtered = dynamic_cn0.ampl_filter(ampl)
    DynamicCN0(dynamic_cn0.cn0s, ampl_filtered, dynamic_cn0.sample_freq, next_time, dynamic_cn0.ampl_filter)
end

@inline function get_amplitude(dynamic_cn0::DynamicCN0, n0)
    dynamic_cn0.ampl
end

@inline function get_amplitude(dynamic_cn0::DynamicCN0)
    dynamic_cn0.ampl
end

@inline function get_jammer_to_noise_ratio(dynamic_cn0::DynamicCN0)
    dynamic_cn0
end


text:
include("day9_solution.jl")
using Test

@test


text:
using ITensors,
      LinearAlgebra, Test

@testset "ITensor Decompositions" begin

  @testset "truncate!" begin
    a = [-0.1, -0.12]
    @test NDTensors.truncate!(a) == (0., 0.)
    @test length(a) == 1
    a = [0.1, 0.01, 1e-13]
    @test NDTensors.truncate!(a,
                             use_absolute_cutoff=true,
                             cutoff=1e-5) == (1e-13, (0.01 + 1e-13)/2)
    @test length(a) == 2
  end

  @testset "factorize" begin
    i = Index(2,"i")
    j = Index(2,"j")
    A = randomITensor(i,j)
    @test_throws ErrorException factorize(A, i, dir="left")
    @test_throws ErrorException factorize(A, i, ortho="fakedir")
  end

  @testset "factorize with eigen_perturbation" begin
    l = Index(4,"l")
    s1 = Index(2,"s1")
    s2 = Index(2,"s2")
    r = Index(4,"r")

    phi = randomITensor(l,s1,s2,r)

    drho = randomITensor(l',s1',l,s1)
    drho += swapprime(drho,0,1)
    drho .*= 1E-5

    U,B = factorize(phi,(l,s1),ortho="left",eigen_perturbation=drho)
    @test norm(U*B - phi) < 1E-5

    # Not allowed to use eigen_perturbation with which_decomp
    # other than "automatic" or "eigen":
    @test_throws ErrorException factorize(phi,(l,s1),
                                          ortho="left",
                                          eigen_perturbation=drho,
                                          which_decomp="svd")
  end

  @testset "factorize with QR" begin
    l = Index(5,"l")
    s = Index(2,"s")
    r = Index(10,"r")
    A = randomITensor(l,s,r)
    Q,R, = factorize(A,l,s; which_decomp="qr")
    q = commonind(Q,R)
    @test A ≈ Q*R atol=1e-13
    @test Q*dag(prime(Q,q)) ≈ δ(Float64,q,q') atol=1e-13

    R,Q, = factorize(A,l,s; which_decomp="qr", ortho="right")
    q = commonind(Q,R)
    @test A ≈ Q*R atol=1e-13
    @test Q*dag(prime(Q,q)) ≈ δ(Float64,q,q') atol=1e-13
  end

  @testset "eigen" begin
    i = Index(2,"i")
    j = Index(2,"j")
    A = randomITensor(i,i')
    eigA = eigen(A)
    Dt, Ut = eigen(NDTensors.tensor(A))
    eigArr = eigen(array(A))
    @test diag(array(eigA.D), 0) ≈ eigArr.values
    @test diag(array(Dt), 0) == eigArr.values
  end

  @testset "exp function" begin
    At = rand(10, 10)
    k  = Index(10, "k")
    A = itensor(At + transpose(At), k, k')
    @test Array(exp(Hermitian(NDTensors.tensor(A)))) ≈ exp(At + transpose(At))
  end

  @testset "Spectrum" begin
    i = Index(100,"i")
    j = Index(100,"j")

    U,S,V = svd(rand(100,100))
    S ./= norm(S)
    A = itensor(U*ITensors.diagm(0=>S)*V', i,j)

    spec = svd(A,i).spec

    @test eigs(spec) ≈ S .^2
    @test truncerror(spec) == 0.0

    spec = svd(A,i; maxdim=length(S)-3).spec
    @test truncerror(spec) ≈ sum(S[end-2:end].^2)

    @test entropy(Spectrum([0.5; 0.5], 0.0)) == log(2)
    @test entropy(Spectrum([1.0], 0.0)) == 0.0
    @test entropy(Spectrum([0.0], 0.0)) == 0.0

    @test isnothing(eigs(Spectrum(nothing, 1.0)))
    @test_throws ErrorException entropy(Spectrum(nothing, 1.0))
    @test truncerror(Spectrum(nothing, 1.0)) == 1.0
  end

  @testset "Eigen QN flux regression test" begin
    cutoff = 1E-12
    N = 4
    s = siteinds("S=1",N;conserve_qns=true)
    A = randomITensor(QN("Sz",2),s[1],s[2],s[3])

    R = A*dag(prime(A,s[1],s[2]))
    F = eigen(R,(s[1],s[2]),(s[1]',s[2]'))

    @test flux(F.Vt)==QN("Sz",0)
  end
end

nothing


text:
const InfiniteArrMissing = Union{
    AbstractArray{<:Union{Missing,Continuous}},
    AbstractArray{<:Union{Missing,Count}}}

# -----------------------------------------------------------
# MeanAbsoluteError

struct MeanAbsoluteError <: Aggregated end

metadata_measure(MeanAbsoluteError;
                 instances = ["mae", "mav", "mean_absolute_error",
                              "mean_absolute_value"],
                 target_scitype           = InfiniteArrMissing,
                 prediction_type          = :deterministic,
                 orientation              = :loss),

const MAE = MeanAbsoluteError
const MAV = MeanAbsoluteError
@create_aliases MeanAbsoluteError

@create_docs(MeanAbsoluteError,
body=
"""
``\\text{mean absolute error} =  n^{-1}∑ᵢ|yᵢ-ŷᵢ|`` or
``\\text{mean absolute error} = n^{-1}∑ᵢwᵢ|yᵢ-ŷᵢ|``
""",
scitype=DOC_INFINITE)

call(::MeanAbsoluteError, ŷ::ArrMissing{<:Real}, y::ArrMissing{<:Real}) =
    abs.(ŷ .- y) |> skipinvalid |> mean

call(::MeanAbsoluteError,
     ŷ::ArrMissing{<:Real},
     y::ArrMissing{<:Real},
     w::Arr{<:Real}) =
         abs.(ŷ .- y) .* w |> skipinvalid |> mean

# ----------------------------------------------------------------
# RootMeanSquaredError

struct RootMeanSquaredError <: Aggregated end

metadata_measure(RootMeanSquaredError;
                 instances                = ["rms", "rmse",
                                             "root_mean_squared_error"],
                 target_scitype           = InfiniteArrMissing,
                 prediction_type          = :deterministic,
                 orientation              = :loss,
                 aggregation              = RootMeanSquare())

const RMS = RootMeanSquaredError
@create_aliases RootMeanSquaredError

@create_docs(RootMeanSquaredError,
body=
"""
``\\text{root mean squared error} = \\sqrt{n^{-1}∑ᵢ|yᵢ-ŷᵢ|^2}`` or
``\\text{root mean squared error} = \\sqrt{\\frac{∑ᵢwᵢ|yᵢ-ŷᵢ|^2}{∑ᵢwᵢ}}``
""",
scitype=DOC_INFINITE)

call(::RootMeanSquaredError, ŷ::ArrMissing{<:Real}, y::ArrMissing{<:Real}) =
    (y .- ŷ).^2 |> skipinvalid |> mean |> sqrt

call(::RootMeanSquaredError,
     ŷ::ArrMissing{<:Real},
     y::ArrMissing{<:Real},
     w::Arr{<:Real}) =
         (y .- ŷ).^2 .* w |> skipinvalid |> mean |> sqrt

# -------------------------------------------------------------------------
# R-squared (coefficient of determination)

struct RSquared <: Aggregated end

metadata_measure(RSquared;
                 instances               = ["rsq", "rsquared"],
                 target_scitype          = InfiniteArrMissing,
                 prediction_type         = :deterministic,
                 orientation             = :score,
                 supports_weights        = false)

const RSQ = RSquared
@create_aliases RSquared

@create_docs(RSquared,
body=
"""
The R² (also known as R-squared or coefficient of determination) is suitable for interpreting linear regression analysis (Chicco et al., [2021](https://doi.org/10.7717/peerj-cs.623)).

Let ``\\overline{y}`` denote the mean of ``y``, then

``\\text{R^2} = 1 - \\frac{∑ (\\hat{y} - y)^2}{∑ \\overline{y} - y)^2}.``
""",
scitype=DOC_INFINITE)

function call(::RSquared, ŷ::ArrMissing{<:Real}, y::ArrMissing{<:Real})
    num = (ŷ .- y).^2 |> skipinvalid |> sum
    mean_y = mean(y)
    denom = (mean_y .- y).^2 |> skipinvalid |> sum
    return 1 - (num / denom)
end

# -------------------------------------------------------------------
# LP

struct LPLoss{T<:Real} <: Unaggregated
    p::T
end

LPLoss(; p=2.0) = LPLoss(p)

metadata_measure(LPLoss;
                 instances = ["l1", "l2"],
                 target_scitype           = InfiniteArrMissing,
                 prediction_type          = :deterministic,
                 orientation              = :loss)

const l1 = LPLoss(1)
const l2 = LPLoss(2)

@create_docs(LPLoss,
body=
"""
Constructor signature: `LPLoss(p=2)`. Reports
`|ŷ[i] - y[i]|^p` for every index `i`.
""",
scitype=DOC_INFINITE)

single(m::LPLoss, ŷ::Real, y::Real) =  abs(y - ŷ)^(m.p)

# ----------------------------------------------------------------------------
# RootMeanSquaredLogError

struct RootMeanSquaredLogError <: Aggregated end

metadata_measure(RootMeanSquaredLogError;
                 instances = ["rmsl", "rmsle", "root_mean_squared_log_error"],
                 target_scitype           = InfiniteArrMissing,
                 prediction_type          = :deterministic,
                 orientation              = :loss,
                 aggregation              = RootMeanSquare())

const RMSL = RootMeanSquaredLogError
@create_aliases RootMeanSquaredLogError

@create_docs(RootMeanSquaredLogError,
body=
"""
``\\text{root mean squared log error} =
n^{-1}∑ᵢ\\log\\left({yᵢ \\over ŷᵢ}\\right)``
""",
footer="See also [`rmslp1`](@ref).",
scitype=DOC_INFINITE)

call(::RootMeanSquaredLogError, ŷ::ArrMissing{<:Real}, y::ArrMissing{<:Real}) =
    (log.(y) - log.(ŷ)).^2 |> skipinvalid |> mean |> sqrt

call(::RootMeanSquaredLogError,
      ŷ::ArrMissing{<:Real},
      y::ArrMissing{<:Real},
      w::Arr{<:Real}) =
          (log.(y) - log.(ŷ)).^2 .* w |> skipinvalid |> mean |> sqrt

# ---------------------------------------------------------------------------
#  RootMeanSquaredLogProportionalError

struct RootMeanSquaredLogProportionalError{T<:Real} <: Aggregated
    offset::T
end

RootMeanSquaredLogProportionalError(; offset=1.0) =
    RootMeanSquaredLogProportionalError(offset)

metadata_measure(RootMeanSquaredLogProportionalError;
                 instances                = ["rmslp1", ],
                 target_scitype           = InfiniteArrMissing,
                 prediction_type          = :deterministic,
                 orientation              = :loss,
                 aggregation              = RootMeanSquare())

const RMSLP = RootMeanSquaredLogProportionalError
@create_aliases RootMeanSquaredLogProportionalError

@create_docs(RootMeanSquaredLogProportionalError,
body=
"""
Constructor signature: `RootMeanSquaredLogProportionalError(; offset = 1.0)`.

``\\text{root mean squared log proportional error} =
n^{-1}∑ᵢ\\log\\left({yᵢ + \\text{offset} \\over ŷᵢ + \\text{offset}}\\right)``
""",
footer="See also [`rmsl`](@ref). ",
scitype=DOC_INFINITE)

call(m::RMSLP, ŷ::ArrMissing{<:Real}, y::ArrMissing{<:Real}) =
    (log.(y .+ m.offset) - log.(ŷ .+ m.offset)).^2 |>
    skipinvalid |> mean |> sqrt

call(m::RMSLP, ŷ::ArrMissing{<:Real}, y::ArrMissing{<:Real}, w::Arr{<:Real}) =
    (log.(y .+ m.offset) - log.(ŷ .+ m.offset)).^2 .* w |>
    skipinvalid |> mean |> sqrt

# --------------------------------------------------------------------------
# RootMeanSquaredProportionalError

struct RootMeanSquaredProportionalError{T<:Real} <: Aggregated
    tol::T
end

RootMeanSquaredProportionalError(; tol=eps()) =
    RootMeanSquaredProportionalError(tol)

metadata_measure(RootMeanSquaredProportionalError;
    instances                = ["rmsp", ],
    target_scitype           = InfiniteArrMissing,
    prediction_type          = :deterministic,
    orientation              = :loss,
    aggregation              = RootMeanSquare())

const RMSP = RootMeanSquaredProportionalError
@create_aliases RMSP

@create_docs(RootMeanSquaredProportionalError,
body=
"""
Constructor keyword arguments: `tol` (default = `eps()`).

``\\text{root mean squared proportional error} =
m^{-1}∑ᵢ \\left({yᵢ-ŷᵢ \\over yᵢ}\\right)^2``

where the sum is over indices such that `abs(yᵢ) > tol` and `m` is the number
of such indices.

""", scitype=DOC_INFINITE)

function call(m::RootMeanSquaredProportionalError,
               ŷ::ArrMissing{<:Real},
               y::ArrMissing{T},
               w::Union{Nothing,Arr{<:Real}}=nothing) where T <: Real
    ret = zero(T)
    count = 0
    @inbounds for i in eachindex(y)
        (isinvalid(y[i]) || isinvalid(ŷ[i])) && continue
        ayi = abs(y[i])
        if ayi > m.tol
            dev = ((y[i] - ŷ[i]) / ayi)^2
            ret += dev
            ret = _scale(ret, w, i)
            count += 1
        end
    end
    return sqrt(ret / count)
end

# -----------------------------------------------------------------------
# MeanAbsoluteProportionalError

struct MeanAbsoluteProportionalError{T} <: Aggregated
    tol::T
end

MeanAbsoluteProportionalError(; tol=eps()) = MeanAbsoluteProportionalError(tol)

metadata_measure(MeanAbsoluteProportionalError;
    instances                = ["mape", ],
    target_scitype           = InfiniteArrMissing,
    prediction_type          = :deterministic,
    orientation              = :loss)

const MAPE = MeanAbsoluteProportionalError
@create_aliases MAPE

@create_docs(MeanAbsoluteProportionalError,
body=
"""
Constructor key-word arguments: `tol` (default = `eps()`).

``\\text{mean absolute proportional error} =  m^{-1}∑ᵢ|{(yᵢ-ŷᵢ) \\over yᵢ}|``

where the sum is over indices such that `abs(yᵢ) > tol` and `m` is the number
of such indices.
""", scitype=DOC_INFINITE)

function call(m::MeanAbsoluteProportionalError,
              ŷ::ArrMissing{<:Real},
              y::ArrMissing{T},
              w::Union{Nothing,Arr{<:Real}}=nothing) where T <: Real
    ret = zero(T)
    count = 0
    @inbounds for i in eachindex(y)
        (isinvalid(y[i]) || isinvalid(ŷ[i])) && continue
        ayi = abs(y[i])
        if ayi > m.tol
        #if y[i] != zero(eltype(y))
            dev = abs((y[i] - ŷ[i]) / ayi)
            ret += dev
            ret =_scale(ret, w, i)
            count += 1
        end
    end
    return ret / count
end

# -------------------------------------------------------------------------
# LogCoshLoss

struct LogCoshLoss <: Unaggregated end

metadata_measure(LogCoshLoss;
    instances                = ["log_cosh", "log_cosh_loss"],
    target_scitype           = InfiniteArrMissing,
    prediction_type          = :deterministic,
    orientation              = :loss)

const LogCosh = LogCoshLoss
@create_aliases LogCoshLoss

@create_docs(LogCoshLoss,
             body="Reports ``\\log(\\cosh(ŷᵢ-yᵢ))`` for each index `i`. ",
             scitype=DOC_INFINITE)

_softplus(x::T) where T<:Real = x > zero(T) ? x + log1p(exp(-x)) : log1p(exp(x))
_log_cosh(x::T) where T<:Real = x + _softplus(-2x) - log(convert(T, 2))

single(::LogCoshLoss, ŷ::Real, y::Real) = _log_cosh(ŷ - y)


text:
using JSON
using Requests
import Requests: URI, post

type CoinigyHandler
  base_url::String
  api::Dict
  debug::Bool
  max_attempts::Int

  function CoinigyHandler(;
                    base_url = "https://api.coinigy.com/api/v1",
                    public_key = "",
                    private_key = "",
                    debug = false,
                    max_attempts = 6)
      api = Dict("apiKey" => public_key, "apiSecret" => private_key)
      endswith(base_url, "/") && (base_url = chop(base_url))
      new(base_url, api, debug, max_attempts)
  end
end

function getResult(handler::CoinigyHandler, url::String, data = Dict())
  headers = Dict("Content-Type" => "application/json",
  "X-API-KEY" => handler.api["apiKey"],
  "X-API-SECRET" => handler.api["apiSecret"])
  attempt = 1
  handler.debug && (@show data)
  resp = post(url; headers = headers, data = JSON.json(data))
  handler.debug && (@show resp)
  # Retry if service temporarily not available
  while resp.status == 503 && attempt < handler.max_attempts
    sleep(attempt) #increase wait time by one second for each failed attempt
    attempt += 1
    resp = post(url; headers = headers, data = JSON.json(data))
  end
  if resp.status == 503
    error("Aborting request with url: $url, data: $data after $(handler.max_attempts) attempts")
  elseif resp.status != 200
    error("$(resp.status): Error executing the request with url: $url, data: $data - $resp")
  else
    try
      parsedresp = Requests.json(resp)
      handler.debug && (@show parsedresp)
      if "error" in keys(parsedresp)
        println("Error parsing response to request with url: $url, data: $data - $(parsedresp["error"])")
      end
      # return Array{Dict}
      return parsedresp
    catch e
      error("Error parsing response to request with url: $url, data: $data - $resp - $parsedresp")
    end
  end
end

getResultData(handler::CoinigyHandler, url::String, data = Dict()) = map(d -> d, getResult(handler, url, data)["data"])

getUserInfo(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/userInfo")
getActivityLog(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/activity")
getNotifications(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/pushNotifications")
getAccounts(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/accounts")
getBalances(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/balances")
getBalanceHistory(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/balanceHistory")
getOrders(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/orders")
getOrderTypes(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/orderTypes")
getAlerts(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/alerts")
getUserWatchList(handler) = getResultData(handler, "$(handler.base_url)/userWatchList")
getNewsFeed(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/newsFeed")
getExchanges(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/exchanges")
getMarkets(handler::CoinigyHandler) = getResultData(handler, "$(handler.base_url)/markets")
getTradeHistory(handler::CoinigyHandler, exchange::String, market::String) = getResultData(handler, "$(handler.base_url)/data", Dict("exchange_code" => exchange, "exchange_market" => market, "type" => "history"))
getAsks(handler::CoinigyHandler, exchange::String, market::String)  = getResultData(handler, "$(handler.base_url)/data", Dict("exchange_code" => exchange, "exchange_market" => market, "type" => "asks"))
getBids(handler::CoinigyHandler, exchange::String, market::String)  = getResultData(handler, "$(handler.base_url)/data", Dict("exchange_code" => exchange, "exchange_market" => market, "type" => "bids"))
getOrders(handler::CoinigyHandler, exchange::String, market::String)  = getResultData(handler, "$(handler.base_url)/data", Dict("exchange_code" => exchange, "exchange_market" => market, "type" => "orders"))
getMarketData(handler::CoinigyHandler, exchange::String, market::String)  = getResultData(handler, "$(handler.base_url)/data", Dict("exchange_code" => exchange, "exchange_market" => market, "type" => "all"))
getTicker(handler::CoinigyHandler, exchange::String, market::String)  = getResultData(handler, "$(handler.base_url)/ticker", Dict("exchange_code" => exchange, "exchange_market" => market))
activateTradingKey(handler::CoinigyHandler, authId) = getResult(handler, "$(handler.base_url)/activateTradingKey", Dict("auth_id" => authId, "auth_trade" => 1))
addOrder(handler::CoinigyHandler, authId, exchId, mktId, orderTypeId, priceTypeId, price, quantity) = getResult(handler, "$(handler.base_url)/addOrder", Dict("auth_id" => authId, "exch_id" => exchId, "mkt_id" => mktId, "order_type_id" => orderTypeId, "price_type_id" => priceTypeId, "limit_price" => price, "order_quantity" => quantity))


text:
using Documenter, Stan

DOC_ROOT = joinpath(dirname(pathof(Stan)), "..", "docs")

makedocs( root = DOC_ROOT,
  modules = [Stan],
  sitename = "StanJulia/Stan.jl",
  authors = "Rob J Goedman",
  pages = Any[
      "Home" => "INTRO.md",
      "Installation" => "INSTALLATION.md",
      "Walkthrough" => "WALKTHROUGH.md",
      "Feature cross reference" => "EXAMPLES.md",
      "Versions" => "VERSIONS.md",
      "Index" => "index.md"
  ]
)

deploydocs(
  root = DOC_ROOT,
  repo = "github.com/StanJulia/Stan.jl.git",
 )


text:
# author: Ziyi (Francis) Yin
# This script trains a Fourier Neural Operator which maps 2D permeability distribution to time-varying CO2 concentration snapshots.
# The PDE is in 2D while FNO requires 3D FFT

using DrWatson
@quickactivate "FNO"
import Pkg; Pkg.instantiate()

using PyPlot
using BSON
using Flux, Random, FFTW, Zygote, NNlib
using MAT, Statistics, LinearAlgebra
using CUDA
using ProgressMeter

try
    CUDA.device()
    global gpu_flag=true
catch e
    println("CUDA.device() found no GPU device on this machine.")
    global gpu_flag=false
end

CUDA.culiteral_pow(::typeof(^), a::Complex{Float32}, b::Val{2}) = real(conj(a)*a)
CUDA.sqrt(a::Complex) = cu(sqrt(a))
Base.broadcasted(::typeof(sqrt), a::Base.Broadcast.Broadcasted) = Base.broadcast(sqrt, Base.materialize(a))

include("utils.jl")

Random.seed!(3)

mutable struct SpectralConv3d_fast{T,N}
    weights1::AbstractArray{T,N}
    weights2::AbstractArray{T,N}
    weights3::AbstractArray{T,N}
    weights4::AbstractArray{T,N}
end

@Flux.functor SpectralConv3d_fast

# Constructor
function SpectralConv3d_fast(in_channels::Integer, out_channels::Integer, modes1::Integer, modes2::Integer, modes3::Integer)
    scale = (1f0 / (in_channels * out_channels))
    if gpu_flag
        weights1 = scale*randn(Complex{Float32}, modes1, modes2, modes3, in_channels, out_channels) |> gpu
        weights2 = scale*randn(Complex{Float32}, modes1, modes2, modes3, in_channels, out_channels) |> gpu
        weights3 = scale*randn(Complex{Float32}, modes1, modes2, modes3, in_channels, out_channels) |> gpu
        weights4 = scale*randn(Complex{Float32}, modes1, modes2, modes3, in_channels, out_channels) |> gpu
    else
        weights1 = scale*rand(Complex{Float32}, modes1, modes2, modes3, in_channels, out_channels)
        weights2 = scale*rand(Complex{Float32}, modes1, modes2, modes3, in_channels, out_channels)
        weights3 = scale*rand(Complex{Float32}, modes1, modes2, modes3, in_channels, out_channels)
        weights4 = scale*rand(Complex{Float32}, modes1, modes2, modes3, in_channels, out_channels)
    end
    return SpectralConv3d_fast{Complex{Float32}, 5}(weights1, weights2, weights3, weights4)
end

function compl_mul3d(x::AbstractArray{Complex{Float32}}, y::AbstractArray{Complex{Float32}})
    # complex multiplication
    # x in (modes1, modes2, modes3, input channels, batchsize)
    # y in (modes1, modes2, modes3, input channels, output channels)
    # output in (modes1,modes2,modes3,output channels,batchsize)
    x_per = permutedims(x,[5,4,1,2,3]) # batchsize*in_channels*modes1*modes2*modes3
    y_per = permutedims(y,[4,5,1,2,3]) # in_channels*out_channels*modes1*modes2*modes3
    x_resh = reshape(x_per,size(x_per,1),size(x_per,2),:) # batchsize*in_channels*(modes1*modes2*modes3)
    y_resh = reshape(y_per,size(y_per,1),size(y_per,2),:) # in_channels*out_channels*(modes1*modes2*modes3)
    out_resh = batched_mul(x_resh,y_resh) # batchsize*out_channels*(modes1*modes2*modes3)
    out_per = reshape(out_resh,size(out_resh,1),size(out_resh,2),size(x,1),size(x,2),size(x,3)) # batchsize*out_channels*modes1*modes2*modes3
    out = permutedims(out_per,[3,4,5,2,1])
    return out
end

function (L::SpectralConv3d_fast)(x::AbstractArray{Float32})
    # x in (size_x, size_y, time, channels, batchsize)
    x_ft = rfft(x,[1,2,3])      ## full size FFT
    modes1 = size(L.weights1,1)
    modes2 = size(L.weights1,2)
    modes3 = size(L.weights1,3)
    ### only keep low frequency coefficients
    out_ft = cat(cat(cat(compl_mul3d(x_ft[1:modes1, 1:modes2, 1:modes3, :,:], L.weights1),
                0f0im .* view(x_ft, 1:modes1, 1:modes2, 1:size(x_ft,3)-2*modes3, :, :),
                compl_mul3d(x_ft[1:modes1, 1:modes2, end-modes3+1:end,:,:], L.weights2),dims=3),
                0f0im .* view(x_ft, 1:modes1, 1:size(x_ft, 2)-2*modes2, :, :, :),
                cat(compl_mul3d(x_ft[1:modes1, end-modes2+1:end, 1:modes3,:,:], L.weights3),
                0f0im .* view(x_ft, 1:modes1, 1:modes2, 1:size(x_ft,3)-2*modes3, :, :),
                compl_mul3d(x_ft[1:modes1, end-modes2+1:end, end-modes3+1:end,:,:], L.weights4),dims=3)
                ,dims=2),
                0f0im .* view(x_ft, 1:size(x_ft,1)-modes1, :, :, :, :),dims=1)
    out_ft = irfft(out_ft, size(x,1),[1,2,3])
end

mutable struct SimpleBlock3d
    fc0::Conv
    conv0::SpectralConv3d_fast
    conv1::SpectralConv3d_fast
    conv2::SpectralConv3d_fast
    conv3::SpectralConv3d_fast
    w0::Conv
    w1::Conv
    w2::Conv
    w3::Conv
    bn0::BatchNorm
    bn1::BatchNorm
    bn2::BatchNorm
    bn3::BatchNorm
    fc1::Conv
    fc2::Conv
end

@Flux.functor SimpleBlock3d

function SimpleBlock3d(modes1::Integer, modes2::Integer, modes3::Integer, width::Integer)
    block = SimpleBlock3d(
        Conv((1, 1, 1), 4=>width),
        SpectralConv3d_fast(width, width, modes1, modes2, modes3),
        SpectralConv3d_fast(width, width, modes1, modes2, modes3),
        SpectralConv3d_fast(width, width, modes1, modes2, modes3),
        SpectralConv3d_fast(width, width, modes1, modes2, modes3),
        Conv((1, 1, 1), width=>width),
        Conv((1, 1, 1), width=>width),
        Conv((1, 1, 1), width=>width),
        Conv((1, 1, 1), width=>width),
        BatchNorm(width, identity; ϵ=1.0f-5, momentum=.1f0),
        BatchNorm(width, identity; ϵ=1.0f-5, momentum=.1f0),
        BatchNorm(width, identity; ϵ=1.0f-5, momentum=.1f0),
        BatchNorm(width, identity; ϵ=1.0f-5, momentum=.1f0),
        Conv((1, 1, 1), width=>128),
        Conv((1, 1, 1), 128=>1)
    )
    return block
end

function (B::SimpleBlock3d)(x::AbstractArray{Float32})
    x = B.fc0(x)
    x1 = B.conv0(x)
    x2 = B.w0(x)
    x = B.bn0(x1+x2)
    x = relu.(x)
    x1 = B.conv1(x)
    x2 = B.w1(x)
    x = B.bn1(x1+x2)
    x = relu.(x)
    x1 = B.conv2(x)
    x2 = B.w2(x)
    x = B.bn2(x1+x2)
    x = relu.(x)
    x1 = B.conv3(x)
    x2 = B.w3(x)
    x = B.bn3(x1+x2)
    x = B.fc1(x)
    x = relu.(x)
    x = B.fc2(x)
    return x
end

mutable struct Net3d
    conv1::SimpleBlock3d
end

@Flux.functor Net3d

function Net3d(modes::Integer, width::Integer)
    return Net3d(SimpleBlock3d(modes,modes,modes,width))
end

function (NN::Net3d)(x::AbstractArray{Float32})
    x = NN.conv1(x)
    x = dropdims(x,dims=4)
end

ntrain = 1000
ntest = 100

batch_size = 10
learning_rate = 1f-4

epochs = 200
step_size = 100
gamma = 5f-1

modes = 4
width = 20

n = (64,64)
#d = (15f0,15f0) # dx, dy in m
d = (1f0/64, 1f0/64)

s = 4

nt = 51
#dt = 20f0    # dt in day
dt = 1f0/nt

# Define raw data directory
mkpath(datadir("data"))
perm_path = datadir("data", "perm.mat")
conc_path = datadir("data", "conc.mat")

# Download the dataset into the data directory if it does not exist
if isfile(perm_path) == false
    run(`wget https://www.dropbox.com/s/xzicjq9fnessdif/'
        'perm.mat -q -O $perm_path`)
end
if isfile(conc_path) == false
    run(`wget https://www.dropbox.com/s/s7ph9gf2xwlu5mb/'
        'conc.mat -q -O $conc_path`)
end

perm = matread(perm_path)["perm"]
conc = matread(conc_path)["conc"]

x_train_ = convert(Array{Float32},perm[1:s:end,1:s:end,1:ntrain])
x_test_ = convert(Array{Float32},perm[1:s:end,1:s:end,end-ntest+1:end])

y_train_ = convert(Array{Float32},conc[:,1:s:end,1:s:end,1:ntrain])
y_test_ = convert(Array{Float32},conc[:,1:s:end,1:s:end,end-ntest+1:end])

y_train_ = permutedims(y_train_,[2,3,1,4])
y_test = permutedims(y_test_,[2,3,1,4])

x_normalizer = UnitGaussianNormalizer(x_train_)
x_train_ = encode(x_normalizer,x_train_)
x_test_ = encode(x_normalizer,x_test_)

y_normalizer = UnitGaussianNormalizer(y_train_)
y_train = encode(y_normalizer,y_train_)

x = reshape(collect(range(d[1],stop=n[1]*d[1],length=n[1])), :, 1)
z = reshape(collect(range(d[2],stop=n[2]*d[2],length=n[2])), 1, :)

grid = zeros(Float32,n[1],n[2],2)
grid[:,:,1] = repeat(x',n[2])'
grid[:,:,2] = repeat(z,n[1])

x_train = zeros(Float32,n[1],n[2],nt,4,ntrain)
x_test = zeros(Float32,n[1],n[2],nt,4,ntest)

for i = 1:nt
    x_train[:,:,i,1,:] = deepcopy(x_train_)
    x_test[:,:,i,1,:] = deepcopy(x_test_)
    for j = 1:ntrain
        x_train[:,:,i,2,j] = grid[:,:,1]
        x_train[:,:,i,3,j] = grid[:,:,2]
        x_train[:,:,i,4,j] .= (i-1)*dt
    end

    for k = 1:ntest
        x_test[:,:,i,2,k] = grid[:,:,1]
        x_test[:,:,i,3,k] = grid[:,:,2]
        x_test[:,:,i,4,k] .= (i-1)*dt
    end
end

# value, x, y, t

train_loader = Flux.Data.DataLoader((x_train, y_train); batchsize = batch_size, shuffle = true)
test_loader = Flux.Data.DataLoader((x_test, y_test); batchsize = batch_size, shuffle = false)

if gpu_flag
    y_normalizer.mean_ = y_normalizer.mean_ |> gpu
    y_normalizer.std_ = y_normalizer.std_   |> gpu
    y_normalizer.eps_ = y_normalizer.eps_   |> gpu
    NN = Net3d(modes, width) |> gpu
else
    NN = Net3d(modes, width)
end

w = Flux.params(NN)
Flux.trainmode!(NN, true)
opt = Flux.Optimise.ADAMW(learning_rate, (0.9f0, 0.999f0), 1f-4)

Loss = zeros(Float32,epochs*Int(ntrain/batch_size))

prog = Progress(ntrain * epochs)

iter = 0
for ep = 1:epochs
    Base.flush(Base.stdout)
    for (x,y) in train_loader
        global iter = iter + 1
        grads = gradient(w) do
            if gpu_flag
                x = x |> gpu
                y = y |> gpu
            end
            out = decode(y_normalizer,NN(x))
            y_n = decode(y_normalizer,y)
            global loss = Flux.mse(out,y_n;agg=sum)
            return loss
        end
        Loss[iter] = loss
        for p in w
            Flux.Optimise.update!(opt, p, grads[p])
        end
        ProgressMeter.next!(prog; showvalues = [(:loss, loss), (:epoch, ep)])
    end
end

if gpu_flag
    NN = NN |> cpu
    w = convert.(Array,w) |> cpu
end

# Define result directory
mkpath(datadir("TrainedNet"))
BSON.@save "data/TrainedNet/2phasenet_$epochs.bson" NN w batch_size Loss modes width learning_rate epochs gamma step_size s n d nt dt


text:
using Base.Test, Jewel

@test Jewel.getthing("Base.fft") == fft
@test Jewel.getthing(Base, [:fft]) == fft

@test Jewel.filemodule(Pkg.dir("Jewel", "src", "module.jl")) == "Jewel"

include("utils.jl")
include("scope.jl")
include("display.jl")


text:
module Comm

using Sockets
using JSON

export send

const PORT = 2001

mutable struct SocClient
    socket::TCPSocket
    chan::Channel{String}

    function SocClient()
        o = new()

        try
            o.socket = connect(PORT)
            o.chan = Channel{String}(10)
        catch ex
            return nothing
        end

        o
    end

end

# The task dealing with the socket to server com.
function listen(soc::SocClient)
    println("Start listening to socket.")

    # **Socket Task**
    @async while true
        if isopen(soc.socket)
            # println("Socket is open. Reading line...")
            data = readline(soc.socket)
            msg = String(data)
            # println("Msg from server socket: ", msg)
            # println("Routing to channel")
            put!(soc.chan, msg)
            # println("Routing done")
        else
            println("Socket closed!")
            break;
        end
    end

    println("Listening...")
end

# Called from Gui loop
function read_channel(soc::SocClient)
    if isready(soc.chan)
        # println("Taking...")
        json = take!(soc.chan)
        # println("Msg from Channel: ", json)

        if json ≠ ""
            data = JSON.parse(json)

            response = Main.handle_msg(data)

            if response ≠ nothing
                println(soc.socket, JSON.json(response))
            end
        end
    else
        # The **Socket Task** above needs time to run so we yield.
        yield()
        # No need to sleep because this function is called in the Gui loop
    end
end

function send(soc::SocClient, data::Dict{String,Any})
    # Put on channel for task to take!
    chan_msg = JSON.json(data)
    put!(soc.chan, chan_msg)
    # println("Put on channel: ", chan_msg)
end

end # Module ------------------------


text:
# ---
# title: 437. Path Sum III
# id: problem437
# author: Indigo
# date: 2021-06-21
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/path-sum-iii/description/>
# hidden: true
# ---
#
# You are given a binary tree in which each node contains an integer value.
#
# Find the number of paths that sum to a given value.
#
# The path does not need to start or end at the root or a leaf, but it must go
# downwards (traveling only from parent nodes to child nodes).
#
# The tree has no more than 1,000 nodes and the values are in the range
# -1,000,000 to 1,000,000.
#
# **Example:**
#
#
#
#     root = [10,5,-3,3,2,null,11,3,-2,null,1], sum = 8
#
#           10
#          /  \
#         **5**   **-3**
#        **/** **\**    **\**
#       **3**   **2**   **11**
#      / \   **\**
#     3  -2   **1**
#
#     Return 3. The paths that sum to 8 are:
#
#     1.  5 -> 3
#     2.  5 -> 2 -> 1
#     3. -3 -> 11
#
#
#
## @lc code=start
using LeetCode

function tree_path_sum(root::TreeNode{Int}, sm::Int)::Int
    res = 0
    tsum(::Nothing, sm::Int) = nothing
    function tsum(root::TreeNode{Int}, sm::Int)::Nothing
        sm -= root.val
        sm == 0 && (res += 1)
        tsum(root.left, sm)
        tsum(root.right, sm)
    end
    tsum(root, sm)
    res += tree_path_sum(root.left, sm)
    res += tree_path_sum(root.right, sm)
end

tree_path_sum(root::Nothing, sum::Int) = 0
## @lc code=end


text:
const DEFAULT_REGISTRY = "General"

const skip_lists = Dict{String,Vector{String}}()
const retry_lists = Dict{String,Vector{String}}()

read_registries() = TOML.parsefile(registries_file())

"""
    prepare_registry([name])

Prepare the given registry for use by this package and update it if `update` is true.
`name` must correspond to an existing stanza in the `deps/Registries.toml` file.
"""
function prepare_registry(name=DEFAULT_REGISTRY; update::Bool=false)
    reg = read_registries()[name]
    regspec = RegistrySpec(name = name, url = reg["url"], uuid = UUID(reg["uuid"]))

    # clone and update the registry
    if !any(existing_regspec -> existing_regspec.name == name,
            Pkg.Registry.reachable_registries())
        Pkg.Types.clone_or_cp_registries([regspec])
    elseif update
        Pkg.Registry.update(name)
    end

    # read some metadata
    skip_lists[name] = get(reg, "skip", String[])
    retry_lists[name] = get(reg, "retry", String[])

    return
end

"""
    read_pkgs([pkgs::Vector{String}]; [registry::String])

Read all packages from a registry and return them as a vector of tuples containing the
package name and registry, its UUID, and a path to it. If `pkgs` is given, only collect
packages matching the names in `pkgs`
"""
function read_pkgs(pkg_names::Vector{String}=String[]; registry=DEFAULT_REGISTRY)
    pkg_names = Set(pkg_names)
    want_all = isempty(pkg_names)

    pkg_data = []
    registry_instance = only(filter(ri->ri.name == "General",
                             Pkg.Registry.reachable_registries()))
    for (uuid, pkg) in registry_instance
        if !want_all
            pkg.name in pkg_names || continue
            delete!(pkg_names, pkg.name)
        end

        push!(pkg_data, pkg)
    end

    if !want_all && !isempty(pkg_names)
        @warn """did not find the following packages in the $registry registry:\n $("  - " .* join(pkg_names, '\n'))"""
    end

    return pkg_data
end


text:
module Webservice
using HTTP, JSON, UUIDs, Distributed, Sockets
using ..Data, ..Models, ..Utils
import ..Utils: log_message, log_progress, PROJECT_ROOT

include("utils.jl")
include("request_handlers.jl")
include("sample_pipeline.jl")
include("server.jl")
include("routes.jl")

export start_server, stop_server
end


text:
@testset "$(rpad("Lagrange Basis Tests",80))" begin

    b = LagrangeBasis([0.0, 1.0])

    @test evaluate(b, 1, 0.0) == 1.0
    @test evaluate(b, 1, 0.5) == 0.5
    @test evaluate(b, 1, 1.0) == 0.0

    @test evaluate(b, 2, 0.0) == 0.0
    @test evaluate(b, 2, 0.5) == 0.5
    @test evaluate(b, 2, 1.0) == 1.0

    @test derivative(b, 1, 0.0) == -1.0
    @test derivative(b, 1, 0.5) == -1.0
    @test derivative(b, 1, 1.0) == -1.0

    @test derivative(b, 2, 0.0) == 1.0
    @test derivative(b, 2, 0.5) == 1.0
    @test derivative(b, 2, 1.0) == 1.0

    @test derivative(b, 1, 1) == -1.0
    @test derivative(b, 1, 2) == -1.0

    @test derivative(b, 2, 1) == 1.0
    @test derivative(b, 2, 2) == 1.0

    @test integral(b, 1, 0.0) == 0.0
    @test integral(b, 1, 0.5) == 0.375
    @test integral(b, 1, 1.0) == 0.5

    @test integral(b, 2, 0.0) == 0.0
    @test integral(b, 2, 0.5) == 0.125
    @test integral(b, 2, 1.0) == 0.5

    @test integral(b, 1, 1) == 0.0
    @test integral(b, 1, 2) == 0.5

    @test integral(b, 2, 1) == 0.0
    @test integral(b, 2, 2) == 0.5

end


text:
# NOTE: Luxor is using JuliaGraphics/Luxor.jl/pull/48
using Luxor, BitBasis

"""
   btable(n, center=O; len=ndigits(n), width=20, title=false)

Draw a table of binary digits with `n` which can be [`BitStr`](@ref) or an `Integer`.
"""
function btable(n::Integer, center=O; len=ndigits(n; base=2), width=20, title=false)
   sethue("black")
   if title
      t = Table(2, len, width, width, center)
   else
      t = Table(1, len, width, width, center)
   end

   bstr = string(n, base=2, pad=len)
   fontsize(width ÷ 2)

   # start plotting
   col = 1
   if title
      for k in 1:size(t, 2)
         text(string(len-k+1), t[1, k], halign=:center, valign=:middle)
      end
      col += 1
   end

   for j in 1:size(t, 2)
      box(t[col, j], width, width, :stroke)
      text(string(bstr[j]), t[col, j], halign=:center, valign=:middle)
   end

   return t
end

function btable(f::Function, n::Integer, center=O;
   len=ndigits(n; base=2), width=20, title=false)
   gsave()
   # before
   before = btable(n, O - (0, 50); width=50, title=true)
   sethue("blue")
   arrow(Point(0, 10), Point(0, 50); linewidth=2)
   # after
   after = btable(f(n), O+(0, 80); width=50, title=false)
   grestore()
   return (before, after)
end

function btable(f::Function, n::BitStr, center=O; width=50, title=false)
   gsave()
   # before
   before = btable(n, O - (0, 50); width=50, title=true)
   sethue("blue")
   arrow(Point(0, 10), Point(0, 50); linewidth=2)
   # after
   after = btable(f(n), O+(0, 80); width=50, title=false)
   grestore()
   return (before, after)
end

btable(n::BitStr, center=O; width=50, title=false) =
   btable(n.val, center; len=length(n), width=width, title=title)

# TODO: push this to upstream: Luxor
function swaparrow(left::Point, right::Point, action=:fill;
    linewidth=1.0, headlength=10, headangle=pi/8, target_angle=pi/3,
    fraction=1/5, debug=false)
    gsave()
    v = right - left
    delta = fraction * v.x
    C1 = Point(left.x + delta, left.y - delta * tan(target_angle))
    C2 = Point(right.x - delta, right.y - delta * tan(target_angle))
    arrow(left, C1, C2, right, action;
        linewidth=linewidth, rightarrow=true, leftarrow=true)
    if debug
        line(left, C1, :stroke)
        line(right, C2, :stroke)
    end
    grestore()
end

@png begin
   t, _ = btable(breflect, bit"1011"; width=50, title=true)
   lineheight = 20
   sethue("red")
   swaparrow(t[1, 1] - (0, lineheight), t[1, 4] - (0, lineheight), :stroke; linewidth=2, fraction=0.25, target_angle=deg2rad(60))
   swaparrow(t[1, 2] - (0, lineheight), t[1, 3] - (0, lineheight), :stroke; linewidth=2, target_angle=deg2rad(80))
end 400 400 "docs/src/assets/breflect"


text:
struct GType

end


text:
module Fakerjl
include("Person.jl")
include("Communication.jl")
include("utils.jl")

mutable struct Generator
    resrc::Dict{String, String}
    lang::String
end

curr_path = @__DIR__

function init(src::String=curr_path, lang::String="en_US")
    resource = Dict{String, String}()
    all_files = readdir("$curr_path/data/$lang")
    [resource[file]="$curr_path/data/$lang/$file" for file = all_files]
    global gen = Generator(resource, lang)
    return gen
end

function setlanguage!(newlang::String)
    gen.lang = newlang
end

function setsource!(newsource::String)
    empty!(gen.resrc)
    resource = Dict{String, String}()
    all_files = readdir(newsource)
    [resource[file]="$newsource/$file" for file = all_files]
    gen.resrc = resource
end

function readsource(file::String)
    try
        numlines = readlines(gen.resrc[file])
        return length(numlines) > 0 ? numlines : error("File is empty!")
    catch SystemError
        return error("File does not exist!")
    end
end

function genericgenerator(pre::String)
    x = readsource(pre)
    randomline = getrandom(x)
    regex = r"{{.*?}}"
    matches = eachmatch(regex, randomline)
    for m in matches
        tokentype = typeoftoken(m.match[3:end-2])
        if tokentype == "function"
            expr = Meta.parse(m.match[4:end-3])
            evl = eval(expr)
            randomline = replace(randomline, regex => evl, count=1)
        else
            file = String(m.match[3:end-2])
            output = genericgenerator(file)
            randomline = replace(randomline, regex => output, count=1)
        end
    end
    return randomline
end

end # end of module


text:
#######################################################################
###    The s-orbital Toy model                                      ###
#######################################################################

export ToyTBModel
export ToyHamiltonian

using JuLIP.Potentials: PairPotential, Morse, SWCutoff, ZeroSitePotential

"""
`ToyHamiltonian`: constructs a simple s-orbital tight binding hamiltonian.
The hopping function
is given by any pair potential, but the default is the morse potential

  h(r) = e0 ⋅ e^{- α (r/r0 - 1)} ⋅ fcut(r; rcut)

where fcut is a Stillinger-Weber type cutoff.

See also `ToyTBModel`.

### Constructors:

* `ToyHamiltonian(V)`  where `V <: PairPotential`
* `ToyHamiltonian(; kwargs...)` where admissible kwargs with defaults are
   `alpha = 2.0, r0 = 1.0, rcut = 2.5, e0 = 10.0`.
"""
mutable struct ToyHamiltonian{VT <: PairPotential} <: SKHamiltonian{ORTHOGONAL, 1}
   V::VT
end

ToyHamiltonian(; alpha = 2.0, r0 = 1.0, rcut = 2.5, e0 = 10.0) =
   ToyHamiltonian(Morse(e0=e0, A=alpha, r0=r0) * SWCutoff(rcut, 1.0))

cutoff(H::ToyHamiltonian) = cutoff(H.V)

hop(H::ToyHamiltonian, r::Real, ::Integer) = H.V(r)

onsite!(::ToyHamiltonian, _1, _2, H_nn) = setindex!(H_nn, 0.0, 1)

onsite_grad!(::ToyHamiltonian, _1, _2, dH_nn) = fill!(dH_nn, 0.0)

"""
`ToyTBModel`: constructs a simple 1-orbital SK-type tight binding model,
with Hamiltonian given by `ToyHamiltonian`. It doesn't model anything but can
be used for quick tests, e.g. in conjunction
with the MaterialsScienceTools.TriangularLattice module.

### Keyword Parameters

* alpha = 2.0, r0 = 1.0, rcut = 2.7  : Morse potential parameters
* beta = 1.0 : Fermi temperature
* fixed_eF = true : if true, then the chemical potential is fixed (default at 0.0)
* eF = 0.0 : chemical potential (if fixed)
* hfd = 1e-6 : finite difference step for computing hessians
"""
function ToyTBModel(; beta = 1.0, fixed_eF = true, eF = 0.0,
             hfd = 1e-6, bzquad = GammaPoint(), kwargs...)
   H = ToyHamiltonian(;kwargs...)
   TBModel( H, ZeroSitePotential(),
            FermiDiracSmearing(beta, eF, 0, fixed_eF), bzquad, hfd )
end


text:
# --------------module PhyExam-----------
# huangrzh@icloud.com
#---------------------------------------------------

module PhyExam

# package code goes here
export ED_1D_hessenberg

type Ops
  eye::Array{Float64,2}
  sp::Array{Float64,2}
  sz::Array{Float64,2}
  heff::Array{Float64,2}
end

type SuperMatrix
  diml_::Int64
  dimr_::Int64
  dim_::Int64
  opl_::Ops
  opr_::Ops
end

import Base.issymmetric
function issymmetric(x::SuperMatrix)
  return true
end

import Base.eltype
function eltype(x::SuperMatrix)
  return eltype(rand(2));
end

import Base.A_mul_B!
function A_mul_B!(v2_, x::SuperMatrix, v1_)
  v1 = reshape(v1_, x.diml_, x.dimr_);
	v2 = x.opl_.heff * v1;
	v2 = v2 + v1 * x.opr_.heff;

	v_1 = 0.5 * x.opl_.sp * v1 * x.opr_.sp;
  v_2 = 0.5 * x.opl_.sp' * v1 * x.opr_.sp';
	v2 = v2 + v_1 + v_2 + x.opl_.sz * v1 * x.opr_.sz;
  v2 = reshape(v2, x.dim_);
  v2_[:] = v2[:];
end

import Base.size
function size(x::SuperMatrix)
  return x.dim_,x.dim_
end

function ED_1D_hessenberg(latl::Int64, jz::Float64)
sp_ = [0. 1.;0. 0.];
sz_ = [0.5*sqrt(jz) 0.;0. -0.5*sqrt(jz)];
heff_ = zeros(2,2);
eye_ = eye(2,2);
op_l = Ops(eye_, sp_, sz_, heff_);
half_l = convert(Int64, latl/2);
dim = 2;
for i_l in 1:1:half_l-1
	dim = 2^(i_l + 1);
	h_i = 0.5*kron(op_l.sp, sp_');
	op_l.heff = kron(op_l.heff,eye(2,2)) + h_i + h_i' + kron(op_l.sz,sz_);
	op_l.sp = kron(op_l.eye, sp_);
	op_l.sz = kron(op_l.eye, sz_);
	op_l.eye = eye(dim, dim);
end
op_r = op_l;

# ED
ED_object = SuperMatrix(dim, dim, dim*dim, op_l, op_r);
E_gs = eigs(ED_object, nev=1, which=:SR);
return E_gs[1];
end

end # module


text:
using .NavAbilityCaesarExtensions
using NavAbilitySDK
using DistributedFactorGraphs
using RoME
using JSON2
using Base64

import DistributedFactorGraphs: addFactor!

function gqlFactorToDfg(dfg, gqlFactor)
  base64JsonData = get(gqlFactor,"data",nothing)
  if base64JsonData !== nothing
    gqlFactor["data"] = JSON2.read(String(base64decode(base64JsonData)))
  end
  return unpackFactor(dfg, gqlFactor)
end

function addFactor!(dfg::NavAbilityDFG, factor::AbstractDFGFactor; maxWaitSeconds::Int = 150)
  context = Client(dfg.userId,dfg.robotId,dfg.sessionId)
  sdkFactor = packFactor(dfg, factor)
  requestId = addPackedFactor(dfg.navabilityClient, context, sdkFactor)
  label = string(getLabel(factor))
  for i in 1:maxWaitSeconds
    @info i
    if exists(dfg, label)
      break
    end
    sleep(1)
  end
  return requestId
end

function getFactor(dfg::NavAbilityDFG, label::Union{Symbol, String})
  context = Client(dfg.userId,dfg.robotId,dfg.sessionId)
  gqlFactor = NavAbilitySDK.getFactor(dfg.navabilityClient, context, string(label)) |> fetch
  return gqlFactorToDfg(dfg, gqlFactor)
end

function getFactors(dfg::NavAbilityDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[], solvable::Int=0)
  context = Client(dfg.userId,dfg.robotId,dfg.sessionId)
  gqlFactors = NavAbilitySDK.getFactors(dfg.navabilityClient, context; detail=FULL) |> fetch
  dfgFactors = map(f -> gqlFactorToDfg(dfg, f),gqlFactors)
  # TODO: Implement regexFilter
  # TODO: Implement tags
  # TODO: Implement solver filter
  return dfgFactors
end


text:
using CollegeEntry, ModelObjectsLH, ModelParams
using Test, TestSetExtensions

include("test_helpers.jl")

@testset "All" begin
    include("helpers_test.jl")
    include("admissions_test.jl");
    include("admission_prob_test.jl");
    include("student_rankings_test.jl")
    include("entry_test.jl");
    include("entry_decisions_test.jl")
    include("entry_results_test.jl")
end

# ----------


text:
type MyType
   x::Array{Float64,1}
   str::ASCIIString
   MyType() = new()
   MyType(a::Array{Float64,1}, b::ASCIIString) = new(a,b)
end

mkexp(str,x...) = :($(parse(str))($(x...)))
# mkexp(ty,str,x...) = :(ty$(parse(str))($(x...)))
strT = "MyType"
@show ex1 = mkexp(strT);
@show typeof(ex1);
mt1 = eval(ex1)

ex2 = mkexp(strT);
mt1 = eval(ex2)


text:
mutable struct SuffixArray
    data::Vector{Int}
    lcparray::Vector{Int}
end

Base.getindex(sa::SuffixArray, i::Int) = sa.data[i]

function SuffixArray(text::Vector{Int})
    sa = sais(text)
    lcparray = StringMatch.lcparray(sa, text)
    SuffixArray(sa, lcparray)
end

"""
    lcparray

Kasai's algorithm for linear-time construction of LCP array from Suffix Array
"""
function lcparray(sa::Vector{Int}, text::Vector{T}) where T<:Integer
    n = length(sa)
    lcps = similar(sa)
    rank = similar(sa)
    for i = 1:n
        rank[sa[i]] = i
    end

    lcp = 0
    for i = 1:n
        rank[i] > 1 || continue
        j = sa[rank[i]-1]
        while i+lcp <= n && j+lcp <= n && text[i+lcp] == text[j+lcp]
            lcp += 1
        end
        lcps[rank[i]] = lcp
        lcp > 0 && (lcp -= 1)
    end
    lcps[1] = 0
    lcps
end


text:
abstract Any_Area_is_an_Area

type AnArea_Regular <: Any_Area_is_an_Area
  data_based_on :: ASCIIString
  lon  :: Array{Float64,1}
  lat  :: Array{Float64,1}
  plon :: Array{Float64,2}
  plat :: Array{Float64,2}
  plon_selected_on_by :: Array{Bool,2}
  data_selected_on_plon_by :: Array{Bool,2}
end

# Select time series
function Selected_within{N}(p_area :: AnArea_Regular,
                            sst :: Array{Float64,N}; map_to :: Symbol = :default)
  plon_s = p_area.plon_selected_on_by
  data_s = p_area.data_selected_on_plon_by
  sst1 = reshape(sst ,length(plon_s), reduce(*,size(sst,[3:ndims(sst);]...)))

  if (map_to == :default )
    sst1 = sst1[find(plon_s[:]),:];
    sst1[find(!data_s[:]),:] = NaN
    sst1 = reshape(sst1,size(data_s)...,size(sst,(3:ndims(sst))...)...)
  elseif(map_to == :whole )
    Filter1= find(plon_s[:])
    Filter2 = Filter1[find(data_s[:])]
    pFilter = fill(false,plon_s); pFilter[Filter2] = true;
    sst1[find(!pFilter[:]),:] = NaN
    sst1 = reshape(sst1, size(sst))
  else
    error("No such 'map_to' stratagy yet")
  end
  return sst1
end

function Selected_within(p_area :: AnArea_Regular,
                         sst :: Array{Float64,2}; map_to :: Symbol = :default)
  plon_selected_on_by = p_area.plon_selected_on_by
  data_selected_on_plon_by = p_area.data_selected_on_plon_by
  if (map_to == :default )
    sst1 = reshape(sst[plon_selected_on_by],size(data_selected_on_plon_by))
    sst1[!data_selected_on_plon_by] = NaN
  elseif(map_to == :whole )
    sst1 = copy(sst)
    dreal = find(plon_selected_on_by)
    dreal = dreal[find(data_selected_on_plon_by)]
    mask = fill(true,size(plon_selected_on_by) )
    mask[dreal] = false
    sst1[mask] = NaN
  else
    error("No such 'map_to' stratagy yet")
  end
  return sst1
end

export Any_Area_is_an_Area, AnArea_Regular
export Selected_within
include("deal_with_models/Regular_Area_Select.jl")
include("special_areas.jl")


text:
get_input() = open("./inputs/day15.txt") do f
	readlines(f)[1] |> l -> split(l, ",") .|> string .|> String |> l -> parse.(Int, l)
end

function run_game(upto)
	nums = get_input()

	prev = pop!(nums)

	memo = Dict((num, i) for (i, num) in enumerate(nums))

	turn = length(nums) + 1

	while turn != upto
		if haskey(memo, prev)
			memo[prev], prev = turn, turn - memo[prev]
		else
			memo[prev] = turn
			prev = 0
		end

		turn += 1
	end

	return prev
end

part1() = run_game(2020)
part2() = run_game(30000000)

println("day 15:\n\tpart 1: ", part1(), "\n\tpart 2: ", part2())


text:
#=
P21 (*) Insert an element at a given position into a list.
Example:
* (insert-at 'alfa '(a b c d) 2)
(A ALFA B C D)
=#

function insert_at(item, coll, index)
    vcat(coll[1:index - 1], [item], coll[index:end])
end

@assert insert_at(:alfa, [:a, :b, :c, :d], 2) == [:a, :alfa, :b, :c, :d]
@assert insert_at(0, [1, 2, 3], 1) == [0, 1, 2, 3]

println("Tests passed: JL-21.jl")


text:
function kernel(x::Float64)
    y = rand(:y, Normal(x, 1.0))
    return y
end
test_markov = () -> markov(:k, kernel, 5, 1.0)
test_plate = () -> plate(:k, kernel, [1.0, 2.0, 3.0, 4.0, 5.0])

@testset "Simulate" begin
    ret, cl = simulate(test_plate)
    for i in 1:5
        @test haskey(cl, (:k, i, :y))
    end
    ret, cl = simulate(test_markov)
    for i in 1:5
        @test haskey(cl, (:k, i, :y))
    end
end

@testset "Constrained generate" begin
    sel = selection([(:k, 3, :y) => 5.0])
    ret, cl, _ = generate(sel, test_plate)
    for i in 1:5
        @test haskey(cl, (:k, i, :y))
    end
    @test get_ret(cl[:k, 3, :y]) == 5.0
    ret, cl, _ = generate(sel, test_markov)
    for i in 1:5
        @test haskey(cl, (:k, i, :y))
    end
    @test get_ret(cl[:k, 3, :y]) == 5.0
end

@testset "Update" begin
    @testset "Vectorized plate" begin
        ret, cl = simulate(test_plate)
        original_score = get_score(cl)
        stored_at_y = get_ret(cl[:k, 3, :y])
        sel = selection([(:k, 3, :y) => 5.0])
        ret, cl, w, rd, d = update(sel, cl)
        # Discard should be original :x
        #@test d[:k => 3 => :y] == stored_at_y
        # New should be equal to constraint.
        #@test cl[:k => 3 => :y] == 5.0
        @test get_score(cl) - w ≈ original_score
    end

    @testset "Vectorized markov" begin
        ret, cl = simulate(test_markov)
        stored_at_y = get_ret(cl[:k, 3, :y])
        sel = selection([(:k, 3, :y) => 5.0])
        ret, cl, w, retdiff, d = update(sel, cl)
        # Discard should be original :x
        #@test d[:k => 3 => :y] == stored_at_y
        # New should be equal to constraint.
        @test get_ret(cl[:k, 3, :y]) == 5.0
    end
end

@testset "Regenerate" begin
    @testset "Vectorized plate" begin
        ret, cl = simulate(test_plate)
        original_score = get_score(cl)
        stored = get_ret(cl[:k, 3, :y])
        sel = selection([(:k, 3, :y)])
        ret, cl, w, retdiff, d = regenerate(sel, cl)
        # Discard should be original :x
        #@test d.query[:x] == stored_at_x
        # New should not be equal to original.
        @test get_ret(cl[:k, 3, :y]) != stored
        @test get_score(cl) - w ≈ original_score
    end

    @testset "Vectorized markov" begin
        ret, cl = simulate(test_markov)
        original_score = get_score(cl)
        stored = get_ret(cl[:k, 3, :y])
        sel = selection([(:k, 3, :y)])
        ret, cl, w, retdiff, d = regenerate(sel, cl)
        # Discard should be original :x
        #@test d.query[:x] == stored_at_x
        # New should not be equal to original.
        @test get_ret(cl[:k, 3, :y]) != stored
        @test get_score(cl) - w ≈ original_score
    end
end

@testset "Propose" begin
    ret, cl, w = propose(test_markov)
    @test get_score(cl) != 0.0
    @test get_score(cl) == w
    ret, cl, w = propose(test_plate)
    @test get_score(cl) != 0.0
    @test get_score(cl) == w
end

@testset "Score" begin
    sel = selection([(:k, 1, :y) => 5.0,
                     (:k, 2, :y) => 5.0,
                     (:k, 3, :y) => 5.0,
                     (:k, 4, :y) => 5.0,
                     (:k, 5, :y) => 5.0])
    ret, sc = score(sel, test_markov)
    @test sc != 0.0
    ret, sc = score(sel, test_plate)
    @test sc != 0.0
end


text:
"""
    img = nodegrid(T, nodes)
    img = nodegrid(nodes)
    img = nodegrid([T], ϕ)

Returns an image `img` which has value 1 at points that are "on the
nodes."  If `nodes`/`ϕ` is two-dimensional, the image will
consist of grid lines; for three-dimensional inputs, the image will
have grid planes.

The meaning of "on the nodes" is that `x[d] == nodes[d][i]` for some
dimension `d` and node index `i`. An exception is made for the edge
values, which are brought inward by one pixel to prevent complete loss
under subpixel deformations.

Optionally specify the element type `T` of `img`.

See also [`warpgrid`](@ref).
"""
function nodegrid(::Type{T}, nodes::NTuple{N,AbstractRange}) where {T,N}
    @assert all(r->first(r)==1, nodes)
    inds = map(r->Base.OneTo(ceil(Int, last(r))), nodes)
    img = Array{T}(undef, map(length, inds))
    fill!(img, zero(T))
    for idim = 1:N
        indexes = Any[inds...]
        indexes[idim] = map(x->clamp(round(Int, x), first(inds[idim])+1, last(inds[idim])-1), nodes[idim])
        img[indexes...] .= one(T)
    end
    img
end
nodegrid(::Type{T}, ϕ::GridDeformation) where {T} = nodegrid(T, ϕ.nodes)
nodegrid(arg) = nodegrid(Bool, arg)

"""
    img = warpgrid(ϕ; [scale=1, showidentity=false])

Returns an image `img` that permits visualization of the deformation
`ϕ`.  The output is a warped rectangular grid with nodes centered on
the control points as specified by the nodes of `ϕ`.  If `ϕ` is a
two-dimensional deformation, the image will consist of grid lines; for
a three-dimensional deformation, the image will have grid planes.

`scale` multiplies `ϕ.u`, effectively making the deformation stronger
(for `scale > 1`).  This can be useful if you are trying to visualize
subtle changes. If `showidentity` is `true`, an RGB image is returned
that has the warped grid in magenta and the original grid in green.

See also [`nodegrid`](@ref).
"""
function warpgrid(ϕ; scale=1, showidentity::Bool=false)
    img = nodegrid(ϕ)
    if scale != 1
        ϕ = GridDeformation(scale*ϕ.u, ϕ.nodes)
    end
    wimg = warp(img, ϕ)
    if showidentity
        n = ndims(img)+1
        return reshape(reinterpret(RGB{Float32}, permutedims(cat(wimg, img, wimg, dims=n), (n,1:ndims(img)...))),(size(img)...,))
    end
    wimg
end


text:
nt = 1000
seed = 99
m = ACQR_simula(nt,seed)
y = m["y"]

############################### starting point
wi = collect(range(0.1, length=2, stop=3))
Ai,Ci,Qi,Ri,Si,m1i,P1i = ACQR_em_starting(y,wi,1,true)

#######################
max_iter=10
tol=1e-6
txo=false
print("Testing ACQR_em: ")
A,C,Q,R,m1,P1,loglik,aic = ACQR_em(y,Ai,Ci,Qi,Ri,m1i,P1i,max_iter,tol,txo)
if sum(A.^2) > 1e-6
    println("OK")
end


text:
module MPIClusterManagers

export MPIManager, launch, manage, kill, procs, connect, mpiprocs, @mpi_do, TransportMode, MPI_ON_WORKERS, TCP_TRANSPORT_ALL, MPI_TRANSPORT_ALL

using Distributed, Serialization
import MPI
const mpiexec = isdefined(MPI, :mpiexec) ? MPI.mpiexec : "mpiexec"

include("mpimanager.jl")

end # module


text:
using ControlSystemIdentification, ControlSystems, Statistics, Random
Random.seed!(0)
wtest = exp10.(LinRange(-3, log10(pi), 30))
freqresptest(G, model) = quantile(
    vec(log10.(abs2.(freqresp(model, wtest))) - log10.(abs2.(freqresp(G, wtest)))),
    0.9,
)

freqresptest(G, model, tol) = freqresptest(G, model) < tol

N = 3000
h = 0.05
t = range(0, length = N, step = h)

ω1, ω2 = 10, 20
f1, f2 = ω1 / 2π, ω2 / 2π
ωvec = exp10.(LinRange(-0, log10(pi / h), 300))

G1 = c2d(tf(ω1, [1, 2 * 0.02 * ω1, ω1^2]), h)
G2 = c2d(tf(ω2, [1, 2 * 0.02 * ω2, ω2^2]), h)
# bodeplot([G1, G2], ωvec)

e = randn(N)
y1, y2 = lsim(G1, e, t)[1], lsim(G2, e, t)[1]
d = iddata(y1 + y2, e, h)

# Gh = arx(d,2,2)
# bodeplot([G1, G2, Gh], ωvec)

##

# P = mt_pgram(vec(d.y), fs=1/h)

H1 = Bandstop(0.9 * f1, 1.1 * f1, fs = d.fs)
H2 = Bandstop(0.9 * f2, 1.1 * f2, fs = d.fs)

Gf1 = arx(d, 2, 2, estimator = weighted_estimator(H1))
Gf2 = arx(d, 2, 2, estimator = weighted_estimator(H2))

@test norm(pole(G1) - pole(Gf1)) / norm(pole(G1)) < 0.05
@test norm(pole(G2) - pole(Gf2)) / norm(pole(G1)) < 0.05

@test freqresptest(G1, Gf1) < 0.22
@test freqresptest(G2, Gf2) < 0.22

isinteractive() && bodeplot(
    [G1, G2, Gf1, Gf2],
    ωvec,
    plotphase = false,
    lab = ["G1" "G2" "Focus f1" "Focus f2"],
)
# plot!(2pi .* P.freq, sqrt.(P.power), l=(:black, 2, :dash), xlims=extrema(ωvec))

Random.seed!(0)
N = 2000
e = randn(N)
t = range(0, length = N, step = h)
y1, y2 = lsim(G1, e, t)[1], lsim(G2, e, t)[1]
d = iddata(y1 + y2, e, h)

# d.y .+= 0.1 .* randn.()

H1 = Bandstop(0.5 * f1, 1.2 * f1, fs = 1 / h)
H2 = Bandstop(0.9 * f2, 1.5 * f2, fs = 1 / h)

Gf1 = n4sid(d, 2, Wf = H1)
Gf2 = n4sid(d, 2, Wf = H2)
Gf1 = Gf1.sys
Gf2 = Gf2.sys

@test norm(pole(G1) - pole(Gf1)) / norm(pole(G1)) < 0.05
@test norm(pole(G2) - pole(Gf2)) / norm(pole(G1)) < 0.05

@test freqresptest(G1, Gf1) < 1
@test freqresptest(G2, Gf2) < 1

# isinteractive() && bodeplot([G1, G2, Gf1, Gf2], ωvec, plotphase=false, lab=["G1" "G2" "Focus f1" "Focus f2"])