text:
using Statistics
using Plots
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using JLD2
using Random
using Distributions
using FITSIO
using LineSearches
using Flux
using StatsBase
using SparseArrays

push!(LOAD_PATH, pwd()*"../main")
using DHC_2DUtils
push!(LOAD_PATH, pwd()*"../scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs

using IJulia
IJulia.notebook(dir="/home/nayantara/Desktop/GalacticDustProject/WST/DHC/scratch_NM/NewWrapper/4-28/")
#=
#Reg, RegCoeff
logbool = false
apdbool=true
isobool = false
Nx=64
filter_hash = fink_filter_hash(1, 8, nx=Nx, t=1, wd=1, Omega=true)
sfdall = readsfd_fromsrc("data/dust10000.fits", Nx, logbool=logbool)
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool)
dbnocffs = get_dbn_coeffs(sfdall, filter_hash, dhc_args)
fmean = mean(dbnocffs, dims=1)
scov= (dbnocffs .- fmean)' * (dbnocffs .- fmean) ./(size(dbnocffs)[1] -1)
cond(scov)
scovinv = invert_covmat(scov, 1e-5)
cholupp = cholesky(scovinv).U
dbnmean =
whitened_ms = (cholupp * (dbnocffs .-  mean(dbnocffs, dims=1))')'
selidx = [3, 4]
p = scatter(whitened[:, selidx[1]], whitened[:, selidx[2]], title="Real im, Not log coeffs")
xlabel!("J=1, L=0")
ylabel!("J=2, L=0")
filter_hash["J_L"]
selidx = [3, 7]
p = scatter(whitened[:, selidx[1]], whitened[:, selidx[2]], title="Real im, Not log coeffs")
xlabel!("J=1, L=0")
ylabel!("J=1, L=1")

#Reg, LogCoeff
logdbn = log.(dbnocffs)
flmean = mean(logdbn, dims=1)
scov= (logdbn .- flmean)' * (logdbn .- flmean) ./(size(logdbn)[1] -1)
cond(scov)
scovinv = invert_covmat(scov, 1e-10)
cholupp = cholesky(scovinv).U
whitened = (cholupp * logdbn')'
selidx = [3, 4]
p = scatter(whitened[:, selidx[1]], whitened[:, selidx[2]], title="Real im, Log coeffs")
xlabel!("J=1, L=0")
ylabel!("J=2, L=0")
filter_hash["J_L"]
selidx = [3, 7]
p = scatter(whitened[:, selidx[1]], whitened[:, selidx[2]], title="Real im, Log coeffs")
xlabel!("J=1, L=0")
ylabel!("J=1, L=1")

#Log
logbool = true
apdbool=true
isobool = false
Nx=64
filter_hash = fink_filter_hash(1, 8, nx=Nx, t=1, wd=1, Omega=true)
sfdall = readsfd_fromsrc("data/dust10000.fits", Nx, logbool=logbool)
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool)
dbnocffs = get_dbn_coeffs(sfdall, filter_hash, dhc_args)
fmean = mean(dbnocffs, dims=1)
scov= (dbnocffs .- fmean)' * (dbnocffs .- fmean) ./(size(dbnocffs)[1] -1)
cond(scov)
scovinv = invert_covmat(scov, 1e-5)
=#

## Reg, RegCoeff
logbool = false
apdbool=true
isobool = true
Nx=64
filter_hash = fink_filter_hash(1, 8, nx=Nx, t=1, wd=1, Omega=true)
sfdall = readsfd_fromsrc("scratch_NM/data/dust10000.fits", Nx, logbool=logbool)
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool)
dbnocffs = get_dbn_coeffs(sfdall, filter_hash, dhc_args)

#Subsection
fmean = mean(dbnocffs, dims=1)
scov= (dbnocffs .- fmean)' * (dbnocffs .- fmean) ./(size(dbnocffs)[1] -1)
cond(scov)
scovinv = invert_covmat(scov, 1e-5)
cholupp = cholesky(scovinv).U
whitened_ms = (cholupp * (dbnocffs .-  mean(dbnocffs, dims=1))')'
selidx = [3, 4]

p = plot(scatter(whitened_ms[:, selidx[1]], whitened_ms[:, selidx[2]], title="Real im, Not log coeffs", xlabel="J=1", ylabel = "J=2"), size=(600, 600))

eps = randn(156, 10000)
covsq = cholesky(scov+Diagonal(ones(156).*1e-5)).U
qsamps = fmean .+ (covsq*eps)'

p = scatter(qsamps[1:1000, selidx[1]], qsamps[1:1000, selidx[2]], label="Samples from Gaussian", color="red", xlabel="J=1", ylabel = "J=2", xlim=(-100, 100), ylim=(-100, 100))
scatter!(dbnocffs[1:1000, selidx[1]], dbnocffs[1:1000, selidx[2]], label="True Samples", color="black", xlabel="J=1", ylabel = "J=2", xlim=(-100, 100), ylim=(-100, 100))

##Using Phi cut to remove bad images
function distribution_percentiles(dbn, idx, perc)
    idxperc = zeros(size(idx))
    idxperc .= (x-> percentile(dbn[:, x], perc)).(idx)
    return idxperc
end

phipval = distribution_percentiles(dbnocffs, [7], 99)
cmask = (dbnocffs[:, 7] .> phipval)
badimg = findall(!iszero, cmask)
bidx = rand(1:length(badimg))
heatmap(sfdall[:, :, badimg[bidx]], title="Phi Value=" * string(dbnocffs[badimg[bidx], 7]))

##Removing bad images based on phi cut
goodimg = findall(iszero, cmask)
subdbn = dbnocffs[goodimg, :]
fmean = mean(subdbn, dims=1)
scov= (subdbn .- fmean)' * (subdbn .- fmean) ./(size(subdbn)[1] -1)
cond(scov)
scovinv = invert_covmat(scov, 1e-5)
cholupp = cholesky(scovinv).U
whitened_ms = (cholupp * (dbnocffs .-  mean(dbnocffs, dims=1))')'
selidx = [3, 4]

p = plot(scatter(whitened_ms[:, selidx[1]], whitened_ms[:, selidx[2]], title="Real im, Not log coeffs", xlabel="J=1", ylabel = "J=2"), size=(600, 600))

eps = randn(156, 10000)
covsq = cholesky(scov+Diagonal(ones(156).*1e-5)).U
qsamps = fmean .+ (covsq*eps)'

p = scatter(qsamps[1:1000, selidx[1]], qsamps[1:1000, selidx[2]], label="Samples from Gaussian", color="red", xlabel="J=1", ylabel = "J=2", xlim=(-100, 100), ylim=(-100, 100))
scatter!(subdbn[1:1000, selidx[1]], subdbn[1:1000, selidx[2]], label="True Samples", color="black", xlabel="J=1", ylabel = "J=2", xlim=(-100, 100), ylim=(-100, 100))

##wd=2
logbool = false
apdbool=true
isobool = false
Nx=64
filter_hash = fink_filter_hash(1, 8, nx=Nx, t=1, wd=2, Omega=true)
sfdall = readsfd_fromsrc("scratch_NM/data/dust10000.fits", Nx, logbool=logbool)
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool)
dbnocffs = get_dbn_coeffs(sfdall, filter_hash, dhc_args)

#Subsection
fmean = mean(dbnocffs, dims=1)
scov= (dbnocffs .- fmean)' * (dbnocffs .- fmean) ./(size(dbnocffs)[1] -1)
cond(scov)
scovinv = invert_covmat(scov, 1e-5)
cholupp = cholesky(scovinv).U
whitened_ms = (cholupp * (dbnocffs .-  mean(dbnocffs, dims=1))')'
selidx = [3, 4]

p = plot(scatter(whitened_ms[:, selidx[1]], whitened_ms[:, selidx[2]], title="Real im, Not log coeffs", xlabel="J=1", ylabel = "J=2"), size=(600, 600))

eps = randn(1192, 10000)
covsq = cholesky(scov+Diagonal(ones(1192).*1e-5)).U
qsamps = fmean .+ (covsq*eps)'

p = scatter(qsamps[1:1000, selidx[1]], qsamps[1:1000, selidx[2]], label="Samples from Gaussian", color="red", xlabel="J=1", ylabel = "J=2", xlim=(-100, 100), ylim=(-100, 100))
scatter!(dbnocffs[1:1000, selidx[1]], dbnocffs[1:1000, selidx[2]], label="True Samples", color="black", xlabel="J=1", ylabel = "J=2", xlim=(-100, 100), ylim=(-100, 100))

dbnocffs= log.(dbnocffs)

fmean = mean(dbnocffs, dims=1)
scov= (dbnocffs .- fmean)' * (dbnocffs .- fmean) ./(size(dbnocffs)[1] -1)
cond(scov)
scovinv = invert_covmat(scov, 1e-5)
cholupp = cholesky(scovinv).U
whitened_ms = (cholupp * (dbnocffs .-  mean(dbnocffs, dims=1))')'
selidx = [3, 4]

p = plot(scatter(whitened_ms[:, selidx[1]], whitened_ms[:, selidx[2]], title="Real im, Log coeffs", xlabel="J=1", ylabel = "J=2"), size=(600, 600))

eps = randn(1192, 10000)
covsq = cholesky(scov+Diagonal(ones(1192).*1e-5)).U
qsamps = fmean .+ (covsq*eps)'

p = scatter(qsamps[1:1000, selidx[1]], qsamps[1:1000, selidx[2]], label="Samples from Gaussian", color="red", xlabel="J=1", ylabel = "J=2", xlim=(-100, 100), ylim=(-100, 100))
scatter!(dbnocffs[1:1000, selidx[1]], dbnocffs[1:1000, selidx[2]], label="True Samples", color="black", xlabel="J=1", ylabel = "J=2", xlim=(-100, 100), ylim=(-100, 100))

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module SerialPorts

export SerialPort, SerialException, setDTR, list_serialports,
       check_serial_access

using Compat, PyCall
VERSION < v"0.4-" && using Docile

const PySerial = PyCall.PyNULL()
const PySerialListPorts = PyCall.PyNULL()

type SerialException <: Base.Exception end

immutable SerialPort <: IO
    port::ASCIIString
    baudrate::Int
    bytesize::Int
    parity::ASCIIString
    stopbits::Int
    timeout
    xonxoff::Bool
    rtscts::Bool
    dsrdtr::Bool
    python_ptr::PyObject
end

function __init__()
    copy!(PySerial, pyimport_conda("serial", "pyserial"))
    copy!(PySerialListPorts, pyimport("serial.tools.list_ports"))
end

function serialport(port, baudrate)
    py_ptr = PySerial[:Serial](port, baudrate)
    SerialPort(port,
               baudrate,
               py_ptr[:bytesize],
               py_ptr[:parity],
               py_ptr[:stopbits],
               py_ptr[:timeout],
               py_ptr[:xonxoff],
               py_ptr[:rtscts],
               py_ptr[:dsrdtr], py_ptr)
end

if VERSION >= v"0.4-"
    function Base.call(::Type{SerialPort}, port, baudrate)
        serialport(port, baudrate)
    end
end

function Base.open(serialport::SerialPort)
    serialport.python_ptr[:open]()
    return serialport
end

function Base.close(serialport::SerialPort)
    serialport.python_ptr[:close]()
    return serialport
end

function Base.isreadable(ser::SerialPort)
    ser.python_ptr[:isreadable]()
end

function Base.iswritable(ser::SerialPort)
    ser.python_ptr[:iswritable]()
end

function Base.write(serialport::SerialPort, data::@compat UInt8)
    serialport.python_ptr[:write](data)
end

function Base.write(serialport::SerialPort, data::ASCIIString)
    serialport.python_ptr[:write](data)
end

function Base.read(ser::SerialPort, bytes::Integer)
    ser.python_ptr[:read](bytes)
end

function Base.nb_available(ser::SerialPort)
    ser.python_ptr[:inWaiting]()
end

function Base.readavailable(ser::SerialPort)
    read(ser, nb_available(ser))
end

function setDTR(ser::SerialPort, val)
    ser.python_ptr[:setDTR](val)
end

function _valid_linux_port(x)
    startswith(x, "ttyS") || startswith(x, "ttyUSB") || startswith(x, "ttyACM")
end

function _valid_darwin_port(x)
    startswith(x, "tty.") || startswith(x, "cu.")
end

@doc """
List available serialports on the system.
""" ->
function list_serialports()
    @unix_only begin
        ports = readdir("/dev/")
        f = @osx ? _valid_darwin_port : _valid_linux_port
        filter!(f, ports)
        return [string("/dev/", port) for port in ports]
    end
    @windows_only begin
        [i[1] for i in collect(PySerialListPorts[:comports]())]
    end
end

@doc """
Check if there are permission issues with accessing serial ports on the current
system.
""" ->
function check_serial_access()
    @linux_only begin
        current_user = ENV["USER"]
        in_dialout() || warn("""User $current_user is not in the 'dialout' group.
                                They can be added with:
                                'usermod -a -G dialout $current_user'""")
    end
end

@linux_only @doc """
On Linux, test if the current user is in the 'dialout' group.
""" ->
function in_dialout()
    "dialout" in split(readall(`groups`))
end

# Submodules

include("Arduino.jl")

end # module

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 9941173e-f2dc-11ea-0f91-4ff7401d7d2c
using Pkg, DrWatson

# ╔═╡ 0b2c7f06-f2dc-11ea-0434-cdc3e04decae
md"## Clip-00-01-03s.jl"

# ╔═╡ 9948d1f4-f2dc-11ea-225d-992c7e2f825f
md"### snippet 0.1"

# ╔═╡ 99496f10-f2dc-11ea-21f9-7fd3bcbdf477
"All models are wrong, but some are useful."

# ╔═╡ 9956c73c-f2dc-11ea-1620-13786d3d33e4
md"### snippet 0.2"

# ╔═╡ 9957701a-f2dc-11ea-18d6-55adda75b7d4
md"##### This is a StepRange, not a vector."

# ╔═╡ 9960bc1a-f2dc-11ea-3a7b-4181c9fa8fe8
x1 = 1:3

# ╔═╡ 996156d4-f2dc-11ea-10e4-4fc246187dd3
md"##### Below still preserves the StepRange."

# ╔═╡ 996b40ea-f2dc-11ea-3814-15033114a487
x2 = x1*10

# ╔═╡ e04af122-f2dc-11ea-2a23-33f8632dbed9
typeof(x2)

# ╔═╡ fbc13db2-f2dc-11ea-1368-9387b64be726
x2[end]

# ╔═╡ 996be34c-f2dc-11ea-2f8f-c389bca0896e
md"##### `Broadcast` log to steprange elements in x2, this returms a vector! Notice the log.(x2) notation."

# ╔═╡ 99780bae-f2dc-11ea-0af3-6157b75c5288
x3 = log.(x2)

# ╔═╡ ed4ea62a-f2dc-11ea-3bed-0ff984265c52
typeof(x3)

# ╔═╡ 9979b86e-f2dc-11ea-3cb9-291cec8627a2
md"##### We can sum the vector x3."

# ╔═╡ 9988b358-f2dc-11ea-1f2d-139af6ac6699
x4 = sum(x3)

# ╔═╡ 998df810-f2dc-11ea-0415-a33278120eac
# Etc.

begin
	x = exp(x4)
	x = x*10
	x = log(x)
	x = sum(x)
	x = exp(x)
end

# ╔═╡ 9991cfc6-f2dc-11ea-0322-f908d09388e7
md"### snippet 0.3"

# ╔═╡ 999f513c-f2dc-11ea-3917-2390217ae1ad
[log(0.01^200) 200 * log(0.01)]

# ╔═╡ 99a096fa-f2dc-11ea-2c15-bd83f1acacf9
md"## End of clip-00-01-03s.jl"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[compat]
DrWatson = "~2.5.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc1"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.DrWatson]]
deps = ["Dates", "FileIO", "LibGit2", "MacroTools", "Pkg", "Random", "Requires", "UnPack"]
git-tree-sha1 = "d6aa02ad618cf31af9bbbf87f87baad632538211"
uuid = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
version = "2.5.0"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "3c041d2ac0a52a12a27af2782b34900d9c3ee68c"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.1"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─0b2c7f06-f2dc-11ea-0434-cdc3e04decae
# ╠═9941173e-f2dc-11ea-0f91-4ff7401d7d2c
# ╟─9948d1f4-f2dc-11ea-225d-992c7e2f825f
# ╟─99496f10-f2dc-11ea-21f9-7fd3bcbdf477
# ╟─9956c73c-f2dc-11ea-1620-13786d3d33e4
# ╟─9957701a-f2dc-11ea-18d6-55adda75b7d4
# ╠═9960bc1a-f2dc-11ea-3a7b-4181c9fa8fe8
# ╟─996156d4-f2dc-11ea-10e4-4fc246187dd3
# ╠═996b40ea-f2dc-11ea-3814-15033114a487
# ╠═e04af122-f2dc-11ea-2a23-33f8632dbed9
# ╠═fbc13db2-f2dc-11ea-1368-9387b64be726
# ╟─996be34c-f2dc-11ea-2f8f-c389bca0896e
# ╠═99780bae-f2dc-11ea-0af3-6157b75c5288
# ╠═ed4ea62a-f2dc-11ea-3bed-0ff984265c52
# ╟─9979b86e-f2dc-11ea-3cb9-291cec8627a2
# ╠═9988b358-f2dc-11ea-1f2d-139af6ac6699
# ╠═998df810-f2dc-11ea-0415-a33278120eac
# ╟─9991cfc6-f2dc-11ea-0322-f908d09388e7
# ╠═999f513c-f2dc-11ea-3917-2390217ae1ad
# ╟─99a096fa-f2dc-11ea-2c15-bd83f1acacf9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
export Ultraspherical

#Ultraspherical Spaces

"""
`Ultraspherical(λ)` is the space spanned by the ultraspherical polynomials
```
    C_0^{(λ)}(x),C_1^{(λ)}(x),C_2^{(λ)}(x),…
```
Note that `λ=1` this reduces to Chebyshev polynomials of the second kind:
`C_k^{(1)}(x) = U_k(x)`.
"""
struct Ultraspherical{T,D<:Domain,R} <: PolynomialSpace{D,R}
    order::T
    domain::D
    Ultraspherical{T,D,R}(m::T,d::D) where {T,D,R} = (@assert m ≠ 0; new(m,d))
    Ultraspherical{T,D,R}(m::Number,d::Domain) where {T,D,R} = (@assert m ≠ 0; new(convert(T,m),convert(D,d)))
    Ultraspherical{T,D,R}(d::Domain) where {T,D,R} = new(one(T),convert(D,d))
    Ultraspherical{T,D,R}(m::Number) where {T,D,R} = (@assert m ≠ 0; new(convert(T,m),D()))
end

Ultraspherical(m::Number,d::Domain) = Ultraspherical{typeof(m),typeof(d),real(prectype(d))}(m,d)
Ultraspherical(m::Number,d) = Ultraspherical(m,Domain(d))
Ultraspherical(m::Number) = Ultraspherical(m,ChebyshevInterval())

order(S::Ultraspherical) = S.order
setdomain(S::Ultraspherical,d::Domain) = Ultraspherical(order(S),d)

convert(::Type{Ultraspherical{T,D,R}}, S::Ultraspherical{T,D,R}) where {T,D,R} = S

canonicalspace(S::Ultraspherical) = Chebyshev(domain(S))
pointscompatible(A::Ultraspherical, B::Chebyshev) = domain(A) == domain(B)
pointscompatible(A::Chebyshev, B::Ultraspherical) = domain(A) == domain(B)

struct UltrasphericalPlan{CT,FT}
    chebplan::CT
    cheb2legplan::FT

    UltrasphericalPlan{CT,FT}(cp,c2lp) where {CT,FT} = new{CT,FT}(cp,c2lp)
end

struct UltrasphericalIPlan{CT,FT}
    chebiplan::CT
    leg2chebplan::FT

    UltrasphericalIPlan{CT,FT}(cp,c2lp) where {CT,FT} = new{CT,FT}(cp,c2lp)
end

function UltrasphericalPlan(λ::Number,vals)
    if λ == 0.5
        cp = plan_transform(Chebyshev(),vals)
        c2lp = FastTransforms.th_cheb2legplan(eltype(vals),length(vals))
        UltrasphericalPlan{typeof(cp),typeof(c2lp)}(cp,c2lp)
    else
        error("Not implemented")
    end
end

function UltrasphericalIPlan(λ::Number,cfs)
    if λ == 0.5
        cp=plan_itransform(Chebyshev(),cfs)
        c2lp=FastTransforms.th_leg2chebplan(eltype(cfs),length(cfs))
        UltrasphericalIPlan{typeof(cp),typeof(c2lp)}(cp,c2lp)
    else
        error("Not implemented")
    end
end

*(UP::UltrasphericalPlan,v::AbstractVector) =
    UP.cheb2legplan*(UP.chebplan*v)
*(UP::UltrasphericalIPlan,v::AbstractVector) =
    UP.chebiplan*(UP.leg2chebplan*v)

plan_transform(sp::Ultraspherical{Int},vals::AbstractVector) = CanonicalTransformPlan(sp,vals)
plan_transform(sp::Ultraspherical,vals::AbstractVector) = UltrasphericalPlan(order(sp),vals)
plan_itransform(sp::Ultraspherical{Int},cfs::AbstractVector) = ICanonicalTransformPlan(sp,cfs)
plan_itransform(sp::Ultraspherical,cfs::AbstractVector) = UltrasphericalIPlan(order(sp),cfs)

## Construction

#domain(S) may be any domain

ones(::Type{T},S::Ultraspherical) where {T<:Number} = Fun(S,fill(one(T),1))
ones(S::Ultraspherical) = Fun(S,fill(1.0,1))

## Fast evaluation

function Base.first(f::Fun{Ultraspherical{Int,D,R}}) where {D,R}
    n = length(f.coefficients)
    n == 0 && return zero(cfstype(f))
    n == 1 && return first(f.coefficients)
    foldr(-,coefficients(f,Chebyshev))
end

Base.last(f::Fun{Ultraspherical{Int,D,R}}) where {D,R} = reduce(+,coefficients(f,Chebyshev))

Base.first(f::Fun{Ultraspherical{O,D,R}}) where {O,D,R} = f(leftendpoint(domain(f)))
Base.last(f::Fun{Ultraspherical{O,D,R}}) where {O,D,R} = f(rightendpoint(domain(f)))

Fun(::typeof(identity), d::Ultraspherical) = Fun(Fun(identity, domain(d)),d)

## Calculus

spacescompatible(a::Ultraspherical,b::Ultraspherical) =
    order(a) == order(b) && domainscompatible(a,b)
hasfasttransform(::Ultraspherical) = true

include("UltrasphericalOperators.jl")
include("DirichletSpace.jl")
include("ContinuousSpace.jl")


text:
using DiffEqFlux, Flux, LinearAlgebra
using DiffEqNoiseProcess
using StochasticDiffEq
using Statistics
using DiffEqSensitivity
using DiffEqBase.EnsembleAnalysis
using Zygote

@testset "Neural SDE" begin
    function sys!(du, u, p, t)
        r, e, μ, h, ph, z, i = p
        du[1] = e * 0.5 * (5μ - u[1]) # nutrient input time series
        du[2] = e * 0.05 * (10μ - u[2]) # grazer density time series
        du[3] = 0.2 * exp(u[1]) - 0.05 * u[3] - r * u[3] / (h + u[3]) * u[4] # nutrient concentration
        du[4] =
            r * u[3] / (h + u[3]) * u[4] - 0.1 * u[4] -
            0.02 * u[4]^z / (ph^z + u[4]^z) * exp(u[2] / 2.0) + i #Algae density
    end

    function noise!(du, u, p, t)
        du[1] = p[end] # n
        du[2] = p[end] # n
        du[3] = 0.0
        du[4] = 0.0
    end

    datasize = 10
    tspan = (0.0f0, 3.0f0)
    tsteps = range(tspan[1], tspan[2], length = datasize)
    u0 = [1.0, 1.0, 1.0, 1.0]
    p_ = [1.1, 1.0, 0.0, 2.0, 1.0, 1.0, 1e-6, 1.0]

    prob = SDEProblem(sys!, noise!, u0, tspan, p_)
    ensembleprob = EnsembleProblem(prob)

    solution = solve(
        ensembleprob,
        SOSRI(),
        EnsembleThreads();
        trajectories = 1000,
        abstol = 1e-5,
        reltol = 1e-5,
        maxiters = 1e8,
        saveat = tsteps,
    )

    (truemean, truevar) = Array.(timeseries_steps_meanvar(solution))

    ann = FastChain(FastDense(4, 32, tanh), FastDense(32, 32, tanh), FastDense(32, 2))
    α = Float64.(initial_params(ann))

    function dudt_(du, u, p, t)
        r, e, μ, h, ph, z, i = p_

        MM = ann(u, p)

        du[1] = e * 0.5 * (5μ - u[1]) # nutrient input time series
        du[2] = e * 0.05 * (10μ - u[2]) # grazer density time series
        du[3] = 0.2 * exp(u[1]) - 0.05 * u[3] - MM[1] # nutrient concentration
        du[4] = MM[2] - 0.1 * u[4] - 0.02 * u[4]^z / (ph^z + u[4]^z) * exp(u[2] / 2.0) + i #Algae density
        return nothing
    end

    function dudt_op(u, p, t)
        r, e, μ, h, ph, z, i = p_

        MM = ann(u, p)

        [e * 0.5 * (5μ - u[1]), # nutrient input time series
         e * 0.05 * (10μ - u[2]), # grazer density time series
         0.2 * exp(u[1]) - 0.05 * u[3] - MM[1], # nutrient concentration
         MM[2] - 0.1 * u[4] - 0.02 * u[4]^z / (ph^z + u[4]^z) * exp(u[2] / 2.0) + i] #Algae density
    end

    function noise_(du, u, p, t)
        du[1] = p_[end]
        du[2] = p_[end]
        du[3] = 0.0
        du[4] = 0.0
        return nothing
    end

    function noise_op(u, p, t)
        [p_[end],
         p_[end],
         0.0,
         0.0]
    end

    prob_nn = SDEProblem(dudt_, noise_, u0, tspan, p = nothing)
    prob_nn_op = SDEProblem(dudt_op, noise_op, u0, tspan, p = nothing)

    function loss(θ)
        tmp_prob = remake(prob_nn, p = θ)
        ensembleprob = EnsembleProblem(tmp_prob)
        tmp_sol = Array(solve(
            ensembleprob,
            EM();
            dt = tsteps.step,
            trajectories = 100,
            sensealg = ReverseDiffAdjoint(),
           ))
        tmp_mean = mean(tmp_sol,dims=3)[:,:]
        tmp_var = var(tmp_sol,dims=3)[:,:]
        sum(abs2, truemean - tmp_mean) + 0.1 * sum(abs2, truevar - tmp_var), tmp_mean
    end

    function loss_op(θ)
        tmp_prob = remake(prob_nn_op, p = θ)
        ensembleprob = EnsembleProblem(tmp_prob)
        tmp_sol = Array(solve(
            ensembleprob,
            EM();
            dt = tsteps.step,
            trajectories = 100,
            sensealg = ReverseDiffAdjoint(),
           ))
        tmp_mean = mean(tmp_sol,dims=3)[:,:]
        tmp_var = var(tmp_sol,dims=3)[:,:]
        sum(abs2, truemean - tmp_mean) + 0.1 * sum(abs2, truevar - tmp_var), tmp_mean
    end

    losses = []
    callback(θ, l, pred) = begin
        push!(losses, l)
        if length(losses)%50 == 0
            println("Current loss after $(length(losses)) iterations: $(losses[end])")
        end
        false
    end

    res1 = DiffEqFlux.sciml_train(
        loss,
        α,
        ADAM(0.1),
        cb = callback,
        maxiters = 200,
    )

    res2 = DiffEqFlux.sciml_train(
        loss_op,
        α,
        ADAM(0.1),
        cb = callback,
        maxiters = 200,
    )
end

@testset "Adaptive neural SDE" begin
    x_size = 2 # Size of the spatial dimensions in the SDE
    v_size = 2 # Output size of the control

    # Define Neural Network for the control input
    input_size = x_size + 1 # size of the spatial dimensions PLUS one time dimensions
    nn_initial = Chain(Dense(input_size,v_size)) # The actual neural network
    p_nn, model = Flux.destructure(nn_initial)
    nn(x,p) = model(p)(x)

    # Define the right hand side of the SDE
    const_mat = zeros(Float64, (x_size, v_size))
    for i = 1:max(x_size,v_size)
        const_mat[i,i] = 1
    end

    function f!(du,u,p,t)
        MM = nn([u;t],p)
        du .= u + const_mat*MM
    end

    function g!(du,u,p,t)
        du .= false*u .+ sqrt(2*0.001)
    end

    # Define SDE problem
    u0 = vec(rand(Float64, (x_size,1)))
    tspan = (0.0, 1.0)
    ts = collect(0:0.1:1)
    prob = SDEProblem{true}(f!, g!, u0, tspan, p_nn)

    W = WienerProcess(0.0,0.0,0.0)
    probscalar = SDEProblem{true}(f!, g!, u0, tspan, p_nn, noise=W)

    # Defining the loss function
    function loss(pars, prob, alg)
        function prob_func(prob, i, repeat)
            # Prepare new initial state and remake the problem
            u0tmp = vec(rand(Float64,(x_size,1)))

            remake(prob, p = pars, u0 = u0tmp)
        end

        ensembleprob = EnsembleProblem(prob, prob_func = prob_func)

        _sol = solve(ensembleprob, alg, EnsembleThreads(), sensealg = BacksolveAdjoint(), saveat = ts, trajectories = 10,
                     abstol=1e-1, reltol=1e-1)
        A = convert(Array, _sol)
        sum(abs2, A .- 1), mean(A)
    end

    # Actually training/fitting the model
    losses = []
    callback(θ, l, pred) = begin
        push!(losses, l)
        if length(losses)%1 == 0
            println("Current loss after $(length(losses)) iterations: $(losses[end])")
        end
        false
    end
    res1 = @time DiffEqFlux.sciml_train((p) -> loss(p,probscalar, LambaEM()), p_nn, ADAM(0.1), cb = callback, maxiters=5)
    res2 = @time DiffEqFlux.sciml_train((p) -> loss(p,probscalar, SOSRI()), p_nn, ADAM(0.1), cb = callback, maxiters=5)

    res1 = @time DiffEqFlux.sciml_train((p) -> loss(p,prob, LambaEM()), p_nn, ADAM(0.1), cb = callback, maxiters=5)
    res2 = @time DiffEqFlux.sciml_train((p) -> loss(p,prob, SOSRA()), p_nn, ADAM(0.1), cb = callback, maxiters=5)
end


text:
"""
	Contains base types from grammar and the related functions

	The library assumes that the grammar structs have at least the following attributes:

	rules::Vector{Any}    # list of RHS of rules (subexpressions)
	types::Vector{Symbol} # list of LHS of rules (types, all symbols)
	isterminal::BitVector # whether rule i is terminal
	iseval::BitVector     # whether rule i is an eval rule
	bytype::Dict{Symbol,Vector{Int}}   # maps type to all rules of said type
	childtypes::Vector{Vector{Symbol}} # list of types of the children for each rule. Empty if terminal
"""

"""
isterminal(rule::Any, types::AbstractVector{Symbol})

Returns true if the rule is terminal, ie does not contain any of the types in the provided vector.
For example, :(x) is terminal, and :(1+1) is terminal, but :(Real + Real) is typically not.
"""
function isterminal(rule::Any, types::AbstractVector{Symbol})
if isa(rule, Expr)
    for arg in rule.args
	if !isterminal(arg, types)
	    return false
	end
    end
end
return rule ∉ types
end

"""
iseval(rule::Any)

Returns true if the rule is the special evaluate immediately function, i.e., _()
"""
iseval(rule) = false
iseval(rule::Expr) = (rule.head == :call && rule.args[1] == :_)

"""
get_childtypes(rule::Any, types::AbstractVector{Symbol})

Returns the child types of a production rule.
"""
function get_childtypes(rule::Any, types::AbstractVector{Symbol})
retval = Symbol[]
if isa(rule, Expr)
    for arg in rule.args
	append!(retval, get_childtypes(arg, types))
    end
elseif rule ∈ types
    push!(retval, rule)
end
return retval
end

abstract type Grammar end

Base.getindex(grammar::Grammar, typ::Symbol) = grammar.bytype[typ]

"""
   append!(grammar1::Grammar, grammar2::Grammar)

Appends the production rules of grammar2 to grammar1.
"""
function Base.append!(grammar1::Grammar, grammar2::Grammar)
    N = length(grammar1.rules)
    append!(grammar1.rules, grammar2.rules)
    append!(grammar1.types, grammar2.types)
    append!(grammar1.isterminal, grammar2.isterminal)
    append!(grammar1.iseval, grammar2.iseval)
    append!(grammar1.childtypes, copy.(grammar2.childtypes))
    for (s,v) in grammar2.bytype
        grammar1.bytype[s] = append!(get(grammar1.bytype, s, Int[]), N .+ v)
    end
    grammar1
end

"""
    nonterminals(grammar::Grammar)

Returns a list of nonterminals in the grammar.
"""
nonterminals(grammar::Grammar) = collect(keys(grammar.bytype))

"""
    return_type(grammar::Grammar, rule_index::Int)

Returns the type of the production rule at rule_index.
"""
return_type(grammar::Grammar, rule_index::Int) = grammar.types[rule_index]

"""
    child_types(grammar::Grammar, rule_index::Int)

Returns the types of the children (nonterminals) of the production rule at rule_index.
"""
child_types(grammar::Grammar, rule_index::Int) = grammar.childtypes[rule_index]

"""
    isterminal(grammar::Grammar, rule_index::Int)

Returns true if the production rule at rule_index is terminal, i.e., does not contain any nonterminal symbols.
"""
isterminal(grammar::Grammar, rule_index::Int) = grammar.isterminal[rule_index]

"""
    iscomplete(grammar::Grammar, node::RuleNode)

    Returns true if the expression given by the node is complete expression, i.e., all leaves are terminal symbols
"""
function iscomplete(grammar::Grammar, node::RuleNode)
	if isterminal(grammar, node)
		return true
	elseif isempty(node.children)
		# if not terminal but has children
		return false
	else
		return all([iscomplete(grammar, c) for c in node.children])
	end
end

"""
    iseval(grammar::Grammar)

Returns true if any production rules in grammar contain the special _() eval function.
"""
iseval(grammar::Grammar) = any(grammar.iseval)

"""
    iseval(grammar::Grammar, rule_index::Int)

Returns true if the production rule at rule_index contains the special _() eval function.
"""
iseval(grammar::Grammar, index::Int) = grammar.iseval[index]

"""
    nchildren(grammar::Grammar, rule_index::Int)

Returns the number of children (nonterminals) of the production rule at rule_index.
"""
nchildren(grammar::Grammar, rule_index::Int) = length(grammar.childtypes[rule_index])

"""
    max_arity(grammar::Grammar)

Returns the maximum arity (number of children) over all production rules in the grammar.
"""
max_arity(grammar::Grammar) = maximum(length(cs) for cs in grammar.childtypes)

"""
    return_types(grammar::Grammar, node::RuleNode)

Returns the return type in the production rule used by node.
"""
return_type(grammar::Grammar, node::RuleNode) = grammar.types[node.ind]

"""
    child_types(grammar::Grammar, node::RuleNode)

Returns the list of child types in the production rule used by node.
"""
child_types(grammar::Grammar, node::RuleNode) = grammar.childtypes[node.ind]

"""
    isterminal(grammar::Grammar, node::RuleNode)

Returns true if the production rule used by node is terminal, i.e., does not contain any nonterminal symbols.
"""
isterminal(grammar::Grammar, node::RuleNode) = grammar.isterminal[node.ind]

"""
    nchildren(grammar::Grammar, node::RuleNode)

Returns the number of children in the production rule used by node.
"""
nchildren(grammar::Grammar, node::RuleNode) = length(child_types(grammar, node))

"""
    contains_returntype(node::RuleNode, grammar::Grammar, sym::Symbol, maxdepth::Int=typemax(Int))

Returns true if the tree rooted at node contains at least one node at depth less than maxdepth
with the given return type.
"""
function contains_returntype(node::RuleNode, grammar::Grammar, sym::Symbol, maxdepth::Int=typemax(Int))
    maxdepth < 1 && return false
    if return_type(grammar, node) == sym
        return true
    end
    for c in node.children
        if contains_returntype(c, grammar, sym, maxdepth-1)
            return true
        end
    end
    return false
end

function Base.show(io::IO, grammar::Grammar)
	for i in eachindex(grammar.rules)
	    println(io, i, ": ", grammar.types[i], " = ", grammar.rules[i])
	end
end


text:
module BinanceF

import HTTP, SHA, JSON, Dates, Printf

# base URL of the Binance API
BINANCE_API_REST = "https://fapi.binance.com/"
BINANCE_API_TICKER = string(BINANCE_API_REST, "fapi/v1/ticker/")
BINANCE_API_DEPTH = string(BINANCE_API_REST, "fapi/v1/depth")
BINANCE_API_INFO = string(BINANCE_API_REST, "/fapi/v1/exchangeInfo")

BINANCE_API_KLINES = string(BINANCE_API_REST, "fapi/v1/klines")
BINANCE_API_USER_DATA_STREAM = string(BINANCE_API_REST, "fapi/v1/userDataStream")

BINANCE_API_WS = "wss://stream.binance.com:9443/ws/"
#BINANCE_API_STREAM = "wss://stream.binance.com:9443/stream/"

function apiKS()
    apiKey = get(ENV, "BINANCE_APIKEY", "")
    apiSecret = get(ENV, "BINANCE_SECRET", "")

    @assert apiKey != "" || apiSecret != "" "BINANCE_APIKEY/BINANCE_APISECRET should be present in the environment dictionary ENV"

    apiKey, apiSecret
end

function dict2Params(dict::Dict)
    params = ""
    for kv in dict
        params = string(params, "&$(kv[1])=$(kv[2])")
    end
    params[2:end]
end

# signing with apiKey and apiSecret
function timestamp()
    Int64(floor(Dates.datetime2unix(Dates.now(Dates.UTC)) * 1000))
end

function hmac(key::Vector{UInt8}, msg::Vector{UInt8}, hash, blocksize::Int=64)
    if length(key) > blocksize
        key = hash(key)
    end

    pad = blocksize - length(key)

    if pad > 0
        resize!(key, blocksize)
        key[end - pad + 1:end] = 0
    end

    o_key_pad = key .⊻ 0x5c
    i_key_pad = key .⊻ 0x36

    hash([o_key_pad; hash([i_key_pad; msg])])
end

function doSign(queryString, apiSecret)
    bytes2hex(hmac(Vector{UInt8}(apiSecret), Vector{UInt8}(queryString), SHA.sha256))
end

# function HTTP response 2 JSON
function r2j(response)
    JSON.parse(String(response))
end

##################### PUBLIC CALL's #####################

# Simple test if binance API is online
function ping()
    r = HTTP.request("GET", string(BINANCE_API_REST, "fapi/v1/ping"))
    r.status
end

# Binance servertime
function serverTime()
    r = HTTP.request("GET", string(BINANCE_API_REST, "fapi/v1/time"))
    r.status
    result = r2j(r.body)

    Dates.unix2datetime(result["serverTime"] / 1000)
end

function get24HR()
    r = HTTP.request("GET", string(BINANCE_API_TICKER, "24hr"))
    r2j(r.body)
end

function getDepth(symbol::String; limit=100) # 500(5), 1000(10)
    r = HTTP.request("GET", string(BINANCE_API_DEPTH, "?symbol=", symbol,"&limit=",limit))
    r2j(r.body)
end

function get24HR(symbol::String)
    r = HTTP.request("GET", string(BINANCE_API_TICKER, "24hr?symbol=", symbol))
    r2j(r.body)
end

function getAllPrices()
    r = HTTP.request("GET", string(BINANCE_API_TICKER, "allPrices"))
    r2j(r.body)
end

function getAllBookTickers()
    r = HTTP.request("GET", string(BINANCE_API_TICKER, "allBookTickers"))
    r2j(r.body)
end

function getExchangeInfo()
    r = HTTP.request("GET", string(BINANCE_API_TICKER, "allBookTickers"))
    r2j(r.body)
end

# binance get candlesticks/klines data
function getKlines(symbol; startDateTime=nothing, endDateTime=nothing, interval="1m")
    query = string("?symbol=", symbol, "&interval=", interval)

    if startDateTime != nothing && endDateTime != nothing
        startTime = Printf.@sprintf("%.0d",Dates.datetime2unix(startDateTime) * 1000)
        endTime = Printf.@sprintf("%.0d",Dates.datetime2unix(endDateTime) * 1000)
        query = string(query, "&startTime=", startTime, "&endTime=", endTime)
    end
    r = HTTP.request("GET", string(BINANCE_API_KLINES, query))
    r2j(r.body)
end

##################### SECURED CALL's NEEDS apiKey / apiSecret #####################
function openOrders()
    query = string("recvWindow=5000&timestamp=", timestamp())
    r = HTTP.request("GET", string(BINANCE_API_REST, "fapi/v1/openOrders?", query, "&signature=", doSign(query)), headers)
    r2j(r.body)
end

function cancelOrder(symbol,origClientOrderId)
    query = string("recvWindow=5000&timestamp=", timestamp(),"&symbol=", symbol,"&origClientOrderId=", origClientOrderId)
    r = HTTP.request("DELETE", string(BINANCE_API_REST, "fapi/v1/order?", query, "&signature=", doSign(query)), headers)
    r2j(r.body)
end

function createOrder(symbol::String, orderSide::String;
    quantity::Float64=0.0, orderType::String = "LIMIT",
    price::Float64=0.0, stopPrice::Float64=0.0,
    icebergQty::Float64=0.0, newClientOrderId::String="")

      if quantity <= 0.0
          error("Quantity cannot be <=0 for order type.")
      end

      println("$orderSide => $symbol q: $quantity, p: $price ")

      order = Dict("symbol"           => symbol,
                      "side"             => orderSide,
                      "type"             => orderType,
                      "quantity"         => Printf.@sprintf("%.8f", quantity),
                      "newOrderRespType" => "FULL",
                      "recvWindow"       => 10000)

      if newClientOrderId != ""
          order["newClientOrderId"] = newClientOrderId;
      end

      if orderType == "LIMIT" || orderType == "LIMIT_MAKER"
          if price <= 0.0
              error("Price cannot be <= 0 for order type.")
          end
          order["price"] =  Printf.@sprintf("%.8f", price)
      end

      if orderType == "STOP_LOSS" || orderType == "TAKE_PROFIT"
          if stopPrice <= 0.0
              error("StopPrice cannot be <= 0 for order type.")
          end
          order["stopPrice"] = Printf.@sprintf("%.8f", stopPrice)
      end

      if orderType == "STOP_LOSS_LIMIT" || orderType == "TAKE_PROFIT_LIMIT"
          if price <= 0.0 || stopPrice <= 0.0
              error("Price / StopPrice cannot be <= 0 for order type.")
          end
          order["price"] =  Printf.@sprintf("%.8f", price)
          order["stopPrice"] =  Printf.@sprintf("%.8f", stopPrice)
      end

      if orderType == "TAKE_PROFIT"
          if price <= 0.0 || stopPrice <= 0.0
              error("Price / StopPrice cannot be <= 0 for STOP_LOSS_LIMIT order type.")
          end
          order["price"] =  Printf.@sprintf("%.8f", price)
          order["stopPrice"] =  Printf.@sprintf("%.8f", stopPrice)
      end

      if orderType == "LIMIT"  || orderType == "STOP_LOSS_LIMIT" || orderType == "TAKE_PROFIT_LIMIT"
          order["timeInForce"] = "GTC"
      end

      order
  end

# account call contains balances
function account(apiKey::String, apiSecret::String)
    headers = Dict("X-MBX-APIKEY" => apiKey)

    query = string("recvWindow=5000&timestamp=", timestamp())

    r = HTTP.request("GET", string(BINANCE_API_REST, "fapi/v1/account?", query, "&signature=", doSign(query, apiSecret)), headers)

    if r.status != 200
        println(r)
        return status
    end

    return r2j(r.body)
end

function executeOrder(order::Dict, apiKey, apiSecret; execute=false)
    headers = Dict("X-MBX-APIKEY" => apiKey)
    query = string(dict2Params(order), "&timestamp=", timestamp())
    body = string(query, "&signature=", doSign(query, apiSecret))
    println(body)

    uri = "fapi/v1/order/test"
    if execute
        uri = "fapi/v1/order"
    end

    r = HTTP.request("POST", string(BINANCE_API_REST, uri), headers, body)
    r2j(r.body)
end

# returns default balances with amounts > 0
function balances(apiKey::String, apiSecret::String; balanceFilter = x -> parse(Float64, x["free"]) > 0.0 || parse(Float64, x["locked"]) > 0.0)
    acc = account(apiKey,apiSecret)
    balances = filter(balanceFilter, acc["balances"])
end

# Websockets functions

function wsFunction(channel::Channel, ws::String, symbol::String)
    HTTP.WebSockets.open(string(BINANCE_API_WS, lowercase(symbol), ws); verbose=false) do io
      while !eof(io);
        put!(channel, r2j(readavailable(io)))
    end
  end
end

function wsTradeAgg(channel::Channel, symbol::String)
    wsFunction(channel, "@aggTrade", symbol)
end

function wsTradeRaw(channel::Channel, symbol::String)
    wsFunction(channel, "@trade", symbol)
end

function wsDepth(channel::Channel, symbol::String; level=5)
    wsFunction(channel, string("@depth", level), symbol)
end

function wsDepthDiff(channel::Channel, symbol::String)
    wsFunction(channel, "@depth", symbol)
end

function wsTicker(channel::Channel, symbol::String)
    wsFunction(channel, "@ticker", symbol)
end

function wsTicker24Hr(channel::Channel)
    HTTP.WebSockets.open(string(BINANCE_API_WS, "!ticker@arr"); verbose=false) do io
      while !eof(io);
        put!(channel, r2j(readavailable(io)))
    end
  end
end

function wsKline(channel::Channel, symbol::String; interval="1m")
  #interval => 1m 3m 5m 15m 30m 1h 2h 4h 6h 8h 12h 1d 3d 1w 1M
    wsFunction(channel, string("@kline_", interval), symbol)
end

function wsKlineStreams(channel::Channel, symbols::Array, interval="1m")
  #interval => 1m 3m 5m 15m 30m 1h 2h 4h 6h 8h 12h 1d 3d 1w 1M
    allStreams = map(s -> string(lowercase(s), "@kline_", interval), symbols)
    error = false;
    while !error
        try
            HTTP.WebSockets.open(string(BINANCE_API_WS,join(allStreams, "/")); verbose=false) do io
            while !eof(io);
                put!(channel, String(readavailable(io)))
            end
      end
        catch e
            println(e)
            error=true;
            println("error occured bailing wsklinestreams !")
        end
    end
end

function wsKlineStreams(callback::Function, symbols::Array; interval="1m")
    #interval => 1m 3m 5m 15m 30m 1h 2h 4h 6h 8h 12h 1d 3d 1w 1M
      allStreams = map(s -> string(lowercase(s), "@kline_", interval), symbols)
      @async begin
        HTTP.WebSockets.open(string("wss://stream.binance.com:9443/ws/",join(allStreams, "/")); verbose=false) do io
            while !eof(io)
                    data = String(readavailable(io))
                    callback(data)
            end
        end
    end
end

function openUserData(apiKey)
    headers = Dict("X-MBX-APIKEY" => apiKey)
    r = HTTP.request("POST", BINANCE_API_USER_DATA_STREAM, headers)
    return r2j(r.body)["listenKey"]
end

function keepAlive(apiKey, listenKey)
    if length(listenKey) == 0
        return false
    end

    headers = Dict("X-MBX-APIKEY" => apiKey)
    body = string("listenKey=", listenKey)
    r = HTTP.request("PUT", BINANCE_API_USER_DATA_STREAM, headers, body)
    return true
end

function closeUserData(apiKey, listenKey)
    if length(listenKey) == 0
        return false
    end
    headers = Dict("X-MBX-APIKEY" => apiKey)
    body = string("listenKey=", listenKey)
    r = HTTP.request("DELETE", BINANCE_API_USER_DATA_STREAM, headers, body)
   return true
end

function wsUserData(channel::Channel, apiKey, listenKey; reconnect=true)

    function keepAlive()
        keepAlive(apiKey, listenKey)
    end

    Timer(keepAlive, 1800; interval = 1800)

    error = false;
    while !error
        try
            HTTP.WebSockets.open(string(Binance.BINANCE_API_WS, listenKey); verbose=false) do io
                while !eof(io);
                    put!(channel, r2j(readavailable(io)))
                end
            end
        catch x
            println(x)
            error = true;
        end
    end

    if reconnect
        wsUserData(channel, apikey, openUserData(apiKey))
    end

end

# helper
filterOnRegex(matcher, withDictArr; withKey="symbol") = filter(x -> match(Regex(matcher), x[withKey]) != nothing, withDictArr);

end


text:
export load_ts_file, load_wf_file, spk_read_channel, spk_read_header,
    spk_channel_labels, spk_read_par

import JSON

# ============================================================================ #
# SpkHeader type for internal use only
# ============================================================================ #
mutable struct SpkHeader
    nchan::UInt16
    offset::Vector{UInt64}
    file::AbstractString
end
# ============================================================================ #
function load_ts_file(ifile::AbstractString)
    ext = splitext(ifile)[2]
    if ext == ".txt"
        return float(open(readlines, ifile))
    elseif ext != ".ts"
        warn("Input file does not have the proper extension, results may vary")
    end

    open(ifile, "r") do io
        seekend(io)
        nbyte = position(io)
        nel = floor(Int64, nbyte / sizeof(Float64))
        seekstart(io)

        return read(io, Float64, nel)
    end
end
# ============================================================================ #
function load_wf_file(ifile::AbstractString)
    open(ifile, "r") do io
        row = read(io, UInt64)
        col = read(io, UInt64)
        x = read(io, Float64, row, col)

        return x
    end
end
# ============================================================================ #
function spk_read_channel(ifile::AbstractString, label::AbstractString)
    hdr = spk_read_header(ifile)
    labels = spk_channel_labels(hdr)

    label in labels || error("Given channel does not exist")

    kchan = findfirst(isequal(label), labels)

    open(ifile, "r") do io

        seek(io, hdr.offset[kchan] + 64)

        nts = read(io, UInt64)
        wf_pts = read(io, UInt64)
        nmrk = read(io, UInt64)

        data = Dict{Symbol, Any}()

        data[:ts] = Vector{Float64}(undef, nts)
        read!(io, data[:ts])

        if wf_pts > 0
            data[:wf] = Matrix{Float64}(undef, wf_pts, nts)
            read!(io, data[:wf])
        else
            data[:wf] = Matrix{Float64}(undef, 0, 0)
        end
        data[:mrk] = Vector{Float64}(undef, nmrk)
        read!(io, data[:mrk])

        return data
    end
end
# ============================================================================ #
function spk_read_header(ifile::AbstractString)
    isfile(ifile) || error("File \"$(ifile)\" does not exist")

    open(ifile, "r") do io
        nchan = read(io, UInt16)
        offset = Vector{UInt64}(undef, Int64(nchan))
        read!(io, offset)

        return SpkHeader(nchan, offset, ifile)
    end
end
# ============================================================================ #
function spk_channel_labels(ifile::AbstractString)
    return spk_channel_labels(spk_read_header(ifile))
end
# ---------------------------------------------------------------------------- #
function spk_channel_labels(hdr::SpkHeader)
    open(hdr.file, "r") do io
        labels = Vector{AbstractString}(undef, hdr.nchan)
        for k = 1:hdr.nchan
            seek(io, hdr.offset[k])
            tmp = Vector{UInt8}(undef, 64)
            read!(io, tmp)
            labels[k] = strip(String(tmp), '\0')
        end

        return labels
    end
end
# ============================================================================ #
function spk_read_par(ifile::String)
    file, ext = splitext(ifile)
    parfile = file * ".json"
    if !isfile(parfile)
        error("Failed to locate PAR/JSON file: \"$(parfile)\"")
    end
    return JSON.parse(read(parfile, String))
end
# ============================================================================ #


text:
export rerun

import Base: show

"""
    rerun(count::Int = -1)

Returns an Observable that mirrors the source Observable with the exception of an error.
If the source Observable calls error, this method will resubscribe to the source Observable for a maximum of `count`
resubscriptions (given as a number parameter) rather than propagating the error call.

# Arguments:
- `count::Int`: Number of retry attempts before failing. Optional. Default is `-1`.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from(1:3) |> safe() |> map(Int, (d) -> d > 1 ? error("Error") : d) |> rerun(3)

subscribe!(source, logger())
;

# output
[LogActor] Data: 1
[LogActor] Data: 1
[LogActor] Data: 1
[LogActor] Data: 1
[LogActor] Error: ErrorException("Error")
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`catch_error`](@ref), [`logger`](@ref), [`safe`](@ref)
"""
rerun(count::Int = -1) = RerunOperator(count)

struct RerunOperator <: InferableOperator
    count :: Int
end

function on_call!(::Type{L}, ::Type{L}, operator::RerunOperator, source) where L
    return proxy(L, source, RerunProxy(operator.count))
end

operator_right(operator::RerunOperator, ::Type{L}) where L = L

struct RerunProxy <: SourceProxy
    count :: Int
end

source_proxy!(::Type{L}, proxy::RerunProxy, source::S) where { L, S } = RerunSource{L, S}(proxy.count, source)

mutable struct RerunActorProps
    count        :: Int
    subscription :: Teardown

    RerunActorProps(count::Int) = new(count, voidTeardown)
end

struct RerunInnerActor{L, S, A} <: Actor{L}
    source :: S
    actor  :: A
    props  :: RerunActorProps
end

getsubscription(actor::RerunInnerActor)         = actor.props.subscription
setsubscription!(actor::RerunInnerActor, value) = actor.props.subscription = value

getcount(actor::RerunInnerActor)         = actor.props.count
setcount!(actor::RerunInnerActor, value) = actor.props.count = max(-1, value)

function on_next!(actor::RerunInnerActor{L}, data::L) where L
    next!(actor.actor, data)
end

function on_error!(actor::RerunInnerActor, err)
    unsubscribe!(getsubscription(actor))

    count = getcount(actor)
    if count === 0
        error!(actor.actor, err)
    else
        setcount!(actor, count - 1)
        setsubscription!(actor, subscribe!(actor.source, actor))
    end
end

function on_complete!(actor::RerunInnerActor)
    complete!(actor.actor)
end

struct RerunSource{L, S} <: Subscribable{L}
    count  :: Int
    source :: S
end

function on_subscribe!(source::RerunSource{L, S}, actor::A) where { L, S, A }
    inner = RerunInnerActor{L, S, A}(source.source, actor, RerunActorProps(source.count))

    setsubscription!(inner, subscribe!(source.source, inner))

    return RerunSubscription(inner)
end

struct RerunSubscription{I} <: Teardown
    inner :: I
end

as_teardown(::Type{<:RerunSubscription}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::RerunSubscription)
    return unsubscribe!(getsubscription(subscription.inner))
end

Base.show(io::IO, ::RerunOperator)             = print(io, "RerunOperator()")
Base.show(io::IO, ::RerunProxy)                = print(io, "RerunProxy()")
Base.show(io::IO, ::RerunSource{L})    where L = print(io, "RerunSource($L)")
Base.show(io::IO, ::RerunSubscription)         = print(io, "RerunSubscription()")


text:
"""
    fit(template::objective_function_template, data::Any, theta::Vector; estimation_method::String = "M", br_method::String = "implicit_trace", optim_method = LBFGS(), optim_options = Optim.Options())

Fit an [`objective_function_template`](@ref) on `data` using M-estimation ([keyword argument](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `estimation_method = "M"`; default) or RBM-estimation (reduced-bias M estimation; [Kosmidis & Lunardon, 2020](http://arxiv.org/abs/2001.03786); [keyword argument](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `estimation_method = "RBM"`). Bias reduction is either through the maximization of the bias-reducing penalized objective in Kosmidis & Lunardon (2020) ([keyword argument](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `br_method = "implicit_trace"`; default) or by subtracting an estimate of the bias from the M-estimates ([keyword argument](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `br_method = "explicit_trace"`). The bias-reducing penalty is constructed internally using automatic differentiation (using the [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package), and the bias estimate using a combination of automatic differentiation (using the [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package) and numerical differentiation (using the [FiniteDiff](https://github.com/JuliaDiff/FiniteDiff.jl) package).

The maximization of the objective or the penalized objective is done using the [**Optim**](https://github.com/JuliaNLSolvers/Optim.jl) package. Optimization methods and options can be supplied directly through the [keyword arguments](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `optim_method` and `optim_options`, respectively. `optim_options` expects an object of class `Optim.Options`. See the [Optim documentation](https://julianlsolvers.github.io/Optim.jl/stable/#user/config/#general-options) for more details on the available options.
"""
function fit(template::objective_function_template,
             data::Any,
             theta::Vector;
             estimation_method::String = "M",
             br_method::String = "implicit_trace",
             optim_method = LBFGS(),
             optim_options = Optim.Options())
    if (estimation_method == "M")
        br = false
    elseif (estimation_method == "RBM")
        if (br_method == "implicit_trace")
            br = true
        elseif (br_method == "explicit_trace")
            br = false
        else
            error(br_method, " is not a recognized bias-reduction method")
        end
    else
        error(estimation_method, " is not a recognized estimation method")
    end
    obj = beta -> -objective_function(beta, data, template, br)
    out = optimize(obj, theta, optim_method, optim_options)
    if (estimation_method == "M")
        theta = out.minimizer
    elseif (estimation_method == "RBM")
        if (br_method == "implicit_trace")
            theta = out.minimizer
        elseif (br_method == "explicit_trace")
            quants = obj_quantities(out.minimizer, data, template, true)
            jmat_inv = quants[2]
            adjustment = FiniteDiff.finite_difference_gradient(beta -> obj_quantities(beta, data, template, true)[1], out.minimizer)
            theta = out.minimizer + jmat_inv * adjustment
            br = true
        end
    end
    GEEBRA_results(out, theta, data, template, br, true, br_method)
end

"""
    fit(template::estimating_function_template, data::Any, theta::Vector; estimation_method::String = "M", br_method::String = "implicit_trace", nlsolve_arguments...)

Fit an [`estimating_function_template`](@ref) on `data` using M-estimation ([keyword argument](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `estimation_method = "M"`; default) or RBM-estimation (reduced-bias M estimation; [Kosmidis & Lunardon, 2020](http://arxiv.org/abs/2001.03786); [keyword argument](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `estimation_method = "RBM"`). Bias reduction is either through the solution of the empirically adjusted estimating functions in Kosmidis & Lunardon (2020) ([keyword argument](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `br_method = "implicit_trace"`; default) or by subtracting an estimate of the bias from the M-estimates ([keyword argument](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1) `br_method = "explicit_trace"`). The bias-reducing adjustments and the bias estimate are constructed internally using automatic differentiation (using the [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl) package).

The solution of the estimating equations or the adjusted estimating equations is done using the [**NLsolve**](https://github.com/JuliaNLSolvers/NLsolve.jl) package. Arguments can be passed directly to `NLsolve.nlsolve` through [keyword arguments](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1). See the [NLsolve README](https://github.com/JuliaNLSolvers/NLsolve.jl) for more information on available options.
"""
function fit(template::estimating_function_template,
             data::Any,
             theta::Vector;
             estimation_method::String = "M",
             br_method::String = "implicit_trace",
             nlsolve_arguments...)
    if (estimation_method == "M")
        br = false
    elseif (estimation_method == "RBM")
        if (br_method == "implicit_trace")
            br = true
        elseif (br_method == "explicit_trace")
            br = false
        else
            error(br_method, " is not a recognized bias-reduction method")
        end
    else
        error(estimation_method, " is not a recognized estimation method")
    end
    ef = get_estimating_function(data, template, br)
    out = nlsolve(ef, theta; nlsolve_arguments...)
    if (estimation_method == "M")
        theta = out.zero
    elseif (estimation_method == "RBM")
        if (br_method == "implicit_trace")
            theta = out.zero
        elseif (br_method == "explicit_trace")
            quants = ef_quantities(out.zero, data, template, true)
            adjustment = quants[1]
            jmat_inv = quants[2]
            theta = out.zero + jmat_inv * adjustment
            br = true
        end
    end
    GEEBRA_results(out, theta, data, template, br, false, br_method)
end


text:
export lens
export micro_lens_array

"""
    lens(f, x, y; λ=550e-9, n=1)

Return complex transmission value at position `x` and `y`
for a lens.
The lens has a focal length `f` and the surrounding medium has
refractive index `n`. Wavelength is `λ`.
"""
function lens(f::Number, x::Number, y::Number; radius=Inf, λ=550e-9, n=1)
    κ = calc_κ(λ, n)
    r2 = @. x ^ 2 + y ^ 2

    out = @. circ(radius, x, y) * exp(-1im * π * κ / f * (x^2 + y^2))
    return out
end

"""
    lens(f, L, size=(512, 512); λ=550e-9, n=1)

Return complex lens transmission value for a field of `size`
and a length of `L`
The lens has a focal length `f` and the surrounding medium has
refractive index `n`. Wavelength is `λ`.

"""
function lens(f::Number, L::Number, size::Tuple=(512, 512); radius=Inf, λ=550e-9, n=1)
    # initialize output array with complex numbers
    out = zeros(ComplexF64, size)

    # iterate over values
    for (j, x) in enumerate(fftpos(L, size[2]))
        for (i, y) in enumerate(fftpos(L, size[1]))
            out[i, j] = circ(radius, x, y) * lens(f, x, y; radius=radius, λ=λ, n=n)
        end
    end

    return out
end

"""
    micro_lens_array(f, L, pitch, size=(512, 512); centered=true, λ=550e-9, n=1, aperture=true,
                     dtype=ComplexF32)

Return a micro lens array (MLA) with pitch size `pitch` in the field size of `L`.
All lenslets have the same focal length `f`.
As default `centered=true` means that the there is a micro lens centered around
the center of the array.
`aperture=true` indicates that there is a circular aperture cutting the field.
Otherwise the aperture are square shaped of the same size as the `pitch`.
"""
function micro_lens_array(f, L, pitch, size=(512, 512); centered=true, λ=550e-9, n=1, aperture=true, dtype=ComplexF32)
    # pitch size in indices
    mla_size = round.(Int, (size .* pitch) ./  L)
    # ÷ doesn't convert to Int tuple
    mla_size = convert(Tuple{Int, Int}, mla_size)

    # if centered, we embed a odd number of micro lenses
    # into a larger arrays
    # at the end, we do a center extract of this larger array
    # otherwise it is troublesome to calculate the offset and
    # cut off micro lenses
    if centered
        N_mla = ceil(Int, L / pitch)
        if mod(N_mla, 2) == 0
            N_mla = N_mla + 1
        end

        size_out = size
        size = N_mla .* mla_size
    end

    out = zeros(dtype, size)
    # put the micro lenses in the out field
    for j = 1:mla_size[2]:size[2]
        for i = 1:mla_size[1]:size[1]
            # min is required in the case when a full micro lens
            # doesn't fit anymore
            mla_i = min(mla_size[1], size[1] - (i - 1))
            mla_j = min(mla_size[2], size[2] - (j - 1))

            out_i = min(size[1], i - 1 + mla_size[1])
            out_j = min(size[2], j - 1 + mla_size[2])

            # add small microlens
            if aperture
                r = pitch / 2
            else
                r = Inf
            end
            out[i:out_i, j:out_j] = lens(f, pitch, mla_size, λ=λ, n=n, radius=r)[1:mla_i, 1:mla_j]
        end
    end

    if centered
        out = center_extract(out, size_out)
    end

    return out
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
type Model
    ws::Vector
end

function (m::Model)(s::State)
    tokens = s.tokens
    s0 = tokens[s.top]
    s1 = isnull(s.left) ? tokens[s.left.top] : nulltoken
    b0 = s.right <= length(tokens) ? tokens[s.right] : nulltoken
    sl = isnull(s.lch) ? tokens[s.lch.top] : nulltoken
    sr = isnull(s.rch) ? tokens[s.rch.top] : nulltoken

end

function train()
    opt = SGD()
    for epoch = 1:nepochs
        println("epoch: $(epoch)")
        opt.rate = 0.0075 / epoch
        loss = fit(train_x, train_y, model, crossentropy, opt)
        println("loss: $(loss)")

        test_z = map(x -> predict(model,x), test_x)
        acc = accuracy(test_y, test_z)
        println("test acc.: $(acc)")
        println("")
    end
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
using CausalDiscovery.CISC: f

@testset "f" begin
  @test f(0.123) > 0
end


text:
function header(args...; size::Int = 1, kwargs...)
  1 <= size <= 6 || error("Invalid header size - expected 1:6")
  func = getproperty(Genie.Renderer.Html, Symbol("h$size"))
  func(class="text-h$size", args...; kwargs...)
end


text:
"""
   Delegates loading functions & traits for the HPF spectrograph
Author: Eric Ford
Created: Sept 2020
"""

"""
Module providing types and traits and customized functions for the HPF Spectrograph.
"""
module HPF
using RvSpectMLBase
using DataFrames, FITSIO

#type EXPRES <: AbstractInstrument end
struct HPF1D <: AbstractInstrument1D end
struct HPF2D <: AbstractInstrument2D end
const AnyHPF = Union{HPF1D,HPF2D}
export HPF1D, HPF2D, AnyHPF

#include("traits.jl")
export min_order, max_order, min_pixel_in_order, max_pixel_in_order
export orders_to_use_default, min_col_default, max_col_default
export metadata_symbols_default, metadata_strings_default

#include("io.jl")
#export read_data, read_solar_data
#export make_manifest,
# read_header not exported to avoid conflict with FITSIO.read_header

#import RvSpectMLBase: get_inst_module
get_inst_module(::AnyHPF) = HPF

end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using PowerSystems
using PowerSimulations
using JuMP

const PS = PowerSimulations

base_dir = string(dirname(dirname(pathof(PowerSystems))))
println(joinpath(base_dir,"data/data_5bus_uc.jl"))
include(joinpath(base_dir,"data/data_5bus_uc.jl"))

sys5 = PowerSystem(nodes5, generators5, loads5_DA, branches5, nothing, 100.0)

#Generator Active and Reactive Power Variables
@test try
    Net = PS.StandardAC
    m = Model()
    netinjection = PS.instantiate_network(Net, sys5)
    PS.constructdevice!(m, netinjection, RenewableGen, PS.RenewableCurtail, Net, sys5)
true finally end

#Cooper Plate and Dispatch
@test try
    Net = PS.CopperPlatePowerModel
    m = Model();
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, RenewableGen, PS.RenewableCurtail, Net, sys5);
true finally end

#PTDF Plate and Dispatch
@test try
    Net = PS.StandardPTDF
    m = Model();
    netinjection = PS.instantiate_network(Net, sys5);
    PS.constructdevice!(m, netinjection, RenewableGen, PS.RenewableCurtail, Net, sys5);
true finally end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -2.3024137037198598 <= q <= 2.756640266423038)

                     add_NL_constraint(m, :(swish(0.8745011162287222 + -0.4747792679508138*swish(0.20979439694353186 + 0.2523963619940135*$(x[1]) + 0.6039440929888729*$(x[2])) + 0.06822380173369114*swish(-0.192203947089018 + 0.35618645260012016*$(x[1]) + -0.9380253880003773*$(x[2]))) + swish(-0.5284568968024108 + 0.8782860084347535*swish(0.20979439694353186 + 0.2523963619940135*$(x[1]) + 0.6039440929888729*$(x[2])) + 0.990986450920027*swish(-0.192203947089018 + 0.35618645260012016*$(x[1]) + -0.9380253880003773*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module LinearRegression
    using Pluto

    include("linear_regression.jl")

    run_notebook() = Pluto.run(notebook=joinpath(@__DIR__, "../notebook/ice_cream_vs_whiskey.jl"))
end


text:
"""

```
grid3D(solute,solute_atoms,mddf_result,output_file; dmin=1.5, ddax=5.0, step=0.5)
```

This function builds the grid of the 3D density function and fills an array of
mutable structures of type Atom, containing the position of the atoms of
grid, the closest atom to that position, and distance.

`solute` is a `ComplexMixtuers.Selection`, defining the solute. `solute_atoms` is the corresponding
vector of `PDBTools.Atom`s, and `mddf_result` is the result of a `mddf_result` calculation with
the correspondign solute.

`dmin` and `dmax` define the range of distance where the density grid will be built, and `step`
defines how fine the grid must be. Be aware that fine grids involve usually a very large (hundreds
of thousands points).

All parameters can be provides as keyword parameters.

### Example

```julia-repl
julia> using ComplexMixtures, PDBTools

julia> pdb = readPDB("./system.pdb");

julia> R = ComplexMixtures.load("./results.json");

julia> protein = select(pdb,"protein");

julia> solute = ComplexMixtures.Selection(protein,nmols=1);

julia> grid = ComplexMixtures.grid3D(solute=solute, solute_atoms=protein, mddf_result=R, output_file="grid.pdb");

```

`grid` will contain a vector of `Atom`s with the information of the MDDF at each grid point, and the
same data will be written in the `grid.pdb` file. This PDB file can be opened in VMD, for example, and contain
in the `beta` field the contribution of each protein residue to the MDDF at each point in space relative to the
protein, and in the `occupancy` field the distance to the protein. Examples of how this information can be
visualized are provided in the user guide of `ComplexMixtures`.

"""
grid3D(;
    solute = nothing,
    solute_atoms = nothing,
    mddf_result = nothing,
    output_file = nothing,
    dmin = 1.5,
    dmax = 5.0,
    step = 0.5,
) = grid3D(
    solute,
    solute_atoms,
    mddf_result,
    output_file,
    dmin = dmin,
    dmax = dmax,
    step = step,
)

function grid3D(
    solute,
    solute_atoms,
    mddf_result,
    output_file;
    dmin = 1.5,
    dmax = 5.0,
    step = 0.5,
)

    if nothing in (solute, solute_atoms, mddf_result)
        error("grid3D requires `solute`, `solute_atoms` and `mddf_result` definitions.")
    end

    # Simple function to interpolate data
    interpolate(x₁, x₂, y₁, y₂, xₙ) = y₁ + (y₂ - y₁) / (x₂ - x₁) * (xₙ - x₁)

    # Maximum and minimum coordinates of the solute
    lims = maxmin(solute_atoms)
    n = @. ceil(Int, (lims.xlength + 2 * dmax) / step + 1)

    # Building the grid with the nearest solute atom information
    igrid = 0
    grid = PDBTools.Atom[]
    for ix = 1:n[1]
        x = lims.xmin[1] - dmax + step * (ix - 1)
        for iy = 1:n[2]
            y = lims.xmin[2] - dmax + step * (iy - 1)
            for iz = 1:n[3]
                z = lims.xmin[3] - dmax + step * (iz - 1)
                rgrid = -1
                iat, r = PDBTools.closest(x, y, z, solute_atoms)
                if (dmin < r < dmax)
                    if rgrid < 0 || r < rgrid
                        at = solute_atoms[iat]
                        # Get contribution of this atom to the MDDF
                        c = ComplexMixtures.contrib(
                            solute,
                            mddf_result.solute_atom,
                            [at.index],
                        )
                        # Interpolate c at the current distance
                        iright = findfirst(d -> d > r, mddf_result.d)
                        ileft = iright - 1
                        cᵣ = interpolate(
                            mddf_result.d[ileft],
                            mddf_result.d[iright],
                            c[ileft],
                            c[iright],
                            r,
                        )
                        if cᵣ > 0
                            gridpoint = Atom(
                                index = at.index,
                                index_pdb = at.index_pdb,
                                name = at.name,
                                chain = at.chain,
                                resname = at.resname,
                                resnum = at.resnum,
                                x = x,
                                y = y,
                                z = z,
                                occup = r,
                                beta = cᵣ,
                                model = at.model,
                                segname = at.segname,
                            )
                            if rgrid < 0
                                igrid += 1
                                push!(grid, gridpoint)
                            elseif r < rgrid
                                grid[igrid] = gridpoint
                            end
                            rgrid = r
                        end # cᵣ>0
                    end # rgrid
                end # dmin/dmax
            end #iz
        end #iy
    end #ix

    # Now will scale the density to be between 0 and 99.9 in the temperature
    # factor column, such that visualization is good enough
    bmin, bmax = +Inf, -Inf
    for gridpoint in grid
        bmin = min(bmin, gridpoint.beta)
        bmax = max(bmax, gridpoint.beta)
    end
    for gridpoint in grid
        gridpoint.beta = (gridpoint.beta - bmin) / (bmax - bmin)
    end

    if output_file != nothing
        writePDB(grid, output_file)
    end
    return grid
end


text:
module QPIterators

using StaticArrays
using LinearAlgebra
using ..RefShapes: manifdim, IntegRule, quadrature
using ..RefShapes: npts, param_coords, weights
using ..FElements: refshape, ndofsperel, _geometrycarrier
import ..FElements: bfun, bfungradpar, jacjac
using ..FESpaces: edofbfnum

# To do: Create the data structure for an element that has multiple degrees of
# freedom per entity. Also store the basic data for the scalar finite element
# (single degree of freedom per entity).
function __bfundata(fesp, qr)
    bfnum = edofbfnum(fesp)
    nedof = length(bfnum)
    pc = qr.param_coords
    w  =  qr.weights
    npts = qr.npts
    MDIM = manifdim(refshape(fesp.fe))
    # First construct vectors of vectors of basis function values and matrices
    # of base function gradients in parametric coordinates for the scalar
    # element
    scalNs = Vector{Float64}[];
    scalgradNps = Vector{Adjoint{Float64,SArray{Tuple{MDIM},Float64,1,MDIM}}}[];
    for j in 1:npts
        push!(scalNs, bfun(fesp.fe, pc[j,:]))
        push!(scalgradNps, bfungradpar(fesp.fe, pc[j,:]))
    end
    scalgradN = similar(scalgradNps[1])
    # Now extend it to the vector element
    Ns = Vector{Float64}[];
    gradNps = Vector{Adjoint{Float64,SArray{Tuple{MDIM},Float64,1,MDIM}}}[];
    for j in 1:npts
        N = [scalNs[j][bfnum[k]] for k in 1:nedof]
        gN = [scalgradNps[j][bfnum[k]] for k in 1:nedof]
        push!(Ns, N)
        push!(gradNps, gN)
    end
    gradN = similar(gradNps[1])
    # Now do this for the geometry carrier
    gcfe = _geometrycarrier(fesp.fe)
    geomscalNs = Vector{Float64}[];
    geomscalgradNps = Vector{Adjoint{Float64,SArray{Tuple{MDIM},Float64,1,MDIM}}}[];
    for j in 1:npts
        push!(geomscalNs, bfun(gcfe, pc[j,:]))
        push!(geomscalgradNps, bfungradpar(gcfe, pc[j,:]))
    end
    return (bfnum, scalNs, scalgradNps, scalgradN, Ns, gradNps, gradN, geomscalNs, geomscalgradNps)
end

"""
    QPIterator{FES, MDIM}

Type of quadrature-point iterator, parameterized by
- `FES`: the type of the finite element space,
- `MDIM`: the manifold dimension of the finite element.
"""
mutable struct QPIterator{FES, MDIM}
    fesp::FES
    _quadr::IntegRule
    _bfnum::Vector{Int64}
    _scalbfuns::Vector{Vector{Float64}}
    _scalbfungrad_ps::Vector{Vector{Adjoint{Float64,SArray{Tuple{MDIM},Float64,1,MDIM}}}}
    _scalbfungrads::Vector{Adjoint{Float64,SArray{Tuple{MDIM},Float64,1,MDIM}}}
    _bfuns::Vector{Vector{Float64}}
    _bfungrad_ps::Vector{Vector{Adjoint{Float64,SArray{Tuple{MDIM},Float64,1,MDIM}}}}
    _bfungrads::Vector{Adjoint{Float64,SArray{Tuple{MDIM},Float64,1,MDIM}}}
    _geomscalbfuns::Vector{Vector{Float64}}
    _geomscalbfungrad_ps::Vector{Vector{Adjoint{Float64,SArray{Tuple{MDIM},Float64,1,MDIM}}}}
    _pt::Int64
end

"""
    QPIterator(fesp::FES, quadraturesettings) where {FES}

Construct quadrature-point iterator by associating it with a finite element space and supplying quadrature rule settings.
"""
function QPIterator(fesp::FES, quadraturesettings) where {FES}
    _quadr = quadrature(refshape(fesp.fe), quadraturesettings)
    bfnum, scalNs, scalgradNps, scalgradN, Ns, gradNps, gradN, geomscalNs, geomscalgradNps = __bfundata(fesp, _quadr)
    _pt = 0
    return QPIterator{FES, manifdim(refshape(fesp.fe))}(fesp, _quadr, bfnum, scalNs, scalgradNps, scalgradN, Ns, gradNps, gradN, geomscalNs, geomscalgradNps, _pt)
end

"""
    Base.iterate(it::QPIterator, state = 1)

Advance a quadrature point iterator.
"""
function Base.iterate(it::QPIterator, state = 1)
    if state > npts(it._quadr)
        return nothing
    else
        return (_update!(it, state), state+1)
    end
end
Base.length(it::QPIterator)  = npts(it._quadr)

function _update!(it::QPIterator, state)
    it._pt = state
    return it
end

"""
    bfun(it::QPIterator)

Retrieve vector of basis function values for the current quadrature point.
"""
function bfun(it::QPIterator)
    return it._bfuns[it._pt]
end

"""
    bfungradpar(it::QPIterator)

Retrieve vector of basis function gradients with respect to the parametric coordinates for the current quadrature point.
"""
function bfungradpar(it::QPIterator)
    return it._bfungrad_ps[it._pt]
end

"""
    bfungrad(it::QPIterator, Jac)

Retrieve vector of basis function gradients with respect to spatial coordinates
for the current quadrature point.

The Jacobian matrix maps between vectors in the parametric space and the spatial
vectors.
"""
function bfungrad(it::QPIterator, Jac)
    for j in 1:length(it._scalbfungrads)
        it._scalbfungrads[j] = it._scalbfungrad_ps[it._pt][j] / Jac
    end
    for k in 1:length(it._bfungrads)
        it._bfungrads[k] = it._scalbfungrads[it._bfnum[k]]
    end
    return it._bfungrads
end

"""
    weight(it::QPIterator)

Retrieve weight of the current quadrature point.
"""
weight(it::QPIterator) = it._quadr.weights[it._pt]

end


text:
function time_series(fileName)
  label = basename(fileName)
  df = readtable(fileName)
  df[:date] = DateTime[DateTime(tStr) for tStr in df[:date]];
  return df;
end

function time_series_autospec(fileName)
  label = basename(fileName)
  df = readtable(fileName, separator='\t')
  df[:date] = DateTime[DateTime(tStr, "yyyy-mm-dd HH:MM:SS") for tStr in df[:Time]];
  return df
end

function single_spec(fileName; dataOnly=false)
  y, gainStep, t, m0, p_pt = parseDataFileBothRows(fileName)
  row_A = vec(y[:,1])
  row_B = vec(y[:,2])
  gainFactor_A = getGainFactor(gainStep, _rowA)
  gainFactor_B = getGainFactor(gainStep, _rowB)
  pixelGain_A = getPixelGain(t, _rowA, gainStep)
  pixelGain_B = getPixelGain(t, _rowB, gainStep)
  baseline_A = polyFit(row_A, baseline_model)
  baseline_B = polyFit(row_B, baseline_model)

  ionsPerSpectrum_A = (row_A .- baseline_A) ./ pixelGain_A / gainFactor_A * C
  ionsPerSpectrum_B = (row_B .- baseline_B) ./ pixelGain_B / gainFactor_A * C
  yMin_A = 4 * std(ionsPerSpectrum_A[80:180])
  yMin_B = 4 * std(ionsPerSpectrum_B[80:180])

  @show(t)
  @show(gainFactor_A)
  @show(gainFactor_B)
  @show(maximum(ionsPerSpectrum_A))
  @show(yMin_A)

  pI_A, pA_A = findPeaks(ionsPerSpectrum_A, 4, yMin_A)
  pI_B, pA_B = findPeaks(ionsPerSpectrum_B, 4, yMin_B)

  s = Spectrum(gainStep, t, m0, p_pt, row_A, row_B, pixelGain_A,
               pixelGain_B, baseline_A, baseline_B, ionsPerSpectrum_A,
               ionsPerSpectrum_B, pI_A, pI_B, pA_A, pA_B)
  return s
end

################################################################################
function getGainFactor(gainStep,row)
  fileName = ""
  gainFactor = 0.0
  if (row == _rowA)
    fileName = "gainFactor_SPACE_A.txt"
  elseif (row == _rowB)
    fileName = "gainFactor_SPACE_B.txt"
  end

  iFile = open(joinpath(homedir(), ".julia/v0.4/DFMS/InstrumentData", fileName), "r")
  i=0
  while !eof(iFile)
    line = readline(iFile)
    if (i == gainStep)
      gainFactor = parse(Float64, matchall(r"(\d+.\d+[eE][+-]\d+)", line)[1])
    end
    i += 1
  end
  close(iFile)
  return gainFactor
end

function get_gain_factor(gainstep, isFMdata)
  if isFMdata
    gain_factor_fm(gainstep)
  else
    gain_factor_pds(gainstep)
  end
end

function gain_factor_pds(gainstep)
  fileName = ""
  gainFactor = [0.0, 0.0]
  counter = 1
  gF = 0.0
  for myrow in [_rowA, _rowB]
    if (myrow == _rowA)
      fileName = "gainFactor_SPACE_A.txt"
    elseif (myrow == _rowB)
      fileName = "gainFactor_SPACE_B.txt"
    end

    iFile = open(joinpath(homedir(), ".julia/v0.4/DFMS/InstrumentData", fileName), "r")
    #iFile = open(joinpath("C:\\Users\\leroy\\.julia\\v0.4\\DFMS\\InstrumentData", fileName), "r")
    i=0
    while !eof(iFile)
      line = readline(iFile)
      if (i == gainstep)
        gF = parse(Float64, matchall(r"(\d+.\d+[eE][+-]\d+)", line)[1])
      end
      i += 1
    end
    close(iFile)
    gainFactor[counter] = gF
    counter += 1
  end
  return gainFactor
end

function gain_factor_fm(gainstep)
  gainFactor = [0.0, 0.0]
  fileName = "gainFactor_FM.txt"
  iFile = open(joinpath(homedir(), ".julia/v0.4/DFMS/InstrumentData", fileName), "r")
  i=1
  while !eof(iFile)
    line = readline(iFile)
    if i == gainstep
      px, volt, gf_a, gf_b = split(line, '\t')
      gf_a = parse(Float64, gf_a)
      gf_b = parse(Float64, gf_b)
      gainFactor[1] = gf_a
      gainFactor[2] = gf_b
    end
    i+=1
  end
  return gainFactor
end

function getGainFactor(gainStep)
  fileName = ""
  gainFactor = [0.0, 0.0]
  counter = 1
  gF = 0.0
  for myrow in [_rowA, _rowB]
    if (myrow == _rowA)
      fileName = "gainFactor_SPACE_A.txt"
    elseif (myrow == _rowB)
      fileName = "gainFactor_SPACE_B.txt"
    end

    iFile = open(joinpath(homedir(), ".julia/v0.4/DFMS/InstrumentData", fileName), "r")
    #iFile = open(joinpath("C:\\Users\\leroy\\.julia\\v0.4\\DFMS\\InstrumentData", fileName), "r")
    i=0
    while !eof(iFile)
      line = readline(iFile)
      if (i == gainStep)
        gF = parse(Float64, matchall(r"(\d+.\d+[eE][+-]\d+)", line)[1])
      end
      i += 1
    end
    close(iFile)
    gainFactor[counter] = gF
    counter += 1
  end
  return gainFactor
end

function parseDataFile(fileName, row=2)
  iFile = open(fileName, "r")
  i = 1
  y = zeros(Float64, 512)
  t = DateTime(2000)
  gainStep = 0
  m0 = 0
  p_pt = 0.0
  iSkip = 329
  isStartTimeFound = false
  while !eof(iFile)
    line = readline(iFile)
    if (contains(line, "START_TIME") & (isStartTimeFound == false))
      tStr = matchall(r"\d+-\d+-\d+T\d+:\d+:\d+", line)[1]
      t = DateTime(tStr)
      isStartTimeFound = true
    elseif contains(line, "ROSINA_DFMS_SCI_GAIN_OF_SPECTRUM")
      gainStep = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line , "ROSINA_DFMS_SCI_MASS")
      m0 = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line, "ROSINA_DFMS_SCI_P_PT_AT_PEAK")
      p_pt = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line, "SPICE_FILE_NAME")
      iSkip = 341
    end

    if (i >= iSkip)
      y[i-(iSkip-1)] = parse(Float64, matchall(r"(\d+)", line)[row])
    end
    i += 1
  end

  close(iFile)
  return gainStep, t, y, m0, p_pt
end

function get_dfms_data(fileName)
  fileType = fileName[end-2:end]
  if fileType == "DAT"
    return load_fm(fileName)
  elseif fileType == "TAB"
    return load_pds(fileName)
  else
    println("File type unknown:")
    @show(fileType)
  end
end

function load_fm(fileName)
  isFMdata = true
  i = 1
  y = zeros(Float64, 512, 2)
  t = DateTime(2000)
  gainStep = 0
  m0 = 0
  p_pt = 0.0
  iSkip = 130
  isStartTimeFound = false

  iFile = open(fileName, "r")
  while !eof(iFile)
    line = readline(iFile)
    if (contains(line, "Packet time") & (isStartTimeFound == false))
      tStr = matchall(r"\d+/\d+/\d+ \d+:\d+:\d+", line)[1]
      t = DateTime(tStr, "dd/mm/yyyy HH:MM:SS")
      isStartTimeFound = true
    elseif (contains(line, "Gain") & !contains(line, "(pre") & !contains(line, "@"))
      gainStep = round(Int, parse(Float64, matchall(r"(\d+.\d+e\+\d+)", line)[1]))
    elseif contains(line , "Mass     ")
      m0 = round(Int, parse(Float64, matchall(r"(\d+.\d+e\+\d+)", line)[1]))
    elseif contains(line, "p/p_T")
      p_pt = round(Int, parse(Float64, matchall(r"(\d+.\d+e\+\d+)", line)[1]))
    end

    if (i >= iSkip)
      y[i-(iSkip-1), 1] = parse(Float64, matchall(r"(\d+)", line)[2])
      y[i-(iSkip-1), 2] = parse(Float64, matchall(r"(\d+)", line)[3])
    end
    i += 1
  end
  close(iFile)
  return y, gainStep, t, m0, p_pt, isFMdata
end

function load_pds(fileName)
  isFMdata = false
  iFile = open(fileName, "r")
  i = 1
  y = zeros(Float64, 512, 2)
  t = DateTime(2000)
  gainStep = 0
  m0 = 0
  p_pt = 0.0
  iSkip = 329
  isStartTimeFound = false
  while !eof(iFile)
    line = readline(iFile)
    if (contains(line, "START_TIME") & (isStartTimeFound == false))
      tStr = matchall(r"\d+-\d+-\d+T\d+:\d+:\d+", line)[1]
      t = DateTime(tStr)
      isStartTimeFound = true
    elseif contains(line, "ROSINA_DFMS_SCI_GAIN_OF_SPECTRUM")
      gainStep = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line , "ROSINA_DFMS_SCI_MASS")
      m0 = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line, "ROSINA_DFMS_SCI_P_PT_AT_PEAK")
      p_pt = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line, "SPICE_FILE_NAME")
      iSkip = 341
    end

    if (i >= iSkip)
      y[i-(iSkip-1), 1] = parse(Float64, matchall(r"(\d+)", line)[2])
      y[i-(iSkip-1), 2] = parse(Float64, matchall(r"(\d+)", line)[3])
    end
    i += 1
  end

  close(iFile)
  return y, gainStep, t, m0, p_pt, isFMdata
end

function load_h5_data(dataset, y, row)
    gainStep = read(dataset, "gainStep")
    m0 = read(dataset, "commandedMass")
    p_pt = read(dataset, "p_pt")
    rows = read(dataset, "y")
    tStr = read(dataset, "tStart")
    t = DateTime(tStr)
    for k in 1:512
      y[k] = rows[k, row]
    end
    return gainStep, t, m0, p_pt
end

function update!(dataset, varName, var)
  if !has(dataset, varName)
    d_write(dataset, varName, var)
  else
    dataset[varName][:,:] = var
  end
end

function load_h5_data!(dataset, y)
    gainStep = read(dataset, "gainStep")
    m0 = read(dataset, "commandedMass")
    p_pt = read(dataset, "p_pt")
    rows = read(dataset, "y")
    tStr = read(dataset, "tStart")
    t = DateTime(tStr)
    for myRow = 1:2
      for pixel in 1:512
        y[pixel, myRow] = rows[pixel, myRow]
      end
    end
    return gainStep, t, m0, p_pt
end

function parseDataFileBothRows(fileName)
  iFile = open(fileName, "r")
  i = 1
  y = zeros(Int, 512, 2)
  t = DateTime(2000)
  gainStep = 0
  m0 = 0
  p_pt = 0.0
  iSkip = 329
  isStartTimeFound = false
  while !eof(iFile)
    line = readline(iFile)
    if (contains(line, "START_TIME") & (isStartTimeFound == false))
      tStr = matchall(r"\d+-\d+-\d+T\d+:\d+:\d+", line)[1]
      t = DateTime(tStr)
      isStartTimeFound = true
    elseif contains(line, "ROSINA_DFMS_SCI_GAIN_OF_SPECTRUM")
      gainStep = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line , "ROSINA_DFMS_SCI_MASS")
      m0 = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line, "ROSINA_DFMS_SCI_P_PT_AT_PEAK")
      p_pt = parse(Int, matchall(r"(\d+)", line)[1])
    elseif contains(line, "SPICE_FILE_NAME")
      iSkip = 341
    end

    if (i >= iSkip)
      y[i-(iSkip-1), 1] = parse(Int, matchall(r"(\d+)", line)[2])
      y[i-(iSkip-1), 2] = parse(Int, matchall(r"(\d+)", line)[3])
    end
    i += 1
  end

  close(iFile)
  return y, gainStep, t, m0, p_pt
end

function get_pixel_gain(t, gainstep, isFMdata)
  if isFMdata
    pixel_gain_fm(gainstep)
  else
    pixel_gain_pds(t, gainstep)
  end
end

function pixel_gain_fm(gainstep)
  pixelGain = zeros(Float64, 512, 2)
  fileName = joinpath(homedir(), ".julia/v0.4/DFMS/InstrumentData/pixelgain_FM/pixelgain_FM_2015_06.txt")
  iCol = (gainstep - 1) * 6 + 5
  data = readdlm(fileName, '\t', Float64)
  for i=1:size(data,1)
    pixelGain[i,1] = data[i,iCol]
    pixelGain[i,2] = data[i,iCol+1]
  end
  return pixelGain
end

function pixel_gain_pds(t, gainstep)
  fileName = ""
  i = 0
  path = joinpath(homedir(), ".julia/v0.4/DFMS/InstrumentData/pixelgain_SPACE")
  pixelGain = zeros(Float64, 512, 2)
  for row in [_rowA, _rowB]
    # pixel gains for july to 1 october
    if DateTime(2014,2,1) < t < DateTime(2014,10,1)
      fileName = joinpath(path, "pg_space_Aug2014.csv")
      i=2
    # pixel gains from 1 october to 31 december
    elseif DateTime(2014,10,1) <= t < DateTime(2015,1,1)
      if gainstep < 12
        fileName = joinpath(path, "pg_space_Nov2014_GS9_new.csv")
      elseif 12 <= gainstep < 15
        fileName = joinpath(path, "pg_space_Nov2014_GS14_new.csv")
      elseif gainstep >= 15
        fileName = joinpath(path, "pg_space_Nov2014_GS16_new.csv")
      end
      i=1
    # pixel gains from 1 january to april
    elseif DateTime(2015,1,1) <= t < DateTime(2016,4,1)
      if gainstep < 12
        fileName = joinpath(path, "pg_space_Feb2015_GS9_new.csv")
      elseif 12 <= gainstep < 15
        fileName = joinpath(path, "pg_space_Feb2015_GS14_new.csv")
      elseif gainstep >= 15
        fileName = joinpath(path, "pg_space_Feb2015_GS16_new.csv")
      end
      i=1
    end
    iFile = open(fileName, "r")
    while !eof(iFile)
      line = readline(iFile)
      pixelGain[i, row] = parse(Float64, matchall(r"(\d+.\d+)", line)[row+2])
      i += 1
    end
  end
  return pixelGain
end

function getPixelGain(t, gainStep)
  fileName = ""
  i = 0
  path = joinpath(homedir(), ".julia/v0.4/DFMS/InstrumentData/pixelgain_SPACE")
  pixelGain = zeros(Float64, 512, 2)
  for row in [_rowA, _rowB]
    # pixel gains for july to 1 october
    if DateTime(2014,2,1) < t < DateTime(2014,10,1)
      fileName = joinpath(path, "pg_space_Aug2014.csv")
      i=2
    # pixel gains from 1 october to 31 december
    elseif DateTime(2014,10,1) <= t < DateTime(2015,1,1)
      if gainStep < 12
        fileName = joinpath(path, "pg_space_Nov2014_GS9_new.csv")
      elseif 12 <= gainStep < 15
        fileName = joinpath(path, "pg_space_Nov2014_GS14_new.csv")
      elseif gainStep >= 15
        fileName = joinpath(path, "pg_space_Nov2014_GS16_new.csv")
      end
      i=1
    # pixel gains from 1 january to april
    elseif DateTime(2015,1,1) <= t < DateTime(2016,4,1)
      if gainStep < 12
        fileName = joinpath(path, "pg_space_Feb2015_GS9_new.csv")
      elseif 12 <= gainStep < 15
        fileName = joinpath(path, "pg_space_Feb2015_GS14_new.csv")
      elseif gainStep >= 15
        fileName = joinpath(path, "pg_space_Feb2015_GS16_new.csv")
      end
      i=1
    end
    iFile = open(fileName, "r")
    while !eof(iFile)
      line = readline(iFile)
      pixelGain[i, row] = parse(Float64, matchall(r"(\d+.\d+)", line)[row+2])
      i += 1
    end
  end
  return pixelGain
end


text:
function tdms_header!(TDMS::TDMSbuf, io::IO, v::Integer)
  fastskip(io, 4)
  TDMS.flags  = fastread(io, UInt32)
  fastskip(io, 4)
  TDMS.nsos   = fastread(io, UInt64) + LEAD_IN_LENGTH
  TDMS.rdos   = fastread(io, UInt64) + LEAD_IN_LENGTH
  TDMS.n_ch   = fastread(io, UInt32) - 0x00000002
  l_path      = fastread(io, UInt32)
  fastskip(io, 4 + l_path)
  nprops      = fastread(io, UInt32)

  for i = 1:nprops
    # Get variable name and code
    L         = fastread(io, UInt32)
    name      = String(fastread(io, L))

    # Add dictionaries as needed
    D = TDMS.hdr
    if occursin(".", name)
      k = String.(split(name, "."))
      L = length(k)
      for j in 1:L-1
        if haskey(D, k[j])
          D = D[k[j]]
        else
          D[k[j]] = Dict{String, Any}()
          D = D[k[j]]
        end
      end
      name = k[L]
    end
    pcode     = fastread(io, UInt32)
    (pcode == 0x00000000) && continue

    # Undefined/bad codes
    T         = get(tdms_codes, pcode, nothing)
    (T == nothing) && error("Undefined property type!")

    if pcode == 0x00000044
      # TDMS documentation has the order of these two variables reversed
      nf = fastread(io, UInt64)
      ns = fastread(io, Int64)

      dt = tdms_dtos + ns + nf*2^(-64)
      if name == "GPSTimeStamp"
        TDMS.ts = round(Int64, sμ*dt)
      else
        D[name] = u2d(dt)
        (v > 1) && println(name, " = ", D[name])
      end
      # Will be corrected to UTC by adding SystemInfomation.GPS.UTCOffset

    elseif pcode == 0x00000020
      # No documentation here either
      nn = fastread(io, UInt32)
      val = fastread(io, nn)
      if name == "name"
        TDMS.name = String(val)
      else
        D[name] = String(val)
        (v > 1) && println(name, " = ", D[name])
      end

    else
      if name == "SamplingFrequency[Hz]"
        TDMS.fs = fastread(io, T)
      elseif name == "Latitude"
        TDMS.oy = fastread(io, T)
      elseif name == "Longitude"
        TDMS.ox = fastread(io, T)
      elseif name == "Altitude"
        TDMS.oz = fastread(io, T)
      else
        D[name] = fastread(io, T)
        (v > 1) && println(name, " = ", D[name])
      end

      #=
       In which direction is StartPosition[m]?
                  what about Start Distance (m)?

      =#
    end
  end
  return nothing
end

function read_silixa_tdms(file::String, nn::String, s::TimeSpec, t::TimeSpec, chans::ChanSpec, memmap::Bool, v::Integer)

  io = memmap ? IOBuffer(Mmap.mmap(file)) : open(file, "r")
  tdms_header!(TDMS, io, v)

  # not documented; inferred from manufacturer's Matlab script
  # uses general formula: mask = 2^(n-1); N & mask != mask
  #                       N = UInt8(data_flags & 0b00100000)
  #                       n = 6 (6th bit = 1 if data are decimated)
  decimated = Bool(UInt8(TDMS.flags & DECIMATE_MASK) != 2^5)
  # I don't think this is actually part of the standard TDMS data spec

  # data header info
  fastskip(io, fastread(io, UInt32)+8) # Group information
  fastskip(io, fastread(io, UInt32)+4) # first channel path and length
  DataType = fastread(io, UInt32)
  fastskip(io, 4)
  chunk_size = Int64(fastread(io, UInt32))
  T = get(tdms_codes, DataType, nothing)
  (T == nothing) && error("Unsupported data type!")

  seg1_length   = Int64(div(div(TDMS.nsos - TDMS.rdos, TDMS.n_ch), sizeof(T)))
  if v > 0
    println("data type = ", T,
            ", chunk size = ", chunk_size,
            ", channel length = ", seg1_length,
            ", decimated = ", decimated)
  end

  # parse start time to get zero-indexed si, ei in each channel
  if (typeof(s) <: Real) && (typeof(t) <: Real)
    si = max(round(Int64, TDMS.fs*s), 0)
    ei = min(round(Int64, TDMS.fs*t), seg1_length) - 1
  else
    Δ = round(Int64, SeisIO.sμ / TDMS.fs)
    if s == "0001-01-01T00:00:00"
      s = u2d(TDMS.ts*μs)
    elseif isa(s, Real)
      s = u2d(TDMS.ts*μs + s)
    end
    if t == "9999-12-31T12:59:59"
      t = u2d(TDMS.ts*μs + seg1_length/TDMS.fs)
    elseif isa(t, Real)
      t = u2d(TDMS.ts*μs + t)
    end
    (d0, d1) = parsetimewin(s, t)
    t0 = DateTime(d0).instant.periods.value*1000 - dtconst
    t1 = DateTime(d1).instant.periods.value*1000 - dtconst
    (v > 2) && println("t0 = ", t0, ", t1 = ", t1)
    si = max(div(t0-TDMS.ts, Δ), 0)
    ei = min(div(t1-TDMS.ts, Δ), seg1_length) - 1
  end
  si = Int64(si)
  ei = Int64(ei)

  # this prints one-indexed forms
  (v > 1) && println("reading from si = ", si+1, " to ei = ", ei+1)

  # Chunk bounds
  first_chunk = div(si, chunk_size) + 1
  last_chunk = div(ei, chunk_size) + 1
  n_chunks = last_chunk - first_chunk + 1

  # Skip to start of data
  nskip = TDMS.rdos - position(io)
  if first_chunk > 1
    chunk_os = (first_chunk-1)*chunk_size*sizeof(T)*TDMS.n_ch
    nskip += chunk_os
    (v > 2) && println("To skip: ", chunk_os, " bytes of chunks.")
  end
  if nskip != 0
    (v > 0) && println("Skipping ", nskip, " bytes total.")
    fastskip(io, nskip)
  end

  # Read Data
  buf = Array{T, 2}(undef, TDMS.n_ch, chunk_size)
  data = Array{Float32, 2}(undef, ei-si+1, TDMS.n_ch)

  # ===================================================================
  # Read chunks
  # sj = starting index within chunk
  # ej = ending index within chunk
  j = 1
  jmax = div(seg1_length, chunk_size)*chunk_size - si
  for i in first_chunk:last_chunk
    if j > jmax
      if VERSION < v"1.4"
        buf = Array{T, 2}(undef, TDMS.n_ch, rem(seg1_length, chunk_size))
        read!(io, buf)
      else
        vbuf = view(buf, :, 1:rem(seg1_length, chunk_size))
        read!(io, vbuf)
      end
    else
      read!(io, buf)
    end
    sj = (i == first_chunk ? rem(si+1, chunk_size) : 1)
    ej = (i == last_chunk ? rem(ei, chunk_size) + 1 : chunk_size)
    nj = ej-sj+1
    vx = view(data, j:j+nj-1, :)
    vb = view(buf, :, sj:ej)
    (v > 2) && println("chunk #", i, "/", n_chunks, ", sj = ", sj, ", ej = ", ej, ", j = ", j, ":", j+nj-1)

    transpose!(vx, vb) # converts to Float32 if needed
    j += nj
  end

  # Done with file
  close(io)

  # ----------------------------------------------------------------
  # String values for :name, :id
  name = (try
    string(TDMS.hdr["SystemInfomation"]["OS"]["HostName"], "_",
           TDMS.hdr["SystemInfomation"]["Devices0"]["Model"], "_",
           TDMS.hdr["SystemInfomation"]["Devices1"]["Model"])
  catch
    ""
  end)
  net = nn * "."
  cha = string("..O", getbandcode(TDMS.fs), "0")

  # -----------------------------------------------------------------
  # Time values
  utc_os = get(TDMS.hdr, "SystemInfomation.GPS.UTCOffset", zero(Float64))
  ts = TDMS.ts + round(Int64, sμ*utc_os)

  # =================================================================
  # Parse to NodalData
  if isempty(chans)
    chans = 1:TDMS.n_ch
  end
  S = NodalData(data, TDMS.hdr, chans, ts)
  fill!(S.fs, TDMS.fs)
  fill!(S.src, realpath(file))
  fill!(S.units, "m/m")
  S.ox = TDMS.ox
  S.oy = TDMS.oy
  S.oz = TDMS.oz
  for (i,j) in enumerate(chans)
    S.id[i] = string(net, lpad(j, 5, '0'), cha)
    S.name[i] = string(name, "_", j)
  end
  # =================================================================

  return S
end