text:
"""
    FineAcquisitionResults

A struct that stores the fine acquisition fine acquisition
results for both the carrier and FFT based methods.

Fields:

- `prn::Int64`: PRN number processed
- `type::String`: fine acquisition method used (either `"fft"` or `"carrier"`)
- `fd_course::Float64`: course acquired Doppler frequency in Hz
- `fd_rate::Float64`: signal Doppler rate in Hz/s
- `n0_idx_course::Int64`: course acquired code start index
- `t_length::Float64`: coherent integration time in seconds
- `fd_fine::Float64`: fine Doppler frequency in Hz
- `fd_est::Float64`: (fine + course) Doppler frequency in Hz
- `phi_init::Float64`: initial phase estimate in rads
- `M::Int64`: number of integrations performed (carrier method only)
- `P::Array{Float64,2}`: 3x3 diagonal matrix containin initial uncertainties
                         of the phase, Doppler, and Doppler rate estimates
- `R::Array{Float64,1}`: measurement uncertainty in rads
"""
struct FineAcquisitionResults
    prn::Int64
    type::String
    fd_course::Float64
    fd_rate::Float64
    n0_idx_course::Int64
    t_length::Float64
    fd_fine::Float64
    fd_est::Float64
    phi_init::Float64
    M::Int64
    P::Array{Float64,2}
    R::Array{Float64,1}
end

"""
    fineacquisition(data::GNSSSignal, replica::ReplicaSignal, prn, fd_course,
                    n₀_idx_course, type::Val{:fft}; fd_rate=0.,
                    t_length=replica.t_length,
                    σω=10., err_bin_num_ϕ=1, err_bin_num_f=2,
                    M=1)

Performs an FFT based fine acquisition on `data`. Note that `t_length` must
equal `replica.t_length`. `data` can be either a `GNSSData` or `L5QSignal`
struct, however, `data` and `replica` must be two seperate structs.

Required Arguments:

- `data::GNSSSignal`: either `GNSSData` or `ReplicaSignal` struct
- `replica::ReplicaSignal`: struct to use for replica signal generation
    * `replica.t_length` determines the integration time, `T`, used in FFT fine
      acquisition
- `prn::Int`: PRN number to track
- `fd_course`: course acquired Doppler frequency in Hz
- `n₀_idx_course`: initial code start index location
- `type::Val{:fft}`: `Val(:fft)` uses the FFT based fine acquisition method

Optional Arguments:

- `fd_rate`: the Doppler frequency rate in Hz/s `(default = 0Hz)`
- `σω`: expected uncertainty of Doppler rate in Hz/s `(default = 1000Hz/s)`
- `err_bin_num_ϕ`: number of bins to to look around peak to estimate the phase
                   uncertainty `(default = 1)`
- `err_bin_num_f`: number of frequency bins equal to the frequency uncertainty
                   `(default = 2)`
    * filter σ is equal to `err_bin_num_f×(1/t_length)`
- `M`: the number of consecutive integrations and use the one with the highest SNR
    * `(default = 1)`

Returns:

- `FineAcquisitionResults` struct
"""
function fineacquisition(data::GNSSSignal, replica::ReplicaSignal, prn, fd_course,
                         n₀_idx_course, type::Val{:fft}; fd_rate=0.,
                         σω=1000., err_bin_num_ϕ=1, err_bin_num_f=0.25,
                         M=1)
    t_length=replica.t_length
    Δf = 1/replica.t_length
    N = replica.sample_num
    N_over_2 = Int64(floor(N/2))
    pk_vals = Array{Complex{Float64}}(undef, M)
    pk_valabs2s = Array{Float64}(undef, M)
    pk_idxs = Array{Int}(undef, M)
    PNs = Array{Float64}(undef, M)
    SNR_ests = Array{Float64}(undef, M)
    fd_fines = Array{Float64}(undef, M)
    fd_ests = Array{Float64}(undef, M)
    ϕ_inits = Array{Float64}(undef, M)
    ϕ_init_errs = Array{Float64}(undef, M)
    for m in 1:M
        # Generate replica
        # Set signal parameters
        definereplica!(replica;
                    prn=prn, f_d=fd_course,
                    fd_rate=fd_rate, phi=0., f_if=0.,
                    code_start_idx=n₀_idx_course,
                    include_carrier=false,
                    start_t=(m-1)*t_length)
        # Generate signal
        generatereplica!(replica)
        # Wipeoff IF and course Doppler from data and multiply by replica
        @threads for i in 1:replica.sample_num
            j = (m-1)*N + i
            t = calc_t_at_i(j, (m-1)*t_length, data.f_s)
            @inbounds replica.data[i] = data.data[j]*cis(-2π*(data.f_if+fd_course+0.5*fd_rate*t)*t)*replica.data[i]
        end
        # Perform in place FFT of `replica.data`
        fft!(replica.data)
        # Find peak within ±[x]kHz, where `x` is defined by freq_lim
        # From index 1 to N/2: positive frequencies
        # From index N/2 to N: negative frequencies
        # N_lim = Int64(floor(freq_lim/Δf))
        # pos_freq_interval = (1, N_lim)
        # neg_freq_interval = (N-N_lim, N)
        pk_val = 0. + 0im
        pk_valabs2 = 0.
        pk_idx = 1
        PN = 0.
        #### TODO: Make the peak searching loops below functions and multi-thread them. ####
        # Check positive frequencies for peak
        @inbounds for i in 1:N
            data_abs2 = abs2(replica.data[i])
            PN += data_abs2
            if data_abs2 > pk_valabs2
                pk_valabs2 = data_abs2
                pk_val = replica.data[i]
                pk_idx = i
            end
        end
        PN = (PN - pk_valabs2) / (N - 1)
        SNR_est = calc_snr(pk_valabs2, PN)
        # # Check negative frequencies for peak
        # @inbounds for i in neg_freq_interval[1]:neg_freq_interval[2]
        #     data_abs2 = abs2(replica.data[i])
        #     if data_abs2 > pk_valabs2
        #         pk_valabs2 = data_abs2
        #         pk_val = replica.data[i]
        #         pk_idx = i
        #     end
        # end
        # Calculate peak frequency (the Doppler frequency error)
        if pk_idx > N_over_2
            fd_fine = (pk_idx-N-1)*Δf
        else
            fd_fine = (pk_idx-1)*Δf
        end
        fd_est = fd_course + fd_fine
        # Calculate initial phase
        ϕ_init = atan(imag(pk_val), real(pk_val))
        # Calculate the covariance matrix
        # We estimate the error to be ±2 Doppler bin for the frequency error
        # and ±1 for the phase error.
        pk_low_idx = ((pk_idx-err_bin_num_ϕ)+replica.sample_num)%replica.sample_num
        pk_high_idx = ((pk_idx+err_bin_num_ϕ)+replica.sample_num)%replica.sample_num
        if pk_low_idx == 0
            pk_low_idx = replica.sample_num
        end
        if pk_high_idx == 0
            pk_high_idx = replica.sample_num
        end
        pk_low = replica.data[pk_low_idx]
        pk_high = replica.data[pk_high_idx]
        ϕ_low = atan(imag(pk_low)/real(pk_low))
        ϕ_high = atan(imag(pk_high)/real(pk_high))
        ϕ_init_err = mean([abs(ϕ_init-ϕ_low), abs(ϕ_init-ϕ_high)])

        pk_vals[m] = pk_val
        pk_valabs2s[m] = pk_valabs2
        pk_idxs[m] = pk_idx
        PNs[m] = PN
        SNR_ests[m] = SNR_est
        fd_fines[m] = fd_fine
        fd_ests[m] = fd_est
        ϕ_inits[m] = ϕ_init
        ϕ_init_errs[m] = ϕ_init_err
    end
    idx = argmax(SNR_ests)
    pk_val = pk_vals[idx]
    pk_valabs2 = pk_valabs2s[idx]
    pk_idx = pk_idxs[idx]
    PN = PNs[idx]
    SNR_est = SNR_ests[idx]
    fd_fine = fd_fines[idx]
    fd_est = fd_ests[idx]
    ϕ_init = ϕ_inits[idx]
    ϕ_init_err = ϕ_init_errs[idx]
    # replica.isreplica = false
    P = diagm([ϕ_init_err^2, (2π*err_bin_num_f*Δf)^2, (2π*σω)^2])
    R = [ϕ_init_err^2]
    # Return `FineAcquisitionResults` struct
    return FineAcquisitionResults(prn, String(:fft), fd_course, fd_rate, n₀_idx_course,
                                  t_length, fd_fine, fd_est, ϕ_init, M, P, R)
end

"""
    fineacquisition(data::GNSSSignal, replica::ReplicaSignal, prn, fd_course,
                    n₀_idx_course, type::Val{:carrier}; fd_rate=0.,
                    t_length=replica.t_length, freq_lim=50000., M=10,
                    σω=1000.)

Performs an carrier based fine acquisition on `data`.
`replica` decides what signal type to use and the length of each `M` segment.

`M` is the multiple of the length of replica to be used for fine acquisition.

Required Arguments:

- `data::GNSSSignal`: either `GNSSData` or `ReplicaSignal` struct
- `replica::ReplicaSignal`: struct to use for replica signal generation
    * `replica.t_length` determines the integration time, `T`
- `prn::Int`: PRN number to track
- `fd_course`: course acquired Doppler frequency in Hz
- `n₀_idx_course`: initial code start index location
- `type::Val{:fft}`: `Val(:fft)` uses the FFT based fine acquisition method

Optional Arguments:

- `fd_rate`: the Doppler frequency rate in Hz/s `(default = 0Hz)`
- `t_length`: `[NOT USED]` `(default = replica.t_length)`
- `M`: Number of integrations to perform `(default = 10)`
- `σω`: expected uncertainty of Doppler rate in Hz/s `(default = 1000Hz/s)`

Returns:

- `FineAcquisitionResults` struct
"""
function fineacquisition(data::GNSSSignal, replica::ReplicaSignal, prn, fd_course,
                         n₀_idx_course, type::Val{:carrier}; fd_rate=0.,
                         t_length=replica.t_length, M=10, σω=1000.)
    # Check that there are at least M replica signals that can fit in time
    # into data.
    Mblocks = floor(Int, data.sample_num/replica.sample_num)
    if Mblocks < M
        error("Data must be at least $(M*replica.t_length) seconds long.")
    end
    # Samples per `M` data segment
    N = replica.sample_num
    # Set ϕ_init and ϕ
    ϕ_init = 0.
    ϕ = Array{Float64}(undef, M)
    ϕi = 0.
    for i in 1:M
        # Generate replica
        # Set signal parameters
        definereplica!(replica;
                       prn=prn, f_d=fd_course,
                       fd_rate=fd_rate, phi=0, f_if=0.,
                       code_start_idx=n₀_idx_course,
                       include_carrier=false)
        # Generate signal
        generatereplica!(replica)
        # Get a `view` of the current data segment and its corresponding time array
        datasegment = view(data.data, (i-1)*N+1:i*N)
        # Wipeoff IF and course Doppler from data
        @threads for j in 1:replica.sample_num
            t = calc_t_at_i((i-1)*N + j, 0, data.f_s)
            @inbounds replica.data[j] = datasegment[j]*cis(-2π*(data.f_if+fd_course+0.5*fd_rate*t)*t)*replica.data[j]
        end
        # Perform in place fft operation
        fft!(replica.data)
        # Get FFT peak value
        pk_idx = argmax(abs2.(replica.data))
        # Take square root of `pk`
        pk = replica.data[pk_idx]
        # Calculate ϕ for current processing block
        ϕi = atan(imag(pk), real(pk))
        # ϕi = atan(imag(pk)/real(pk))
        @inbounds ϕ[i] = ϕi
        if i == 1
            @inbounds ϕ_init = ϕi
        end
    end
    # Take the difference between ϕs & check and correct phase values if > than π
    dϕ = diff(ϕ)
    for i in 1:M-1
        if dϕ[i] > π
            dϕ[i] -= 2π
        elseif dϕ[i] < -π
            dϕ[i] += 2π
        else
            # Do nothing. Within bounds.
        end
    end
    # Compute variance from average phase change
    # dϕvar = var(dϕ)
    dϕvar = abs2.(dϕ)
    # Find the maxumum variance and omit from fine Doppler fequency calculation
    maxvaridx = argmax(dϕvar)
    # Take the mean while ommitting the largest phase change in `dϕ`
    dϕavg = 0.
    for i in 1:M-1
        if i != maxvaridx
            dϕavg += dϕ[i]
        end
    end
    dϕavg /= (M-2)
    # Calculate the fine Doppler frequency
    f_s = replica.f_s
    # fd_fine = dϕavg/(2π*N/f_s)
    fd_fine = dϕavg/(2π*replica.t_length)
    fd_est = fd_course + fd_fine
    dϕ_init_err = std(([dϕ[1:maxvaridx-1]; dϕ[maxvaridx+1:end]] .- dϕavg))
    fd_err = dϕ_init_err/(2π*replica.t_length)
    # P = diagm([dϕ_init_err^2, 74.7*fd_err^2, σω^2])
    # P = diagm([dϕ_init_err^2, fd_err^2, σω^2])
    P = diagm([dϕ_init_err^2, (2π*replica.t_length)^2, (2π*σω)^2])
    R = [dϕ_init_err^2]
    # Return `FineAcquisitionResults` struct
    return FineAcquisitionResults(prn, String(:carrier), fd_course, fd_rate, n₀_idx_course,
                                  t_length, fd_fine, fd_est, ϕ_init, M, P, R)
end

"""
    fineacquisition(data::GNSSSignal, replica::ReplicaSignal, prn, fd_course,
                    n₀_idx_course, type::Val{:ffl}; fd_rate=0.,
                    t_length=replica.t_length, freq_lim=50000., M=10,
                    σω=1000.)

Performs an carrier based fine acquisition on `data`.
`replica` decides what signal type to use and the length of each `M` segment.

`M` is the multiple of the length of replica to be used for fine acquisition.

Required Arguments:

- `data::GNSSSignal`: either `GNSSData` or `ReplicaSignal` struct
- `replica::ReplicaSignal`: struct to use for replica signal generation
    * `replica.t_length` determines the integration time, `T`
- `prn::Int`: PRN number to track
- `fd_course`: course acquired Doppler frequency in Hz
- `n₀_idx_course`: initial code start index location
- `type::Val{:fll}`: `Val(:fll)` uses the FLL based fine acquisition method

Optional Arguments:

- `fd_rate`: the Doppler frequency rate in Hz/s `(default = 0Hz)`
- `t_length`: `[NOT USED]` `(default = replica.t_length)`
- `M`: Number of integrations to perform `(default = 10)`
- `σω`: expected uncertainty of Doppler rate in Hz/s `(default = 1000Hz/s)`
- `B`: FLL bandwith in Hz `(default = 20)`
- `dampening`: dampening ratio for FLL `(default = 0.707)`

Returns:

- `FineAcquisitionResults` struct
"""
function fineacquisition(data::GNSSSignal, replica::ReplicaSignal, prn, fd_course,
                         n₀_idx_course, type::Val{:fll}; fd_rate=0.,
                         t_length=replica.t_length, M=10, σω=1000.,
                         B=20, dampening=0.707, q_a=1, h₋₂=2e-20)
    ωₙ = B/0.53
    T = t_length
    sig_freq = replica.signal_type.sig_freq
	qω = 2*π^2*h₋₂  # oscillator frequency PSD
    # Check that there are at least M replica signals that can fit in time
    # into data.
    Mblocks = floor(Int, data.sample_num/replica.sample_num)
    if Mblocks < M
        error("Data must be at least $(M*replica.t_length) seconds long.")
    end
    # Samples per `M` data segment
    N = replica.sample_num
    # Initialize FLL
    Δω = Array{Float64}(undef, M-1)
    SNR = Array{Float64}(undef, M)
    x = Array{Float64}(undef, 2, M)
    x[:,1] .= [fd_course, fd_rate] .* 2π
    P = Array{Float64}(undef, 2, M)
    P[:,1] .= ([0.5/t_length, σω] .* 2π).^2
    x⁺ = [x[1,1]; x[2,1]]
    P⁺ = [P[1,1] 0; 0 P[2,1]]
    Δθₖ₋₁ = 0
    K = [2*dampening*ωₙ*T-3*ωₙ^2*T^2/2; ωₙ^2*T]
    K[2] = 0
    H = [1 T/2]
    A = [1 T; 0 1]
    Q = (2π*sig_freq)^2 .* [T*qω+T^3*q_a/(3*c^2) T^2*q_a/(2*c^2);
                            T^2*q_a/(2*c^2)      T*q_a/c^2]
    for i in 1:M
        # Generate replica
        # Set signal parameters
        f_d, fd_rate = x⁺ ./ 2π
        fd_rate = 0
        definereplica!(replica;
                       prn=prn, f_d=fd_course,
                       fd_rate=fd_rate, phi=0, f_if=0.,
                       code_start_idx=n₀_idx_course,
                       include_carrier=false,
                       start_t=(i-1)*T)
        # Generate signal
        generatereplica!(replica)
        # Get a `view` of the current data segment and its corresponding time array
        datasegment = view(data.data, (i-1)*N+1:i*N)
        # Wipeoff IF and course Doppler from data
        @threads for j in 1:replica.sample_num
            t = calc_t_at_i((i-1)*N + j, 0, data.f_s)
            @inbounds replica.data[j] = datasegment[j]*cis(-2π*(data.f_if+fd_course+0.5*fd_rate*t)*t)*replica.data[j]
        end
        # Perform in place fft operation
        fft!(replica.data)
        # Get FFT peak value
        pk_idx = argmax(abs2.(replica.data))
        # Take square root of `pk`
        pk = replica.data[pk_idx]
        # Calculate SNR
        PS = abs2(pk)
        PN = (sum(abs2.(replica.data)) - PS) / (N - 1)
        snr = calc_snr(PS, PN)
        SNR[i] = snr
        # Calculate ϕ for current processing block
        # ϕi = atan(imag(pk), real(pk))
        Δϕₖ = atan(imag(pk)/real(pk))
        if i == 1
            Δωₖ = 0
        else
            dϕ = (Δϕₖ - Δθₖ₋₁)
            if dϕ > π
                dϕ -= 2π
            elseif dϕ < -π
                dϕ += 2π
            else
                # Do nothing. Within bounds.
            end
            Δωₖ = dϕ/T
            Δω[i-1] = Δωₖ
        end
        Δθₖ₋₁ = Δϕₖ
        # x⁺ = x⁺ + K*Δωₖ
        x⁺ = x⁺ .+ Δωₖ
        # println(Δωₖ/2π)
        P⁺ = (I - K*H)*P⁺
        x[:,i] .= x⁺
        P[:,i] .= diag(P⁺)
        # Propagate states
        x⁺ = A*x⁺
        P⁺ = A*P⁺*A' + Q
    end
    return (x, P, A, Q, H, K, SNR, Δω)
    # Return `FineAcquisitionResults` struct
    # return FineAcquisitionResults(prn, String(:carrier), fd_course, fd_rate, n₀_idx_course,
    #                               t_length, fd_fine, fd_est, ϕ_init, "N/A", P, R)
end


text:
nt = 1000
seed = 99
m = ABCDQR_simula(nt,seed)
y = m["y"]
u = m["u"]
A = m["A"]
B = m["B"]
C = m["C"]
D = m["D"]
Q = m["Q"]
R = m["R"]
x10 = m["m1"]
P10 = m["P1"]

###########################
print("Testing ABCDQR_kfilter_s: ")
xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik = ABCDQR_kfilter_s(y,u,A,B,C,D,Q,R,x10)
if sum(xtt.^2) > 1e-6
    println("OK")
end

ext: jl


text:
export readcoverage_bam
export readcoverage_transcript_bam

"""
	readcoverage_bam(path_bam::String, chrom::String, leftpos::Int64, rightpos::Int64;
					output_prefix::String = "")

Calculates a read coverage for a genomic interval.

# Arguments
---------
- `path_bam::String`: the path to a BAM file.
- `chrom::String`: the chromosome name for a genomic interval.
- `leftpos::Int64`: the left position for a genomic interval.
- `rightpos::Int64`: the right position for a genomic interval.
- `output_prefix::String`: the prefix for output files. If this keyword is not specified or set to "" (defalut), no output files are saved.
"""
function readcoverage_bam(path_bam::String, chrom::String, leftpos::Int64, rightpos::Int64; output_prefix::String = "")
	if output_prefix != ""
		path_out = output_prefix * ".readCoverage.txt"
		path_out_plot = output_prefix * ".readCoverage.pdf"
		println(@sprintf "Start calculating read coverage...\n- bam: %s\n- output: %s\n          %s\n" path_bam path_out path_out_plot)
	end

	# File check
	if !isfile(path_bam)
		error(@sprintf "No such file: %s\n" path_bam)
	end
	if output_prefix != ""
		if !uv_access_writable(dirname(output_prefix))
			error(@sprintf "Output files are not writable to: %s\n" (dirname(output_prefix) != "" ? dirname(output_prefix) : "."))
		end
	end

	# Define the output read coverage container
	coverage = zeros(Int, rightpos - leftpos + 1)

	# Caclutate the read coverage
	open(BAM.Reader, path_bam, index=path_bam*".bai") do bam_reader

		# Check wheteher the region exists in BAM
		if ! exists_region_bam(bam_reader, chrom, leftpos, rightpos)
			error("The region does not exist in the BAM file.")
		end

        coverage = readcoverage_bam_base(bam_reader, chrom, leftpos, rightpos)
	end

	if output_prefix != ""

		# Write output (TSV) (Position, Coverage)
		println("Write results...")
		bin_numbers = collect(leftpos:rightpos)
		open(path_out, "w") do io
			writedlm(io, [bin_numbers coverage], '\t')
		end

		# Save plot
		plot_read_coverage(coverage, out_path = path_out_plot)

		# Message
		println(@sprintf "Finished! Check output files:\n- %s\n- %s" path_out path_out_plot)
	end

	return(coverage)
end

"""
	readcoverage_bam_base(bam_reader::BAM.Reader, chrom::String, leftpos::Int64, rightpos::Int64)

Generates the read coverage for a genomic interval.
Both `leftpos` and `rightpos` coordinates are assumed to be 1-based.
"""
function readcoverage_bam_base(bam_reader::BAM.Reader, chrom::String, leftpos::Int64, rightpos::Int64)
	# Define output read coverage container
	coverage = zeros(Int, rightpos - leftpos + 1)

	try
		# Evaluate each BAM record
		for record in eachoverlap(bam_reader, chrom, leftpos:rightpos)

			# Skip a BAM record for an unmapped read (if pileupread.alignment.is_unmapped:continue)
			if ! BAM.ismapped(record)
				continue
			end

			# Compatibility with RSeQC genebody_coverage.py
			## if pileupread.alignment.is_secondary:continue
			if BAM.flag(record) & SAM.FLAG_SECONDARY != 0
				continue
			end

			## if pileupread.alignment.is_qcfail:continue
			if BAM.flag(record) & SAM.FLAG_QCFAIL != 0
				continue
			end

			## if pileupread.alignment.is_duplicate:continue
			if BAM.flag(record) & SAM.FLAG_DUP != 0
				continue
			end

			read_leftpos = BAM.position(record)
			cigarRle = BAM.cigar_rle(record)
			offset = 0

			# Decode CIGAR string to get coverage of split-aligned reads
			for i in 1:length(cigarRle[1])
				if cigarRle[1][i] == OP_MATCH
					if leftpos <= read_leftpos + offset <= rightpos
						for j in 1:min(cigarRle[2][i], rightpos - read_leftpos - offset + 1)
							coverage[read_leftpos + offset + j - leftpos] += 1
						end
					elseif read_leftpos + offset < leftpos
						for j in (leftpos-(read_leftpos + offset)+1):min(cigarRle[2][i], rightpos - read_leftpos - offset + 1)
							coverage[read_leftpos + offset + j - leftpos] += 1
						end
					end
				end
				offset += cigarRle[2][i]
			end
		end
	catch e
		println("Catched and ignored an error: ", typeof(e))
	end
	return(coverage)
end

"""
	readcoverage_transcript_bam(bam_reader::BAM.Reader, t::BED.Record)

Returns read coverage for a transcript.
"""
function readcoverage_transcript_bam(bam_reader::BAM.Reader, t::BED.Record)

	# Check if the region of the transcript exists in BAM
	if ! exists_region_bam(bam_reader, BED.chrom(t), BED.chromstart(t), BED.chromend(t))
		error("The transcript region does not exist in the BAM file.")
	end

	# Get block_sizes (= exon lengths) and block_starts (=start position of exons relative to chromStart of t)
	block_sizes = BED.blocksizes(t)
	transcript_length = sum(block_sizes)
	block_starts = BED.blockstarts(t)

	coverage = zeros(Float64, transcript_length)

	# Get the read coverage on each exon and save it to `cov`
	# Note: Read coverage is saved in `cov` in a 'Left justified' manner,
	#       and left is 5'-end.
	# If the strand of a transcript is '-' STRAND_NEG, read coverage was reversed.
	p_start = 1
	p_end = 0
	s_g = BED.chromstart(t) + 1
	e_g = 0

	for i in 1:length(block_sizes) # For each exon

		# Calculate start and end positions of an exon on the genome
		s_g = block_starts[i] + s_g
		e_g = s_g + block_sizes[i] - 1

		# Calculate start and end positions of an exon on mRNA
		p_start = p_end + 1
		p_end = p_start + block_sizes[i] - 1

		# Get read coverage and save it
		if BED.strand(t) == STRAND_POS
			pos_cov = p_start:p_end
		elseif BED.strand(t) == STRAND_NEG
			pos_cov = (transcript_length-p_start+1):-1:(transcript_length-p_end+1)
		else # STRAND_NA, STRAND_BOTH
			continue
		end

		coverage[pos_cov] = readcoverage_bam_base(bam_reader, BED.chrom(t), s_g, e_g)
	end

	return(coverage)
end

"""
For test
"""
function bamToCoverage_cigarAware_position(bam_reader::BAM.Reader, chrom::String, pos::Int64)
	coverage = 0

	# Evaluate each BAM record
	for record in eachoverlap(bam_reader, chrom, pos:pos)
		# Skip a BAM record for an unmapped read
		if ! BAM.ismapped(record)
			continue
		end

		read_leftpos = BAM.position(record)
		cigarRle = BAM.cigar_rle(record)
		offset = 0

		# Decode CIGAR string to get coverage of split-aligned reads
		for i in 1:length(cigarRle[1])
			if cigarRle[1][i] == OP_MATCH
				if read_leftpos + offset <= pos && pos <= read_leftpos + offset + cigarRle[2][i] - 1
					coverage += 1
				end
			end
			offset += cigarRle[2][i]
		end
	end

	return(coverage)
end

"""
	exists_region_bam(bam_reader::BAM.Reader, chrom::String, leftpos::Int64, rightpos::Int64)

Returns true if a region chrom:leftpos-rightpos exists in a BAM file for which a reader is opened.
"""
function exists_region_bam(bam_reader::BAM.Reader, chrom::String, leftpos::Int64, rightpos::Int64)::Bool
	if in(chrom, bam_reader.refseqnames) &&
			leftpos <= rightpos &&
			leftpos >= 0 &&
			rightpos <= bam_reader.refseqlens[findfirst(isequal(chrom), bam_reader.refseqnames)]
		return true
	else
		return false
	end
end


text:
using CartesianFDM
using LinearAlgebra
using SparseArrays

n = (33, 32)
top = (nonperiodic(), periodic())

ctx = cartesianfdmcontext(top, n)

# square domain (0, 1) x (0, 1)
h = spacing.(top, n)
x = coordinate.(top, n)

# cylinder centered at (0.5, 0.5) with radius 0.25
cylinder(x, y, z...) = 0.25 - (2x - 1) ^ 2 - (2y - 1) ^ 2

# capacities
using Vofinit

h₀ = Cdouble.((h..., sum(h) / 2))
xₑ = Cdouble[0, 0, 0, 0]

using Base.Iterators

# volume
V = map(zip(CartesianIndices(n), product(x...))) do (ind, el)
    x₀ = Cdouble.((el..., 0.))
    prod(h) * getcc(cylinder, x₀, h₀, xₑ)
end

V = map(mask(ctx, reshape(V, :))) do val
    iszero(val) ? eps(val) : val
end

# area
A₁ = map(zip(CartesianIndices(n), product(x...))) do (ind, el)
    x₀ = Cdouble.((el..., 0.))
    h[2] * getcc(x₀, h₀, xₑ) do y...
        cylinder(x₀[1], y[2])
    end
end

A₂ = map(zip(CartesianIndices(n), product(x...))) do (ind, el)
    x₀ = Cdouble.((el..., 0.))
    h[1] * getcc(x₀, h₀, xₑ) do y...
        cylinder(y[1], x₀[2])
    end
end

A = mask(ctx, [reshape(A₁, :), reshape(A₂, :)])

# potential flow
W = reshape(map(CartesianIndices(n)) do ind
                ind[1] ≤ 2 || ind[1] ≥ n[1] - 3 ? 1. : 0.
            end, :)
_, U = potentialflow(ctx, A, V, W)

using Plots

heatmap(reshape(U[1], n...))

#=
T = scalar(:T, n)
D = ones(prod(n))

# secondary
G = gradient(dir, ctx, A, V, T, D)
L = divergence(dir, ctx, A, G)

#
M = map(V) do el
    el ≤ 10eps(el) ? prod(h) : el
end
=#

filename = "heat.jl"

# if file already exists
#include(filename)
# otherwise
#(heat, heat!) = build_function(L ./ 2M, T)
#open(filename, "w") do file
#    write(file, string(heat))
#    write(file, "\n")
#    write(file, string(heat!))
#end

#=
using DifferentialEquations

T₀ = zeros(prod(n))
tspan = (0.0, 1.0)
prob = ODEProblem(heat, T₀, tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

using Plots

heatmap(reshape(sol(0.1), n...))
=#

#=
potential = eval(first(build_function(Δ, Φ)))

using NLsolve

sol = nlsolve(zeros(prod(n))) do res, Φ
    res .= potential(Φ)
end
phi = getproperty(sol, :zero)
=#

#J = Symbolics.jacobian(Δ, Φ)
#J = eval(first(build_function(sparse(Symbolics.jacobian(Δ, Φ)))))()

#phi = rand(prod(n))
#res1 = eval(potential)(phi)
#res2 = J * phi + eval(potential)(zeros(prod(n)))

#@assert norm(res1 - res2) < 1e-14

# also in CYLINDER!
#=
vec = diag(J)
broadcast!(vec, vec) do x
    iszero(x) ? -8one(x) : zero(x)
end

J .+= vec

phi = (-J) \ eval(potential)(zeros(prod(n)))
=#

#=
@variables t p
(rhs, rhs!) = build_function(L, T, p, t)

using DifferentialEquations

tspan = (0.0, 100.0)
T₀ = zeros(prod(n))
prob = ODEProblem(eval(rhs!), T₀, tspan)
sol = solve(prob)

using Plots

heatmap(reshape(sol(0.1), n...))
=#

### Step 1 -- thermal flow
#
# 1. Potential flow ;
# 2. Add convection ;
# 3. Neumann B. C. on the outflow ;
# 4. Robin B. C.
#
### Step 2 -- chemically reacting flow
#
# 1. More scalars with Fick's law ;
# 2. Heterogeneous reactions (with Robin B. C.) ;
# 3. Homogeneous reactions ;
# 4. Fluid-solid interaction with quasi-steady assumption.
#
### Step 2bis -- coupling with Navier-Stokes solver
#
### Step 3 -- senstivity analysis


text:
using Seismic
using Base.Test

# test of Conjugate Gradients
nd = 100
nm = 300
L = 10*randn(nd,nm)
m = 10.0*rand(nm)
d = zeros(nd)
d = L*m
d = d + 2.0*randn(nd)

println(norm(d))
println(norm(m))

m1,cost = ConjugateGradients(d,[MatrixMultiplyOp],[Dict(:matrix=>L)],Niter=20,mu=0.2)

ext = Seismic.Extent(size(d,1),size(d,2),1,1,1,
0,0,0,1,1,
1,1,1,1,1,
"","","","","",
"","","","","",
"")
h = Header[]
for ix = 1:size(d,2)
		push!(h,Seismic.InitSeisHeader())
		h[ix] = Seismic.InitSeisHeader()
		h[ix].tracenum = ix
		h[ix].n1 = size(d,1)
		h[ix].d1 = 1
end

SeisWrite("tmp_d.seis",d,h,ext)

ConjugateGradients("tmp_m.seis","tmp_d.seis",[MatrixMultiplyOp],[Dict(:matrix=>L)],"tmp_cost.txt",Niter=20,mu=0.2)
m2,ext = SeisRead("tmp_m.seis")

# test that quality factor between disk and memory based CG
# is greater than 50 Decibels
quality_factor = 10*log10(norm(m1[:],2)/norm(m2[:]-m1[:],2))
println("Quality factor = ",quality_factor)
@test quality_factor > 10.

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
import Rocket as R

function stateful_flatmap(::Type{ElemT}, fn::Function, seed::SeedT) where { ElemT, SeedT }
    scan_fn(data, (_, state)) = fn(data, state)

    R.scan(Tuple{Vector{<:ElemT}, SeedT}, scan_fn, (ElemT[], seed)) |>
        R.concat_map(ElemT, ((val, _),) -> R.from(val))
end

max_stars_repo_name: miksu/FluidIntonation

max_issues_repo_name: miksu/FluidIntonation

max_forks_repo_name: miksu/FluidIntonation


text:
module UpSampling2DM
    using LoopVectorization

    mutable struct UpSampling2D
        save_layer::Any
        load_layer::Any
        activate::Any
        initialize::Any
        update::Any

        input_shape::Tuple
        output_shape::Tuple
        size::Tuple{Int64, Int64}
        activation_function::Module

        value::Array{Float32}
        output::Array{Float32}
        ∇biases::Array{Float32}
        δ::Array{Float32}

        function UpSampling2D(;input_shape::Tuple, size::Tuple{Int64, Int64}, activation_function::Module)
            output_shape = (input_shape[1]*size[1], input_shape[2]*size[2], input_shape[3])
            new(save_UpSampling2D, load_UpSampling2D, activate_UpSampling2D, init_UpSampling2D, update_UpSampling2D, input_shape, output_shape, size, activation_function)
        end
    end

    function init_UpSampling2D(layer::UpSampling2D, mini_batch::Int64)
        layer.value = zeros(Float32, layer.output_shape..., mini_batch)
        layer.output = zeros(Float32, layer.output_shape..., mini_batch)
        layer.∇biases = zeros(Float32, layer.output_shape..., mini_batch)
        layer.δ = zeros(Float32, layer.input_shape..., mini_batch)
    end

    function activate_UpSampling2D(layer::UpSampling2D, input::Array{Float32})
        x, y = layer.size
        @avxt for i in axes(input, 1), j in axes(input, 2), c in axes(input, 3), b in axes(input, 4)
            for s₁ in 1:layer.size[1], s₂ in 1:layer.size[2]
                layer.value[(i-1)*x+s₁, (j-1)*y+s₂, c, b] = input[i, j, c, b]
            end
        end

        layer.activation_function.func!(layer.output, layer.value)
    end

    function update_UpSampling2D(layer::UpSampling2D, optimizer::String, input::Array{Float32}, δₗ₊₁::Array{Float32}, α::Float64, parameters::Tuple, direction::Int64=0)
        layer.activation_function.get_∇biases!(layer.∇biases, layer.value, δₗ₊₁)
        layer.δ .= 0.0f0

        x, y = layer.size
        @avxt for i in axes(layer.δ, 1), j in axes(layer.δ, 2), c in axes(layer.δ, 3), b in axes(layer.δ, 4)
            for s₁ in 1:layer.size[1], s₂ in 1:layer.size[2]
                layer.δ[i, j, c, b] += layer.∇biases[(i-1)*x+s₁, (j-1)*y+s₂, c, b]
            end
        end
    end

    function save_UpSampling2D(layer::UpSampling2D, file::Any, id::String)
        write(file, id, "UpSampling2D")
        write(file, id*"input_shape", collect(layer.input_shape))
        write(file, id*"size", collect(layer.size))
        write(file, id*"activation_function", layer.activation_function.get_name())
    end

    function load_UpSampling2D(layer::UpSampling2D, file::Any, id::String)
    end

    function get_args(file::Any, id::String)
        input_shape = tuple(read(file, id*"input_shape")...)
        size = tuple(read(file, id*"size")...)
        return (input_shape=input_shape, size=size)
    end
end


text:
module TestRandomAdditions

using ConicHulls, ConicHulls.RefHull

function test(NC::Int, ngen::Int, r::Int)
    eval(ConicHulls.Hulls,:(numfacets=0))
    eval(ConicHulls.Primitives,:(numgenerators=0))
    test(create_hull(NC), ngen, r)
end

function test(hull::ConicHull, ngen::Int, r::Int)
    NC, F, G = nconic(hull), ftype(hull), gtype(hull)

    gs = init_hull!(hull, eye(Int, NC))

    for k=1:ngen
        verify(hull)
        verify_hull(hull, gs)

        x = rand(-r:r, NC-1)
        push!(x, 1-sum(x))
        g = G(x)

        add!(hull, g)
        push!(gs, g)
    end
end

srand(673826715664)
for NC in 3:5
    @time for r=1:8, ngen=1:8, k=1:(1<<(8-NC))/ngen; test(NC, ngen, r); end
end

end # module


text:
using MonteCarloMeasurements: Particles

@testset "from_namedtuple" begin
    rng = MersenneTwister(42)

    nchains, ndraws = 4, 10
    sizes = (x=(), y=(2,), z=(3, 4))
    dims = Dict("y" => ["yx"], "z" => ["zx", "zy"])
    coords = Dict("yx" => ["y1", "y2"], "zx" => 1:3, "zy" => 1:4)

    nts = [
        "NamedTuple" =>
            (; (k => randn(rng, nchains, ndraws, v...) for (k, v) in pairs(sizes))...),
        "Vector{NamedTuple}" => [
            (; (k => randn(rng, ndraws, v...) for (k, v) in pairs(sizes))...) for
            _ in 1:nchains
        ],
        "Matrix{NamedTuple}" => [
            (; (k => randn(rng, v...) for (k, v) in pairs(sizes))...) for
            _ in 1:nchains, _ in 1:ndraws
        ],
        "Vector{Vector{NamedTuple}}" => [
            [(; (k => randn(rng, v...) for (k, v) in pairs(sizes))...) for _ in 1:ndraws] for _ in 1:nchains
        ],
        "Vector{NamedTuple} particles" => [
            (; (k => Particles(randn(rng, ndraws, v...)) for (k, v) in pairs(sizes))...) for _ in 1:nchains
        ],
    ]

    @testset "posterior::$(type)" for (type, nt) in nts
        idata1 = from_namedtuple(nt; dims=dims, coords=coords, library="MyLib")
        test_namedtuple_data(
            idata1,
            :posterior,
            keys(sizes),
            nchains,
            ndraws;
            library="MyLib",
            coords=coords,
            dims=dims,
        )

        idata2 = convert_to_inference_data(nt; dims=dims, coords=coords, library="MyLib")
        test_namedtuple_data(
            idata2,
            :posterior,
            keys(sizes),
            nchains,
            ndraws;
            library="MyLib",
            coords=coords,
            dims=dims,
        )
    end

    @testset "$(group)" for group in [
        :posterior_predictive, :sample_stats, :predictions, :log_likelihood
    ]
        @testset "::$(type)" for (type, nt) in nts
            idata1 = convert_to_inference_data(
                nt; (group => nt,)..., dims=dims, coords=coords, library="MyLib"
            )
            test_namedtuple_data(
                idata1,
                :posterior,
                keys(sizes),
                nchains,
                ndraws;
                library="MyLib",
                coords=coords,
                dims=dims,
            )
            test_namedtuple_data(
                idata1,
                group,
                keys(sizes),
                nchains,
                ndraws;
                library="MyLib",
                coords=coords,
                dims=dims,
            )

            idata2 = convert_to_inference_data(
                nt; (group => keys(sizes),)..., dims=dims, coords=coords, library="MyLib"
            )
            test_namedtuple_data(
                idata2,
                group,
                keys(sizes),
                nchains,
                ndraws;
                library="MyLib",
                coords=coords,
                dims=dims,
            )
        end
    end

    @testset "$(group)" for group in [:prior_predictive, :sample_stats_prior]
        @testset "::$(type)" for (type, nt) in nts
            idata1 = from_namedtuple(;
                prior=nt, (group => nt,)..., dims=dims, coords=coords, library="MyLib"
            )
            test_namedtuple_data(
                idata1,
                :prior,
                keys(sizes),
                nchains,
                ndraws;
                library="MyLib",
                coords=coords,
                dims=dims,
            )
            test_namedtuple_data(
                idata1,
                group,
                keys(sizes),
                nchains,
                ndraws;
                library="MyLib",
                coords=coords,
                dims=dims,
            )

            idata2 = from_namedtuple(;
                prior=nt,
                (group => keys(sizes),)...,
                dims=dims,
                coords=coords,
                library="MyLib",
            )
            test_namedtuple_data(
                idata2,
                group,
                keys(sizes),
                nchains,
                ndraws;
                library="MyLib",
                coords=coords,
                dims=dims,
            )
        end
    end

    @testset "$(group)" for group in
                            [:observed_data, :constant_data, :predictions_constant_data]
        _, nt = nts[1]
        idata = from_namedtuple(
            nt;
            (group => Dict("w" => [1.0, 2.0]),)...,
            dims=Dict("w" => ["wx"]),
            coords=Dict("wx" => 1:2),
            library="MyLib",
        )
        @test idata isa InferenceData
        @test group in ArviZ.groupnames(idata)
        ds = getproperty(idata, group)
        sizes = dimsizes(ds)
        @test length(sizes) == 1
        @test "w" in keys(vardict(ds))
        @test dimdict(ds)["w"] == ("wx",)
        @test "inference_library" in keys(attributes(ds))
        @test attributes(ds)["inference_library"] == "MyLib"

        # ensure that dims are matched to named tuple keys
        # https://github.com/arviz-devs/ArviZ.jl/issues/96
        idata = from_namedtuple(
            nt;
            (group => (w=[1.0, 2.0],),)...,
            dims=Dict("w" => ["wx"]),
            coords=Dict("wx" => 1:2),
            library="MyLib",
        )
        @test idata isa InferenceData
        @test group in ArviZ.groupnames(idata)
        ds = getproperty(idata, group)
        sizes = dimsizes(ds)
        @test length(sizes) == 1
        @test "w" in keys(vardict(ds))
        @test dimdict(ds)["w"] == ("wx",)
        @test "inference_library" in keys(attributes(ds))
        @test attributes(ds)["inference_library"] == "MyLib"
    end
end


text:
using FastGaussQuadrature
using Test, LinearAlgebra, Random, SpecialFunctions

include("test_gausschebyshev.jl")
include("test_gausslegendre.jl")
include("test_gaussjacobi.jl")
include("test_gaussradau.jl")
include("test_gausslobatto.jl")
include("test_besselroots.jl")
include("test_gausslaguerre.jl")
include("test_gausshermite.jl")


text:
"""
    bench_core2core_latency([cpuids = 0:Sys.CPU_THREADS-1; nbench = 5, nsamples::Integer = 100, mode::Symbol = :min])
A tool for measuring core-to-core latency (i.e. inter-core latency) in nanoseconds.

The measured latencies correspond to a full roundtrip between two cores. Divide them by two to obtain an estimate for the time needed to fetch data from another core.

**Important:** At least two Julia threads are required (`julia -t2`)!

**Refs:** Largely inspired by [rigtorp/c2clat](https://github.com/rigtorp/c2clat) and [ajakubek/core-latency](https://github.com/ajakubek/core-latency).
"""
function bench_core2core_latency(cpuids=0:(Sys.CPU_THREADS - 1); nbench=5, kwargs...)
    # check validity of cpuids input
    for c in cpuids
        if c < 0 || c > Sys.CPU_THREADS
            @error("CPU IDs must all be non-negative and ≤ Sys.CPU_THREADS.")
        end
    end
    # backup current thread affinity
    pinning_before = getcpuids()
    # run benchmarks
    ncpuids = length(cpuids)
    latencies = zeros(ncpuids, ncpuids)
    for b in 1:nbench
        for (j, cpu2) in pairs(cpuids)
            for (i, cpu1) in pairs(cpuids)
                @inbounds latencies[i, j] += Core2CoreLatency.run_bench(
                    cpu1, cpu2; kwargs...
                )
            end
        end
    end
    latencies ./= nbench
    # restore previous thread affinity
    pinthreads(pinning_before)
    return latencies
end


text:
using PyPlot, StatsBase, Statistics

include("representable_numbers.jl")

float8 = representable_floats(8,3)
float16 = representable_floats(16,5)
bfloat16 = representable_floats(16,8)

f8_am, f8_wda = wcdp_float(float8)
f_am, f_wda = wcdp_float(float16)
bf_am, bf_wda = wcdp_float(bfloat16)

i_am = vcat(0.5000001,1.49999999:1:(2^15-1))
i_wda = decprec(i_am,round.(i_am))
i_am[1] = 1

## log fixed point
nint = 8
nfrac = 7

# all integer bit being 1 means NaN
lfxp_bf16 = 2 .^ collect((-2.0^(nint-1)):(2.0^-nfrac):(2.0^(nint-1)-1-2.0^-nfrac))

nint = 5
nfrac = 10
lfxp_f16 = 2 .^ collect((-2.0^(nint-1)):(2.0^-nfrac):(2.0^(nint-1)-2.0^-nfrac))

lfxp_bf16_am, lfxp_bf16_wda = wcdp_approx(lfxp_bf16)
lfxp_f16_am, lfxp_f16_wda = wcdp_approx(lfxp_f16)

# include floatmin/2 which gets round to floatmin
# include the no overflow rounding mode
lfxp_f16_am = vcat(lfxp_f16[1]/2,lfxp_f16_am,lfxp_f16[end]*10)
lfxp_f16_wda = vcat(decprec(lfxp_f16[1]/2,lfxp_f16[1]),lfxp_f16_wda,0)

## PLOTTING
ioff()
fig,ax1 = subplots(1,1,figsize=(7,3))

# dashed lines for reprsentable numbers
ax1.fill_between(i_am,-0.1,i_wda,edgecolor="C0",facecolor="none",linestyle="--")
ax1.plot(ones(2)*float16[1],[0,.5],"k",ls="--")
ax1.plot(ones(2)*float16[end],[0,4],"k",ls="--",lw=2)
ax1.plot(ones(2)*float8[1],[0,.5],"C2",ls="--")
ax1.plot(ones(2)*float8[end],[0,2],"C2",ls="--")
ax1.plot(ones(2)*lfxp_f16[1],[0,3.8],"C1",ls="--")
ax1.plot(ones(2)*lfxp_f16[end],[0,4],"C1",ls="--")

# floats and integers
ax1.plot(f_am,f_wda,"k",lw=2)
ax1.plot(bf_am,bf_wda,"grey",lw=1.4)
ax1.plot(f8_am,f8_wda,"C2",lw=1.3)
ax1.plot(i_am,i_wda,"C0",lw=2)

# log fix points
ax1.plot(lfxp_bf16_am,lfxp_bf16_wda,"C3",lw=2)
ax1.plot(lfxp_f16_am,lfxp_f16_wda,"C1",lw=2)

x0,x1 = 1e-16,1e16

ax1.set_xlim(x0,x1)
ax1.set_xscale("log",basex=10)
ax1.set_ylim(0,6)

ax1.set_xlabel("value")
ax1.set_ylabel("decimal places")

ax1.text(6e-3,1,"Float8",color="C2",rotation=64.5,va="bottom")
ax1.text(5e-7,4.1,"Float16",color="k")
ax1.text(1e-15,3.2,"BFloat16",color="grey")
ax1.text(6e4,5,"Int16",color="C0")

ax1.text(2e5,3.7,"LogFixPoint16",color="C1",rotation=0)
ax1.text(1e10,3.2,"BLogFixPoint16",color="C3",rotation=0)

ax1.set_title("Decimal precision",loc="left")

tight_layout()
savefig("/Users/milan/git/LogFixPoint16s.jl/figs/decimal_precision.png",dpi=100)
close(fig)


text:
a = 1
c = -5
limit = 10

mf = SigmoidMF(a, c, limit)

# Evaluation tests
@assert mf.eval(c) == 0.5
@assert mf.eval(c - 5) < mf.eval(c + 5)

# Mean finding tests
@assert mf.mean_at(1) > mf.mean_at(0.6) > mf.mean_at(0.3) > mf.mean_at(0)


text:
using SafeTestsets

@safetestset "is rectangular" begin
    using Test, TwoDots
    import TwoDots: is_rectangular
    game = Game()
    map(x->x.color = "green", game.dots)
    push!(game.selected_dots, game.dots[1,1])
    push!(game.selected_dots, game.dots[2,1])
    push!(game.selected_dots, game.dots[3,1])
    map(x->x.selected = true, game.selected_dots)
    @test !is_rectangular(game)

    game = Game()
    map(x->x.color = "green", game.dots)
    push!(game.selected_dots, game.dots[1,1])
    map(x->x.selected = true, game.selected_dots)
    @test !is_rectangular(game)

    game = Game()
    map(x->x.color = "green", game.dots)
    push!(game.selected_dots, game.dots[1,1])
    push!(game.selected_dots, game.dots[1,2])
    push!(game.selected_dots, game.dots[1,3])
    map(x->x.selected = true, game.selected_dots)
    @test !is_rectangular(game)

    game = Game()
    map(x->x.color = "green", game.dots)
    push!(game.selected_dots, game.dots[1,1])
    push!(game.selected_dots, game.dots[1,2])
    push!(game.selected_dots, game.dots[1,3])
    push!(game.selected_dots, game.dots[2,3])
    map(x->x.selected = true, game.selected_dots)
    @test !is_rectangular(game)

    game = Game()
    map(x->x.color = "green", game.dots)
    push!(game.selected_dots, game.dots[1,1])
    push!(game.selected_dots, game.dots[1,2])
    push!(game.selected_dots, game.dots[1,3])
    push!(game.selected_dots, game.dots[2,3])
    push!(game.selected_dots, game.dots[2,2])
    map(x->x.selected = true, game.selected_dots)
    @test is_rectangular(game)

    game = Game()
    map(x->x.color = "green", game.dots)
    push!(game.selected_dots, game.dots[1,1])
    push!(game.selected_dots, game.dots[1,2])
    push!(game.selected_dots, game.dots[1,3])
    push!(game.selected_dots, game.dots[2,3])
    push!(game.selected_dots, game.dots[3,3])
    push!(game.selected_dots, game.dots[3,2])
    push!(game.selected_dots, game.dots[3,1])
    push!(game.selected_dots, game.dots[2,1])
    map(x->x.selected = true, game.selected_dots)
    @test is_rectangular(game)
end

@safetestset "click submit" begin
    using Test, TwoDots
    import TwoDots: click_submit!

    game = Game(n_rows=6, n_cols=2)
    game.dots[1,1].color = "green"
    game.dots[2,1].color = "red"
    game.dots[3,1].color = "yellow"
    game.dots[4,1].color = "green"
    game.dots[5,1].color = "purple"
    game.dots[6,1].color = "orange"
    map(x->x.color = "green", game.dots[:,2])

    push!(game.selected_dots, game.dots[3,1])
    push!(game.selected_dots, game.dots[5,1])

    click_submit!(game)

    @test isempty(game.selected_dots)

    @test game.dots[3,1].color == "green"
    @test game.dots[4,1].color == "red"
    @test game.dots[5,1].color == "green"
    @test game.dots[6,1].color == "orange"

    @test all(x->x.color == "green", game.dots[:,2])

    game = Game(n_rows=6, n_cols=2)
    game.dots[1,1].color = "green"
    game.dots[2,1].color = "red"
    game.dots[3,1].color = "yellow"
    game.dots[4,1].color = "green"
    game.dots[5,1].color = "purple"
    game.dots[6,1].color = "orange"
    map(x->x.color = "green", game.dots[:,2])

    # opposite order
    push!(game.selected_dots, game.dots[5,1])
    push!(game.selected_dots, game.dots[3,1])

    click_submit!(game)

    @test isempty(game.selected_dots)

    @test game.dots[3,1].color == "green"
    @test game.dots[4,1].color == "red"
    @test game.dots[5,1].color == "green"
    @test game.dots[6,1].color == "orange"

    @test all(x->x.color == "green", game.dots[:,2])
end


text:
abstract type HMCIntegrator end

function _ahmc_step_size(rng::AbstractRNG, integrator::HMCIntegrator, hamiltonian::AdvancedHMC.Hamiltonian, x_init::AbstractVector{<:Real})
    if isnan(integrator.step_size)
        oftype(integrator.step_size, AdvancedHMC.find_good_stepsize(rng, hamiltonian, x_init, max_n_iters = 100))
    else
        integrator.step_size
    end
end

@with_kw struct LeapfrogIntegrator <: HMCIntegrator
    step_size::Float64 = NaN
end

function ahmc_integrator(rng::AbstractRNG, integrator::LeapfrogIntegrator, hamiltonian::AdvancedHMC.Hamiltonian, x_init::AbstractVector{<:Real})
    return AdvancedHMC.Leapfrog(_ahmc_step_size(rng, integrator, hamiltonian, x_init))
end

@with_kw struct JitteredLeapfrogIntegrator <: HMCIntegrator
    step_size::Float64 = NaN
    jitter_rate::Float64 = 1.0
end

function ahmc_integrator(rng::AbstractRNG, integrator::JitteredLeapfrogIntegrator, hamiltonian::AdvancedHMC.Hamiltonian, x_init::AbstractVector{<:Real})
    return AdvancedHMC.JitteredLeapfrog(_ahmc_step_size(rng, integrator, hamiltonian, x_init), integrator.jitter_rate)
end

@with_kw struct TemperedLeapfrogIntegrator <: HMCIntegrator
    step_size::Float64 = NaN
    tempering_rate::Float64 = 1.05
end

function ahmc_integrator(rng::AbstractRNG, integrator::TemperedLeapfrogIntegrator, hamiltonian::AdvancedHMC.Hamiltonian, x_init::AbstractVector{<:Real})
    return AdvancedHMC.TemperedLeapfrog(_ahmc_step_size(rng, integrator, hamiltonian, x_init), integrator.tempering_rate)
end


text:
function logprior(state::StateAM, prior::PriorAM)
  lp = logpdf(prior.gamma, state.gammaC)
  lp += logpdf(prior.gamma, state.gammaT)
  lp += logpdf(prior.eta, state.etaC)
  lp += logpdf(prior.eta, state.etaT)
  lp += sum(logpdf.(prior.mu, state.mu))
  lp += sum(logpdf.(prior.sigma, state.sigma))
  lp += sum(logpdf.(prior.phi, state.phi))
  lp += sum(logpdf.(prior.nu, state.nu))
  return lp
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

function _euler(ai::T, bi::T, select::Integer, FK4::Bool, radians::Bool) where {T<:AbstractFloat}

    if select>6 || select<1
        error("Input for coordinate transformation should be an integer in the range 1:6")
    end

    if FK4
        psi = (0.57595865315, 4.9261918136, 0.0, 0.0, 0.11129056012, 4.7005372834)
        stheta = (0.88781538514, -0.88781538514,0.39788119938, -0.39788119938,
                  0.86766174755, -0.86766174755)
        ctheta = (0.46019978478, 0.46019978478, 0.91743694670, 0.9174369467,
                  0.49715499774, 0.49715499774)
        phi = (4.9261918136,  0.57595865315, 0.0, 0.0, 4.7005372834, 0.11129056012)
    else
        psi = (0.574770433, 4.9368292465, 0.0, 0.0, 0.11142137093, 4.71279419371)
        stheta = (0.88998808748, -0.88998808748, 0.39777715593, -0.39777715593,
                  0.86766622025, -0.86766622025)
        ctheta = (0.45598377618, 0.45598377618, 0.91748206207, 0.91748206207,
                  0.49714719172, 0.49714719172)
        phi = (4.9368292465, 0.574770433, 0.0, 0.0, 4.71279419371, 0.11142137093)
    end

    if !radians
        ai = deg2rad(ai)
        bi = deg2rad(bi)
    end
    sa, ca = sincos(ai - phi[select])
    sb, cb = sincos(bi)
    x = (cb * sa, cb * ca, sb)
    bo = ctheta[select]*x[3] - stheta[select]*x[1]
    ao = mod2pi(atan(ctheta[select]*x[1] + stheta[select]*x[3], x[2]) + psi[select])
    bo = asin(bo)

    if radians
        return (ao, bo)
    end
    return (rad2deg(ao), rad2deg(bo))
end

"""
    euler(ai, bi, select[, FK4=true, radians=true])

### Purpose ###

Transform between Galactic, celestial, and ecliptic coordinates.

### Explanation ###

The function is used by the astro procedure.

### Arguments ###

* `ai`: input longitude, scalar or vector.
* `bi`: input latitude, scalar or vector.
* `select` : integer input specifying type of coordinate
  transformation.
  SELECT   From          To     | SELECT   From       To
     1   RA-Dec (2000) Galactic |   4    Ecliptic   RA-Dec
     2   Galactic      RA-DEC   |   5    Ecliptic   Galactic
     3   RA-Dec        Ecliptic |   6    Galactic   Ecliptic

* `FK4` (optional boolean keyword) : if this keyword is set to `true`,
  then input and output celestial and ecliptic coordinates should be
  given in equinox B1950. When `false`, by default, they should be given in
  equinox J2000.
* `radians` (optional boolean keyword) : if this keyword is set to
  `true`, all input and output angles are in radians rather than degrees.

### Output ###

a 2-tuple `(ao, bo)`:

* `ao`: output longitude in degrees.
* `bo`: output latitude in degrees.

### Example ###

Find the Galactic coordinates of Cyg X-1 (ra=299.590315, dec=35.201604)

```jldoctest
julia> using AstroLib

julia> euler(299.590315, 35.201604, 1)
(71.33498957116959, 3.0668335310640984)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
euler(ai::Real, bi::Real, select::Integer; FK4::Bool=false, radians::Bool=false) =
    _euler(promote(float(ai), float(bi))..., select, FK4, radians)

euler(aibi::Tuple{Real, Real}, select::Integer; FK4::Bool=false, radians::Bool=false) =
    euler(aibi[1], aibi[2], select, FK4=FK4, radians=radians)

function euler(ai::AbstractVector{R}, bi::AbstractVector{<:Real}, select::Integer;
               FK4::Bool=false, radians::Bool=false) where {R<:Real}
    @assert length(ai) == length(bi) "ai and bi arrays should be of the same length"
    typeai = float(R)
    ai_out  = similar(ai,  typeai)
    bi_out = similar(bi, typeai)
    for i in eachindex(ai)
        ai_out[i], bi_out[i] = euler(ai[i], bi[i], select,
                                        FK4=FK4, radians=radians)
    end
    return ai_out, bi_out
end


text:
include("test_lattice.jl")
include("test_model.jl")
include("test_hamiltonian.jl")
include("test_mesh.jl")
# include("test_bandstructure.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ e7af1703-0873-4b3f-8b8f-a8a2c874bcb1
begin
	using PlutoUI, LinearAlgebra, Distributions, SparseArrays
	import Random
end

# ╔═╡ cf46a0fe-43a4-4b6a-b4f8-a0c5fb4e9f03
# For binder, uncomment this cell ...
#=
begin
	import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
		Pkg.PackageSpec(name="PlutoUI"),
        Pkg.PackageSpec(name="Distributions"),
		Pkg.PackageSpec(name="SparseArrays")
    ])
end
=#

# ╔═╡ cfaedd73-8e4c-4481-99d0-c04eafd79fe9
TableOfContents(title="📚 Table of Contents",aside=true)

# ╔═╡ 00839532-0da0-4cbf-8093-01310a3799d0
md"""
# Randomized Linear Algebra

See [Yuxin Chen, Randomized linear algebra](https://www.princeton.edu/~yc5/orf570/randomized_linear_algebra.pdf).

"""

# ╔═╡ a6f4632f-5d64-4efb-a72f-706235cb3957
md"""
# Matrix multiplication

Let $A\in\mathbb{R}^{m\times n}$ and $B\in\mathbb{R}^{n\times p}$. Then,

$$
C=AB=\sum_{i=1}^{n} A_{:,i} B_{i,:}.$$

Assume, for simplicity $m=n=p$.

__Idea:__ approximate $C$ by randomly sampling $r$ rank-one components.

__Algorithm:__ for $l= 1,\cdots ,r$ pick $i_l\in\{1,···,n\}$ i.i.d. with probability $\mathbb{P}\{i_l=k\}=p_k$ and compute

$$M=\sum_{l=1}^r \frac{1}{rp_{i_l}} A_{:,i_l} B_{i_l,:}$$

__Rationale:__ $M$ is an unbiased estimate of $C$,

$$\begin{aligned}
\mathbb{E}[M]&=\sum_{l=1}^r \sum_k \mathbb{P}\{i_l=k\} \frac{1}{r p_k} A_{:,k}B_{k,:}
=\sum_k A_{:,k}B_{k,:}=C.
\end{aligned}$$

__Importance sampling porbabilities__ $p_k$ are

* __Uniform sampling:__ $p_k=\displaystyle\frac{1}{n}$

* __Nonuniform sampling:__

$$
p_k=\frac{\|A_{:,k}\|_2 \|B_{k,:}\|_2}
{\sum_l \|A_{:,l}\|_2 \|B_{l,:}\|_2},  \tag{1}$$

which is computable in one-pass and requires $O(n)$ memory and $O(n^2)$ operations.

__Theorem.__ [Optimality] $\mathbb{E}[\|M-AB\|_F^2]$ is minimized for $p_k$ given by (1).

__Theorem.__ [Error] Choose $p_k\geq \displaystyle\frac{\beta \|A_{:,k}\|_2 \|B_{k,:}\|_2} {\sum_l \|A_{:,l}\|_2 \|B_{l,:}\|_2}$ for some
$0<\beta \leq 1$. If $r\geq \displaystyle\frac{\log n}{\beta}$, then

$$\|M-AB\|_F\leq \sqrt{\frac{\log n}{\beta r}}
\|A\|_F \|B\|_F$$

with probability exceeding $1-O(n^{-10})$.
"""

# ╔═╡ 25a893eb-2c23-40dc-8e24-f45ed8253697
md"
## Full matrix
"

# ╔═╡ 4dc397d7-5e07-4da9-94ea-fa22145d7222
begin
	Random.seed!(6789)
	n=3000
	A=rand(n,n)
	B=rand(n,n)
	β=1.0
	log(n)/β
end

# ╔═╡ 31305f95-1a9d-4981-bdff-bd8f8edc3802
begin
	# Uniform
	r=1000
	iᵣ=rand(1:n,r)
	p=1/n
	@time M=A[:,iᵣ]*B[iᵣ,:]/(r*p)
	@time C=A*B;
end

# ╔═╡ 8422981a-2daf-4b93-b1a9-499f1f1f7eaa
norm(M-C), norm(C)

# ╔═╡ 477d5ccd-b9e8-482e-ad73-ae1c56eb4cfc
begin
	# Nonuniform
	pA=[norm(view(A,:,k)) for k=1:n]
	pB=[norm(view(B,k,:)) for k=1:n]
	s=pA⋅pB
	p₁=[pA[k]*pB[k]/s for k=1:n]
	sum(p₁)
end

# ╔═╡ f0733a4f-5a05-41b4-9313-678ecad74356
begin
	i₁=rand(Categorical(p₁),r)
	@time Mₙ=A[:,i₁]*inv(Diagonal(r*p₁[i₁]))*B[i₁,:]
	norm(Mₙ-C), norm(C), √(log(n)/(β*r))*norm(A)*norm(B)
end

# ╔═╡ b0908c15-d3f1-429e-8165-e40bf58fd991
md"
## Sparse matrix
"

# ╔═╡ 2ec3b0a7-d24f-4fb5-864d-65b7f1ec0e06
begin
	# Sparse, nonuniform
	Random.seed!(9751)
	nₛ=10000
	Aₛ=sprand(nₛ,nₛ,0.1)
	Bₛ=sprand(nₛ,nₛ,0.1)
	@time Cₛ=Aₛ*Bₛ
	βₛ=1.0
	log(nₛ)/βₛ
end

# ╔═╡ 1260c0a9-2e6c-4c46-b84d-533c19cc2f67
# Cₛ is full
nnz(Cₛ)/prod(size(Cₛ))

# ╔═╡ c906c80d-00e0-41da-8c6e-fa0d58c4ae60
begin
	# Nonuniform
	pAₛ=[norm(view(Aₛ,:,k)) for k=1:nₛ]
	pBₛ=[norm(view(Bₛ,k,:)) for k=1:nₛ]
	t=pAₛ⋅pBₛ
	pₛ=[pAₛ[k]*pBₛ[k]/t for k=1:nₛ]
	sum(pₛ)
end

# ╔═╡ 402e6d6d-502d-41bf-acff-1d9c09bd0bf0
begin
	rₛ=3000
	iₛ=rand(Categorical(pₛ),rₛ);
	@time Mₛ=Aₛ[:,iₛ]*inv(Diagonal(rₛ*pₛ[iₛ]))*Bₛ[iₛ,:]
	norm(Mₛ-Cₛ),√(log(nₛ)/(βₛ*rₛ))*norm(Aₛ)*norm(Bₛ), norm(Cₛ)
end

# ╔═╡ 3800b4ce-10fb-4300-9129-09a0ca744a89

# ╔═╡ Cell order:
# ╠═cf46a0fe-43a4-4b6a-b4f8-a0c5fb4e9f03
# ╠═e7af1703-0873-4b3f-8b8f-a8a2c874bcb1
# ╠═cfaedd73-8e4c-4481-99d0-c04eafd79fe9
# ╟─00839532-0da0-4cbf-8093-01310a3799d0
# ╟─a6f4632f-5d64-4efb-a72f-706235cb3957
# ╟─25a893eb-2c23-40dc-8e24-f45ed8253697
# ╠═4dc397d7-5e07-4da9-94ea-fa22145d7222
# ╠═31305f95-1a9d-4981-bdff-bd8f8edc3802
# ╠═8422981a-2daf-4b93-b1a9-499f1f1f7eaa
# ╠═477d5ccd-b9e8-482e-ad73-ae1c56eb4cfc
# ╠═f0733a4f-5a05-41b4-9313-678ecad74356
# ╟─b0908c15-d3f1-429e-8165-e40bf58fd991
# ╠═2ec3b0a7-d24f-4fb5-864d-65b7f1ec0e06
# ╠═1260c0a9-2e6c-4c46-b84d-533c19cc2f67
# ╠═c906c80d-00e0-41da-8c6e-fa0d58c4ae60
# ╠═402e6d6d-502d-41bf-acff-1d9c09bd0bf0
# ╠═3800b4ce-10fb-4300-9129-09a0ca744a89

max_stars_repo_path: src/Module C - Applications/L17 Randomized Linear Algebra.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
# This code shows the calculation of rotor losses in high speed PM machines
# author: Jaime Renedo Anglada (renedo.jaime@ieee.org)

# Define geometry:

mu_0=4*pi*10.0^-7; # [m kg s^-2 A^-2]
p=2; # number of pole pairs

B_r_ext=2.553*10.0^-3; # external field [T]
space_order=9; # space order of the harmonic
q=p*space_order; # harmonic order

n_mech=65000; # mechanical speed [rpm]
k_time=1*6; # time order
f_1=n_mech/60*p; # [Hz]

N_div=1000.0; # number of divisions simplified method

mu_r=1.07;

# Conductivity:
sigma_1=6.7*10.0^6; # [S/m]
sigma_2=0.77*10.0^6; # [S/m]
sigma_3=2.2*10.0^4; # [S/m]
sigma_4=3*10.0^-15; # [S/m]
sigma_5=3*10.0^-15; # [S/m]

# test:
sigma_3=3*10.0^-15; # [S/m]
#sigma_3=8*10^5
#sigma_sleeve=3*10.0^-15; # [S/m]

# Permeability:
mu_1=750*mu_0; # [m kg s^-2 A^-2]
mu_2=mu_r*mu_0; # [m kg s^-2 A^-2]
mu_3=1*mu_0; # [m kg s^-2 A^-2]
mu_4=1*mu_0; # [m kg s^-2 A^-2]
mu_5=5000*mu_0; # [m kg s^-2 A^-2]

# Geometry:
t_sl=2.95*10.0^-3; # [m] sleeve thickness

R_1=20*10.0^-3; # [m] rotor radius
R_2=27.4*10.0^-3; # [m] PM radius
R_3=R_2+t_sl; # [m] sleeve outer radius
R_4=31.15*10.0^-3; # [m] current sheet radius
R_5=50*10.0^-3; # [m] end domain radius

r_wave=R_2+0.25*10.0^-3; # where you measured the harmonics

#B_r_ext=B_r_ext*(R_3/(R_2+0.5*10.0^-3))^(q-1);

L=109*10.0^-3; # [m] axial length

println("R_1 = $R_1")
println("R_2 = $R_2")
println("R_3 = $R_3")
println("R_4 = $R_4")
println("R_5 = $R_5")
println("r_wave = $r_wave")

println("B_r_ext = $(B_r_ext*1000)")

# no load losses

using SpecialFunctions

vec_Br=[5.4033 8.703 2.4408 1.407 2.1913 3.5324 0.9911]*10.0^-3;
vec_so=[5 7 9 11 11 13 15];
vec_time=[1 1 1 1 2 2 2];

vec_loss=zeros(1,length(vec_Br));

for count=1:length(vec_Br)
    space_order=vec_so[count]
    B_r_ext=vec_Br[count]
    k_time=vec_time[count]

    vec_loss[count]=Calc_rotor_loss( n_mech, space_order, p, k_time, B_r_ext,L,R_1,R_2,R_3,R_4,R_5,sigma_1,sigma_2,sigma_3,sigma_4,sigma_5,mu_r,r_wave);

end

count=1;

# Print no load losses:
println("$vec_loss")

# on load losses

using SpecialFunctions

vec_Br=[5.4033+7.8406 8.703+4.86 2.553 1.407 2.1913+1.827 3.5324+1 0.9911]*10.0^-3;
vec_so=[5 7 9 11 11 13 15];
vec_time=[1 1 1 1 2 2 2];

vec_loss=zeros(1,length(vec_Br));

for count=1:length(vec_Br)
    space_order=vec_so[count]
    B_r_ext=vec_Br[count]
    k_time=vec_time[count]

    vec_loss[count]=Calc_rotor_loss( n_mech, space_order, p, k_time, B_r_ext,L,R_1,R_2,R_3,R_4,R_5,sigma_1,sigma_2,sigma_3,sigma_4,sigma_5,mu_r,r_wave);

end

count=1;

# Print on load losses:
println("$vec_loss")


text:
module Cond
using MLStyle

export @cond
@nospecialize
function cond(cases, source, mod::Module)
    @switch cases begin
        @case Expr(:block, cases...)
        default = Expr(
            :call,
            throw,
            "None of the branches have satisfied conditions, at $(string(source)).",
        )
        last_lnode = source
        folded = foldr(cases, init = default) do case, last
            @switch case begin
                @case ::LineNumberNode
                last_lnode = case
                return last
                @case :(_ => $b)
                return b
                @case :($a => $b)
                return Expr(:if, a, b, Expr(:block, last_lnode, last))
                @case _
                throw("Invalid syntax for conditional branches at $last_lnode.")
            end
        end
        return Expr(:block, source, folded)
        @case _
        msg = "Malformed ast template, the second arg should be a block with a series of pairs(`a => b`), at $(string(source))."
        throw(SyntaxError(msg))
    end
end

"""
```julia
@cond begin
    x > 1 => true
    _ => false
end
```
"""
macro cond(cases)
    cond(cases, __source__, __module__) |> esc
end
@specialize
end


text:
# # [Linear Time Discrete System](@id linear_discrete)
#
# We will start by estimating the underlying dynamical system of a time discrete process based on some measurements via [Dynamic Mode Decomposition](https://arxiv.org/abs/1312.0041) on a simple linear system of the form ``u(k+1) = A u(k)``.
#
# At first, we simulate the correspoding system using `OrdinaryDiffEq.jl` and generate a [`DiscreteDataDrivenProblem`](@ref DataDrivenProblem) from the simulated data.

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
#md using Plots

A = [0.9 -0.2; 0.0 0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 11.0)

f(u,p,t) = A*u

sys = DiscreteProblem(f, u0, tspan)
sol = solve(sys, FunctionMap());

# Next we transform our simulated solution into a [`DataDrivenProblem`](@ref). Given that the solution knows its a discrete solution, we can simply write

prob = DataDrivenProblem(sol)

# And plot the solution and the problem

#md plot(sol, label = string.([:x₁ :x₂]))
#md scatter!(prob)

# To estimate the underlying operator in the states ``x_1, x_2``, we `solve` the estimation problem using the [`DMDSVD`](@ref) algorithm for approximating the operator. First, we will have a look at the [`DataDrivenSolution`](@ref)

res = solve(prob, DMDSVD(), digits = 1)
#md println(res) # hide

# We see that the system has been recovered correctly, indicated by the small error and high AIC score of the result. We can confirm this by looking at the resulting [`Basis`](@ref)

system = result(res)
using Symbolics

#md println(system) # hide

# And also plot the prediction of the recovered dynamics

#md plot(res)

# Or a have a look at the metrics of the result

#md metrics(res)

# To have a look at the representation of the operator as a `Matrix`, we can simply call

#md Matrix(system)

# to see that the operator is indeed our initial `A`. Since we have a linear representation, we can gain further insights into the stability of the dynamics via its eigenvalues

#md eigvals(system)

# And plot the stability margin of the discrete System

#md φ = 0:0.01π:2π
#md plot(sin.(φ), cos.(φ), xlabel = "Real", ylabel = "Im", label = "Stability margin", color = :red, linestyle = :dash)
#md scatter!(real(eigvals(system)), imag(eigvals(system)), label = "Eigenvalues", color = :black, marker = :cross)

# Similarly, we could use a sparse regression to derive our system from our data. We start by defining a [`Basis`](@ref)

using ModelingToolkit

@parameters t
@variables x[1:2](t)

basis = Basis(x, x, independent_variable = t, name = :LinearBasis)
#md print(basis) #hide

# Afterwards, we simply `solve` the already defined problem with our `Basis` and a `SparseOptimizer`

sparse_res = solve(prob, basis, STLSQ())
#md println(sparse_res) #hide

# Which holds the same equations
sparse_system = result(sparse_res)
#md println(sparse_system) #hide

# Again, we can have a look at the result

#md plot(
#md     plot(prob), plot(sparse_res), layout = (1,2)
#md )

# Both results can be converted into a `DiscreteProblem`

@named sys = DiscreteSystem(equations(sparse_system), get_iv(sparse_system),states(sparse_system), parameters(sparse_system))
#md println(sys) #hide

# And simulated using `OrdinaryDiffEq.jl` using the (known) initial conditions and the parameter mapping of the estimation.

x0 = [x[1] => u0[1], x[2] => u0[2]]
ps = parameter_map(sparse_res)

discrete_prob = DiscreteProblem(sys, x0, tspan, ps)
estimate = solve(discrete_prob, FunctionMap());

# And look at the result
#md plot(sol, color = :black)
#md plot!(estimate, color = :red, linestyle = :dash)

#md # ## [Copy-Pasteable Code](@id linear_discrete_copy_paste)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

## Test the result #src
for r_ in [res, sparse_res] #src
    @test all(l2error(r_) .<= 1e-10) #src
    @test all(aic(r_) .>= 1e10) #src
    @test all(determination(r_) .≈ 1.0) #src
end  #src
@test Array(sol) ≈ Array(estimate) #src


text:
#=
gaussian_smoothing:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-04-08
=#
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays
# using Pandas
using Pio3d
# using Seaborn

V = [[0 0]; [2 0]; [1 2]; [3 2]; [4 2]; [4 0] ]'
EVch = spzeros(Int8, 7, 6)
EVch[1,1] = 1
EVch[1,2] = 1
EVch[2,2] = 1
EVch[2,3] = 1
EVch[3,1] = 1
EVch[3,3] = 1
EVch[4,3] = 1
EVch[4,4] = 1
EVch[5,2] = 1
EVch[5,4] = 1
EVch[6,2] = 1
EVch[6,6] = 1
EVch[7,4] = 1
EVch[7,6] = 1
Matrix(EVch)

FVch = spzeros(Int8, 3, 6)
FVch[1, 1] = 1
FVch[1, 2] = 1
FVch[1, 3] = 1
FVch[2, 2] = 1
FVch[2, 3] = 1
FVch[2, 4] = 1
# FVch[3, 2] = 1
# FVch[3, 4] = 1
# FVch[3, 6] = 1

FEch = spzeros(Int8, 3, 7)
FEch[1, 1] = 1
FEch[1, 2] = 1
FEch[1, 3] = 1
FEch[2, 2] = 1
FEch[2, 4] = 1
FEch[2, 5] = 1

function eye(sz)
    arr = spzeros(Int8, sz, sz)
    for i=1:sz
        arr[i,i] = 1
    end
    return arr
end

function setIzero!(arr)
    for i=1:size(arr)[1]
        arr[i,i] = 0
    end
    return arr
end

function getDiag(arr)
    sz = size(arr)[1]
    arro = zeros(eltype(arr), sz)
    for i=1:sz
        arro[i] = arr[i,i]
    end
    return arro
end

function make_full(arr)
    arrout = Array{eltype(arr), 2}(undef, size(arr)...)
    for i=1:size(arr)[1]
        for j=1:size(arr)[2]
            arrout[i,j] = arr[i,j]
        end
    end
    return arrout
end
#neighboor vertexes

function smoothing_EV(V, EVch::SparseMatrixCSC, k=0.35)
#     Matrix(setIzero!(EVch' * EVch))

    neighboorNumber = getDiag(EVch' * EVch)
    neighboors = setIzero!(EVch' * EVch)

#     targetV =  neighboors * V' ./ neighboorNumber
#     diff = targetV - V'
#     newV = (V' + k * diff)'

    targetV = V * neighboors ./ neighboorNumber'
    println("targetV shape: ", size(targetV))
    diff = targetV - V
    println("diff shape: ", size(diff))
    newV = V + k * diff

    return make_full(newV)
end

smoothing_EV(V, EVch)

## Read data from file
xystep = 1
zstep = 1
threshold = 4000;
pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")
datap = Pio3d.read3d(pth)

data3d = datap["data3d"]
segmentation = data3d .> threshold

block_size = [5,5,5]
basicmodel, Flin, larmodel = LarSurf.get_surface_grid_per_block(segmentation, block_size; return_all=true)
someV, topology = basicmodel
filtered_bigFV = topology[1]
# someV, filtered_bigFV = basicmodel
bigV = someV
# bigV, tmodel = larmodel
# bigVV, bigEV, bigFV, bigCV = tmodel
# return (bigV,[FVreduced])

Plasm.View(basicmodel)
# Plasm.View((bigV,[bigVV, bigEV, filtered_bigFV]))

FV = filtered_bigFV
function get_EV(FV)
	# @info "getFV " FV
	EV = []
	for f in FV
		push!(EV, [[f[1],f[2]],[f[3],f[4]],  [f[1],f[3]],[f[2],f[4]]])
	end
	# display(cat(EV))
	@info "EV example !" EV EV[1]

	catev = cat(EV)
	@info "cat(EV)" catev
	@info "hcat(EV)" hcat(EV)
	@info "vcat(EV)" vcat(EV)
	doubleedges = Base.sort(catev)

	# doubleedges = Base.sort(cat(EV;dims=2))
	@info "EV double example" doubleedges
	doubleedges = convert(Lar.Cells, doubleedges)
	EV = [doubleedges[k] for k=1:2:length(doubleedges)]
	return EV
end
function smoothing_FV(V::Array, FV::Array{Array{Int64,1},1})
	# @info "smoothing FV " FV
	EV = get_EV(FV)
	aEV = LarSurf.ll2array(EV)

	# nrn4

	# notAllBigEV =

	kEV = LarSurf.characteristicMatrix(aEV, size(bigV)[2])
	# kEV = Lar.characteristicMatrix(EV)
	newBigV = smoothing_EV(V, kEV, 0.6)
	return newBigV
end

@info "smoothing_FV example" bigV FV FV[1]
newBigV = smoothing_FV(bigV, FV)

Plasm.View((newBigV * 100,[filtered_bigFV]))
# Plasm.View((newBigV * 100,[bigVV, bigEV, filtered_bigFV]))


text:
# # Battle Royale
# ![](battle.gif)
#
# This example illustrates how to leverage higher dimensions of a `GridSpace` to identify
# the distance from neighbors not just spatially, but also categorically. We'll also use
# the [`walk!`](@ref) function extensively.

# The `Models` module includes this example as [`Models.battle`](@ref).

# ## Rules of Engagement
# Agents wander around the map looking for opponents. When a grid space is occupied by
# two or more agents there will be a battle. With experience gained from the fight, the
# victor searches for more opponents to crush and losers scurry away defeated or possibly
# even die. This process repeats until there is a single, definitive winner.
#
# For this battle ground to exist, the following rules must be followed:
# - Agents have an experience level, starting at level 1 up to a maximum of 10.
# - Agents will search for the nearest worthy opponent (one with equal or ±1 experience level) and move towards them to attack, so long as something more important doesn't happen, which could be
#   - A tougher opponent (with experience level +2 or higher) is nearby: run!
#   - There are no worthy opponents available, but there are weak ones (with experience level -2 or lower): chase them down.
#   - Capture and taunt a weaker opponent, then kill them.
#   - Notice a tough opponent is occupied, sneak up and kill them.
#   - There is no-one worthy to fight, but also no-one left to taunt. All bets are off: THERE CAN BE ONLY ONE.

# Battles are won by weighted chance - a higher level gives an agent a larger chance of
# winning, but does not guarantee it. When a victor is chosen
#   - The difference in experience between opponents is swapped.
#   - If an agents experience reaches 0, they die.
#
# Captured opponents will be killed once taunted. The captor will gain half of their
# experience. If an opportunist manages to take the captor by surprise, they can gain
# up to half of the captor's experience. This means a level 1 agent may eliminate a
# level 10 captor and jump straight to level 6.
#
# Once all rules of engagement have been exhausted, the final showdown begins. Opponents
# fight their closest adversary regardless of experience level. Winner takes all.

# ## Model Setup

cd(@__DIR__) #src
using Random # hide
using Agents
using Plots
gr() # hide

mutable struct Fighter <: AbstractAgent
    id::Int
    pos::Dims{3}
    has_prisoner::Bool
    capture_time::Int
    shape::Symbol # For plotting
end

# As you can see, the properties of out agent are very simple and contain only two
# parameters that are needed to store context from one time step to the next. All
# other properties needed are stored in the space. `pos` is three-dimensional, two
# for the actual space agents move within, and a third categorical dimension representing
# their level. `shape` is used solely for plotting (well, used once just for convenience).

# Now let's set up the battle field:
function battle(; fighters = 50)
    model = ABM(
        Fighter,
        GridSpace((100, 100, 10); periodic = false);
        scheduler = random_activation,
    )

    n = 0
    while n != fighters
        pos = (rand(1:100, 2)..., 1) # Start at level 1
        if isempty(pos, model)
            add_agent!(pos, model, false, 0, :diamond)
            n += 1
        end
    end

    return model
end

Random.seed!(6547) # hide
model = battle()

# 50 opponents positioned randomly on a 100x100 grid, with no escape
# (`periodic = false`). To leverage categorical dimensions fully, non-periodic chebyshev
# space is necessary.

# ## Game Dynamics
#
# To implement the rules of engagement, only an `agent_step!` function is required,
# along with a few helper functions.

space(agent) = agent.pos[1:2]
level(agent) = agent.pos[3]

# `space` allows us to invoke a number of helpful utilities provided by Agents.jl
# but only operate on our spatial dimensions, `level` is a wrapper to access the agent's
# experience easily.

# Nearest agents that satisfy our search criteria can be identified via Euclidean
# distance solely on the spatial dimensions of our `GridSpace`.

function closest_target(agent::Fighter, ids::Vector{Int}, model::ABM)
    if length(ids) == 1
        closest = ids[1]
    else
        close_id = argmin(map(id -> edistance(space(agent), space(model[id]), model), ids))
        closest = ids[close_id]
    end
    return model[closest]
end

# Since our battles are only between opponents with equal, or as much as one level apart,
# the odds can be set explicitly. Stronger opponents have twice the capacity of winning a
# match.

function battle!(one::Fighter, two::Fighter, model)
    if level(one) == level(two)
        ## Odds are equivalent
        one_winner = rand() < 0.5
    elseif level(one) > level(two)
        ## Odds are in favor of one
        one_winner = 2 * rand() > rand()
    else
        ## Odds are in favor of two
        one_winner = rand() > 2 * rand()
    end

    one_winner ? (up = one; down = two) : (up = two; down = one)

    new_lvl_up = min(level(up) + 1, 10)
    new_pos_up = clamp.(rand(-1:1, 2) .+ space(up), [1, 1], size(model.space)[1:2])
    move_agent!(up, (new_pos_up..., new_lvl_up), model)
    new_lvl_down = level(down) - 1
    if new_lvl_down == 0
        kill_agent!(down, model)
    else
        move_agent!(down, (space(down)..., new_lvl_down), model)
    end
end

# If an agent has a prisoner, it will taunt it for a time, then kill it, so long as
# an opportunist doesn't sneak up on them first!
# Here we use the tuple constructor with [`nearby_ids`](@ref) to look for agents at
# the same position as the captor `(0, 0)`, and *any* level `(..., 10)`. We could also
# use the range constructor in this instance
# `nearby_ids(agent, model, [(1, 0:0), (2, 0:0)])`, meaning which is more performant but not as
# readable.

function captor_behavior!(agent, model)
    close_ids = collect(nearby_ids(agent, model, (0, 0, 10)))
    if length(close_ids) == 1
        ## Taunt prisoner or kill it
        prisoner = model[close_ids[1]]
        if prisoner.capture_time > 10
            agent.shape = :square
            gain = ceil(Int, level(prisoner) / 2)
            new_lvl = min(level(agent) + gain, 10)
            kill_agent!(prisoner, model)
            agent.has_prisoner = false
            move_agent!(agent, (space(agent)..., new_lvl), model)
        end
    else
        ## Someone is here to kill the captor. Could be more than one opponent
        prisoner = [model[id] for id in close_ids if model[id].capture_time > 0][1]
        exploiter = rand([
            model[id]
            for
            id in close_ids if
            model[id].capture_time == 0 && model[id].has_prisoner == false
        ])
        exploiter.shape = :square
        gain = ceil(Int, level(agent) / 2)
        new_lvl = min(level(agent) + rand(1:gain), 10)
        kill_agent!(agent, model)
        move_agent!(exploiter, (space(exploiter)..., new_lvl), model)
        ## Prisoner runs away in the commotion
        prisoner.shape = :utriangle
        prisoner.capture_time = 0
        walk!(prisoner, (rand(-1:1, 2)..., 0), model)
    end
end

# When there are only few fighters standing, the stakes are higher. Prior experience is
# paramount since there is no gain, and fights are to the death.

function endgame!(agent, model)
    origin = space(agent)
    end_ids = collect(Iterators.filter(
        id -> model[id].shape == :circle && id != agent.id,
        allids(model),
    ))
    agent.shape = :circle
    if !isempty(end_ids)
        opponent = closest_target(agent, end_ids, model)
        target = space(opponent)
        if origin == target
            ## Battle
            agent.shape = :square
            opponent.shape = :square
            showdown!(agent, opponent, model)
        else
            walk!(agent, (sign.(target .- origin)..., 0), model)
        end
    end
end

function showdown!(one::Fighter, two::Fighter, model)
    if level(one) == level(two)
        ## Odds are equivalent
        one_winner = rand() < 0.5
    elseif level(one) > level(two)
        ## Odds are in favor of one
        one_winner = level(one) - level(two) * rand() > rand()
    else
        ## Odds are in favor of two
        one_winner = rand() > level(two) - level(one) * rand()
    end

    one_winner ? kill_agent!(two, model) : kill_agent!(one, model)
end

# The rest of our interactions flow down a hierarchy, so we'll place them directly in the
# `agent_step!` function. We use the tuple search for `occupied_ids` here, as we did with
# `close_ids` above. The rest of the searches however use the range search to provide a
# more precise criteria.
#
# The easiest context to explore is `worthy_ids`: all we want to do is find an agent with
# a similar level. If we used the tuple search here, we would have to search
# `(100, 100, 1)` - even though we are not at all interested in the spatial location of
# the neighbors at this time. `(3, -1:1)` is therefore more accurate representation.
#
# A more complex example is that of `strong_ids`. We are looking for agents with a level
# 2-4 points higher withing a distance of `(5, 5)`. The range search becomes a little
# verbose, but precise. An equivalent tuple search is not completely possible however.
# The closest solution is `(5, 5, 4)`, which also looks for *weaker* opponents and must
# be filtered to the correct neighbor set after the fact. In this instance the range
# search has significant performance gains.

function agent_step!(agent, model)
    if agent.capture_time > 0
        ## Captured agents are powerless, but we need to keep track of how long
        ## they have been in this state
        agent.capture_time += 1
    elseif agent.has_prisoner
        captor_behavior!(agent, model)
    else
        origin = space(agent)
        ## Find agents that have captives, they are not focused
        occupied_ids = collect(Iterators.filter(
            id -> model[id].has_prisoner,
            nearby_ids(agent, model, (7, 7, 10)),
        ))
        if !isempty(occupied_ids)
            ## Sneak up behind them
            target = space(closest_target(agent, occupied_ids, model))
            agent.shape = :pentagon
            walk!(agent, (sign.(target .- origin)..., 0), model)
        else
            ## Opponents that are greatly higher in rank that the current agent
            strong_ids = collect(nearby_ids(agent, model, [(1, -5:5), (2, -5:5), (3, 2:4)]))
            if !isempty(strong_ids)
                ## Run away from nearest
                target = space(closest_target(agent, strong_ids, model))
                agent.shape = :utriangle
                walk!(agent, (sign.(origin .- target)..., 0), model)
            else
                ## There are no distractions. Search for the closest worthy opponent
                worthy_ids = collect(nearby_ids(agent, model, [(3, -1:1)]))
                if !isempty(worthy_ids)
                    opponent = closest_target(agent, worthy_ids, model)
                    target = space(opponent)
                    if origin == target
                        ## Battle
                        agent.shape = :square
                        opponent.shape = :square
                        battle!(agent, opponent, model)
                    else
                        ## Move towards worthy opponent
                        agent.shape = :diamond
                        walk!(agent, (sign.(target .- origin)..., 0), model)
                    end
                else
                    ## Find any weak targets in the vicinity
                    weak_ids = collect(nearby_ids(
                        agent,
                        model,
                        [(1, -10:10), (2, -10:10), (3, -4:-2)],
                    ))
                    if !isempty(weak_ids)
                        prisoner = closest_target(agent, weak_ids, model)
                        target = space(prisoner)
                        if origin == target
                            ## Capture and taunt target
                            agent.has_prisoner = true
                            agent.shape = :vline
                            prisoner.capture_time += 1
                            prisoner.shape = :hline
                        else
                            ## Chase down nearest (can move 2 steps at a time!)
                            agent.shape = :star4
                            walk!(agent, (2 .* sign.(target .- origin)..., 0), model)
                        end
                    else
                        ## Abandon honour. This is the end
                        endgame!(agent, model)
                    end
                end
            end
        end
    end

    return nothing
end

# ## Let the Battle Begin
# Plotting is relatively straightforward. [`plotabm`](@ref) cannot be used explicitly (yet)
# since it expects our categorical dimension is actually a third spatial one.
# We start with some custom legends to easier understand the dynamics.

clr(agent) = cgrad(:tab10)[level(agent)]
lvl = plot(
    (1:10)',
    label = ["Level 1" "Level 2" "Level 3" "Level 4" "Level 5" "Level 6" "Level 7" "Level 8" "Level 9" "Level 10"],
    palette = :tab10,
    legend = :top,
    framestyle = :none,
)
state = scatter(
    (1:8)',
    xlims = (-2, -1),
    label = ["Battle" "Run" "Showdown" "Sneak" "Duel" "Captor" "Prisoner" "Chase"],
    markershape = [:square :utriangle :circle :pentagon :diamond :vline :hline :star4],
    color = :black,
    legend = :top,
    framestyle = :none,
)
anim = @animate for i in 0:225
    posn = [space(model[id]) for id in by_id(model)]
    cm = [clr(model[id]) for id in by_id(model)]
    shp = [model[id].shape for id in by_id(model)]
    ll = @layout [a{0.7w} [b; c]]
    battle = scatter(
        posn,
        legend = :none,
        color = cm,
        markersize = 7,
        markershape = shp,
        xlims = (-2, 103),
        ylims = (-2, 103),
        showaxis = false,
        minorgrid = true,
        ticks = (0:10:100, []),
    )
    plot(battle, lvl, state, layout = ll, size = (850, 600))
    step!(model, agent_step!, 1)
end
gif(anim, "battle.gif", fps = 10)

# Some interesting behaviour emerges: sometimes you see a group of diamonds chasing one triangle.
# What ends up happening here is usually a close pair that wishes to fight gets caught
# out by the weaker one of the two running away from an even stronger opponent. Problem
# is that this stronger opponent is chasing the stronger of the pair, but since the
# weakest of the pair is still closer to the newcomer, there is a stalemate.
# This is usually resolved by hitting a boundary or other opponents.


text:
using FiniteDiff, LinearAlgebra, SparseArrays, Test, StaticArrays

function f(x)
  xm1 = [0;x[1:end-1]]
  xp1 = [x[2:end];0]
  xm1 - 2x + xp1
end

function second_derivative_stencil(N)
  A = zeros(N,N)
  for i in 1:N, j in 1:N
      (j-i==-1 || j-i==1) && (A[i,j]=1)
      j-i==0 && (A[i,j]=-2)
  end
  A
end

x = @SVector ones(30)
J = FiniteDiff.finite_difference_jacobian(f,x, Val{:forward}, eltype(x))
@test J ≈ second_derivative_stencil(30)
_J = sparse(J)
J = FiniteDiff.finite_difference_jacobian(f,x, Val{:central}, eltype(x))
@test J ≈ second_derivative_stencil(30)

J = FiniteDiff.finite_difference_jacobian(f,x, Val{:complex}, eltype(x))
@test J ≈ second_derivative_stencil(30)

#1x1 SVector test
x = SVector{1}([1.])
f(x) = x
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:forward}, eltype(x))
@test J ≈ SMatrix{1,1}([1.])
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:central}, eltype(x))
@test J ≈ SMatrix{1,1}([1.])
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:complex}, eltype(x))
@test J ≈ SMatrix{1,1}([1.])


text:
using ExploratoryDataAnalysis
using Test

using StatsBase
using FreqTables

a = [1, 2, 3, 1, 2, 3, 1, 1, 2, 2]
b = [11, 11, 12, 12, 11, 11, 12, 12, 11, 12]

@testset "Cramer's V" begin
    @test round(cramerv(a, b), digits = 7) == 0.4472136
    @test round(cramerv(counts(a, b)), digits = 7) == 0.4472136
    @test round(cramerv(freqtable(a, b)), digits = 7) == 0.4472136
    @test round(cramerv(freqtable(a, b).array), digits = 7) == 0.4472136
end

@testset "infovalue" begin
    @test round(infovalue(counts(a, b)), digits = 7) == 0.8788898
    @test round(infovalue(freqtable(a, b)), digits = 7) == 0.8788898
    @test round(infovalue(freqtable(a, b).array), digits = 7) == 0.8788898
end

@testset "mutual info" begin
    @test round(mutualinfo(a, b), digits = 7) == 0.1046496
    @test round(mutualinfo(counts(a, b)), digits = 7) == 0.1046496
    @test round(mutualinfo(freqtable(a, b)), digits = 7) == 0.1046496
    @test round(mutualinfo(freqtable(a, b).array), digits = 7) == 0.1046496
end