text:
module example_PSD_waveform

# ----------------------------------------------------
# --- Modules
# ----------------------------------------------------
using DigitalComm
# --- External Modules
using Plots
gr();
using Printf
using FFTW

# ----------------------------------------------------
# --- Core functions
# ----------------------------------------------------
""" psdWaveform.m
---
Compute the power spectral density (i.e the spectrum here) of the signal parametrized by the waveform structure waveform, for a number of symbol nbSymb.
The frequency allocation is the one inherited from the waveform structure (i.e waveform.allocatedSubcarriers).
# --- Syntax
( freq,psd )	= psdWaveform(waveform,nbSymb,allocatedSubcarriers);
# ---  Input parameters
- waveform  : Structure associated to transmitted waveform
- nbSymb	  : Number of symbol to be transmitted [Int]
- nbIt	  : Monte carlo parameter for PSD evaluation (should be > 1)
# --- Output parameters
- freq	  : Vector of frequency evaluation (between -0.5 and 0.5). [Array{Float64,L}]
- psd		  : Spectrum evaluated on freq [Array{Complex{Float64}},L]
# --- Input parameters
-
# --- Output parameters
-
# ---
# v 1.0 - Robin Gerzaguet.
"""
function  psdWaveform(waveform,nbSymb,nbIt)
	# ----------------------------------------------------
	# --- PSD calculation
	# ----------------------------------------------------
	# --- Getting frequency allocation
	allocatedSubcarriers  = waveform.allocatedSubcarriers;
	# --- Getting number of bits
	# First, frequency size
	nbSubcarriers	      = length(allocatedSubcarriers);
	# Force a fiven mcs
	mcs				      = 4;	  # QPSK.
	# Deduce number of required bits
	nbBits			      = nbSymb * nbSubcarriers * Int(log2(mcs));
	# --- Init psd evaluator
	psd = 0;
	# --- Iterative PSD calculation
	for iN = 1 : 1 : nbIt
		# --- Binary sequence
		bitSeq	      = genBitSequence(nbBits);
		# Mapping
		qamSeq		  = bitMappingQAM(mcs,bitSeq);
		# --- T/F matrix
		qamMat		  = reshape(qamSeq,nbSubcarriers,nbSymb);
		# --- Signal
		sigPSD		  = genSig(qamMat,waveform);
		# --- Mean PSD:
		psd		  	  = psd .+ 1/nbIt*1/length(sigPSD)*abs.(fftshift(fft(sigPSD))).^2;
	end
	# --- Calculating sampling frequency
	# Returns Nyquist frequency
	fe		= 1;
	Basefe	= (0:(length(psd) .-1))./length(psd)*fe .-fe/2;
	return (Basefe,psd);
end

# ----------------------------------------------------
# --- Main routine
# ----------------------------------------------------
function main()
	# ----------------------------------------------------
	# --- Overall parameters
	# ----------------------------------------------------
	# --- Overall PHY parameters
	nbIt			= 50;			  # --- Iteration number
	nbSymb 			= 14;			  # --- Number of symbols (one frame)
	nFFT 			= 1024;			  # --- Base FFT size
	samplingFreq	= 15.36;		  # --- Frequency value (MHz)
	# --- Frequency allocation
	#allocatedSubcarriers= getLTEAlloc(nFFT);
	#allocatedSubcarriers = (1:12*4);
	# 4 RB alloc. 1 RB space. 4 RB allocated
	allocatedSubcarriers = [1:12*4; 12*5 .+ (1:12*4)];

	# ----------------------------------------------------
	# --- Waveform contender
	# ----------------------------------------------------
	# --- Init OFDM structure
	ofdm  = initOFDM(
					 nFFT,						        # --- nFFT                 : FFT size
					 72,						        # --- nCP                  : CP size
					 allocatedSubcarriers		        # --- allocatedSubcarriers : Subcarrier allocation
					 );
	# --- Init SCFDMA structure
	scfdma  = initSCFDMA(
						 nFFT,						        # --- nFFT                 : FFT size
						 72,						        # --- nCP                  : CP size
						 allocatedSubcarriers,		        # --- allocatedSubcarriers : Subcarrier allocation
						 12;								# --- sizeDFT			   : DFT preprocessing size
						 );
	# --- Init UF-OFDM structure
	ufofdm  = initUFOFDM(
						 nFFT,					        # --- nFFT                 : FFT size
						 73,					        # --- L                    : Filter length (same size +1 due to conv)
						 allocatedSubcarriers,	        # --- allocatedSubcarriers : Subcarrier allocation
						 applyPD=1,				        # --- applyPD              : Do predistortion at Tx stage
						 attenuation=40,		        # --- attenuation          : Filter attenuation in dB
						 );
	# --- Init BF-OFDM structure
	bfofdm	= initBFOFDM(
						 32,				            # --- nFBMC                : PPN size (max number of carriers)
						 64,				            # --- nOFDM                : Precoder size (OFDM sizer)
						 3,				            	# --- K                    : Overlapping factor
						 9,					            # --- GI                   : CP size of precoder
						 0.5,				            # --- δ                    : compression factor
						 allocatedSubcarriers,          # --- allocatedSubcarriers : Subcarrier allocation
						 "gaussian",	            # --- filterName           : Pulse shape name
						 BT=0.36,				        # --- BT                   : Potential BT value for Gaussian
						 filterStopBand = 110,			# --- filterStopBand       : DC stopband value
						 fS=[],				            # --- fS                   : Potential frequency coefficient for FS filter
						 nFFT= 1024,		            # --- nFFT                 : associated FFT value in Rx
						 nCP= 72,			            # --- nCP                  : extended CP size
						 );
	# --- Init WOLA-OFDM structure
	wola  = initWOLA(
					 nFFT,						        # --- nFFT                 : FFT size
					 72,						        # --- nCP                  : CP size
					 allocatedSubcarriers,		        # --- allocatedSubcarriers : Subcarrier allocation
					 "triangle",						# --- Window type @Tx side
					 20,								# --- Window size @Tx side
					 "triangle",						# --- Window type @Rx side
					 20,								# --- Window size @Rx side
					 );
	fbmc  = initFBMC(
					 nFFT,						        # --- nFFT                 : FFT size
					 4,									# --- K					   : Overlapping factor
					 allocatedSubcarriers		        # --- allocatedSubcarriers : Subcarrier allocation
					 );

	# ----------------------------------------------------
	# --- Merging structures
	# ----------------------------------------------------
	# Create  a dictionnary to rule them all
	waveforms 	= initWaveforms(ofdm,
								scfdma,
								ufofdm,
								bfofdm,
								wola,
								fbmc,
								);
	# ----------------------------------------------------
	# --- PSD main calculation
	# ----------------------------------------------------
	# --- Init plot container
	plt	  = plot(reuse=false);
	decim	= 1;	# decimation for light plots
	# --- Iterative PSD generation
	for (name,struc) in waveforms
		# --- Calculate PSD for the configuration
		(fe,psd)  = psdWaveform(struc,nbSymb,nbIt);
		# Plot the result
		plot!(plt,fe[1:decim:end].*samplingFreq,10 .* log10.(psd[1:decim:end]/maximum(psd)),label=name,legend=:topleft);
	end
	# --- Update plot and adding labels
	# Purpose is to zoom out on allocated region.
	scsN 		= (1/1024)*samplingFreq; 		# Subscarrier spacing (normalized)
	rbV 		= (12*12); 		# See several RB for psd fall-off
	ylims!(-120,5);
	xlims!(-rbV*scsN,maximum(allocatedSubcarriers)*scsN+2*12*scsN);
	xlabel!("Frequency [MHz]");
	ylabel!("Spectrum");
	display(plt)

end

end


text:
module Contexts

export withcontext, getuser

using ..Model

mutable struct Context
    user::User
end

function withcontext(f, user::User)
    task_local_storage(:CONTEXT, Context(user)) do
        f()
    end
end

function getuser()
    if haskey(task_local_storage(), :CONTEXT)
        return task_local_storage(:CONTEXT).user
    else
        throw(ArgumentError("no valid context set"))
    end
end

end # module

ext: jl

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
"""
    ngwp_analysis(G::GraphSig, wavelet_packet::Array{Float64,3})

For a GraphSig object `G`, generate the matrix of NGWP expansion coefficients.

# Input Arguments
- `G::GraphSig`: an input GraphSig object
- `wavelet_packet::Array{Float64,3}`: the varimax wavelets packet.

# Output Argument
- `dmatrix::Array{Float64,3}`: the expansion coefficients matrix.

"""
function ngwp_analysis(G::GraphSig, wavelet_packet::Array{Float64,3})
    f = G.f
    fcols = size(f, 2)
    (N, jmax, ) = Base.size(wavelet_packet)

    dmatrix = zeros(N, jmax, fcols)
    dmatrix[:, 1, :] = f

    for j = 2:jmax
        for i = 1:N
            dmatrix[i, j, :] = f' * wavelet_packet[i, j, :]
        end
    end

    return dmatrix
end

"""
    const_proj_wavelets(𝚽,vlist,elist; method = "Modified Gram-Schmidt with Lp Pivoting")

construct projection wavelets, i.e., project δ ∈ vlist onto span({φⱼ| j ∈ elist}).

# Input Arguments
- `𝚽::Matrix{Float64}`: graph Laplacian eigenvectors 𝚽
- `vlist::Array{Int}`: the list of considered node indices.
- `elist::Array{Int}`: the list of considered eigenvector indices.
- `method::Symbol`: default is `:MGSLp`. other options: `:IP` (Iterative-Projection),
    `GS` (Gram Schmidt).

# Output Argument
- `Wav::Matrix{Float64}`: a matrix whose columns are projected wavelet vectors.

"""
function const_proj_wavelets(𝚽, vlist, elist; method = :MGSLp)
    if length(vlist) == 1
        return 𝚽[:, elist]
    end
    N = size(𝚽, 1)
    m = length(vlist)
    Wav = zeros(N, m)

    B = 𝚽[:, elist]

    if method == :IP
        for k in 1:length(vlist)
            wavelet = Proj(spike(vlist[k], N), B)
            Wav[:, k] .= wavelet ./ norm(wavelet)
            B = wavelet_perp_Matrix(wavelet, B)
        end
    elseif method == :GS
        P = B * B'
        for k in 1:length(vlist)
            wavelet = P * spike(vlist[k], N)
            Wav[:,k] .= wavelet ./ norm(wavelet)
        end
        Wav, complement_dim = gram_schmidt(Wav)
        if complement_dim != 0
            complement_space = B * nullspace(Wav' * B)
            Wav = hcat(Wav, complement_space)
        end
    elseif method == :MGSLp
        P = B * B'
        for k in 1:length(vlist)
            wavelet = P * spike(vlist[k], N)
            Wav[:,k] .= wavelet ./ norm(wavelet)
        end
        Wav, complement_dim = mgslp(Wav)
        if complement_dim != 0
            complement_space = B * nullspace(Wav' * B)
            Wav = hcat(Wav, complement_space)
        end
    else
        error("Do not support method ", method)
    end

    return Wav
end

"""
    function NGWP_jkl(GP_star::GraphPart, drow::Int, dcol::Int)

Generate the (j,k,l) indices for the NGWP basis vector corresponding to the coefficient dmatrix(drow,dcol)

### Input Arguments
* `GP_star::GraphPart`: a GraphPart object of the dual grpah
* `drow::Int`: the row of the expansion coefficient
* `dcol::Int`: the column of the expansion coefficient

### Output Argument
* `j`: the level index of the expansion coefficient
* `k`: the subregion in dual graph's index of the expansion coefficient
* `l`: the tag of the expansion coefficient
"""
function NGWP_jkl(GP_star::GraphPart, drow::Int, dcol::Int)
    (j, k, l) = GHWT_jkl(GP_star, drow, dcol)
    return j,k,l
end

"""
    (dvec_ngwp, BS_ngwp) = ngwp_bestbasis(dmatrix::Array{Float64,3}, GP_star::GraphPart;
                                          cfspec::Any = 1.0, flatten::Any = 1.0,
                                          j_start::Int = 1, j_end::Int = size(dmatrix, 2),
                                          useParent::Bool = true)

Select the best basis from the matrix of NGWP expansion coefficients.

# Input Arguments
- `dmatrix::Array{Float64,3}`: the matrix of expansion coefficients
- `GP_star::GraphPart`: an input GraphPart object of the dual graph
- `cfspec::Any`: the specification of cost functional to be used (default = 1.0,
    i.e., 1-norm)
- `flatten::Any`: the method for flattening vector-valued data to scalar-valued
    data (default = 1.0, i.e, 1-norm)
- `useParent::Bool`: the flag to indicate if we update the selected best basis
    subspace to the parent when parent and child have the same cost (default = false)

#  Output Arguments
- `dvec_ngwp::Matrix{Float64}`: the vector of expansion coefficients corresponding
    to the NGWP best basis
- `BS_ngwp::BasisSpec`: a BasisSpec object which specifies the NGWP best basis
"""
function ngwp_bestbasis(dmatrix::Array{Float64,3}, GP_star::GraphPart;
                            cfspec::Any = 1.0, flatten::Any = 1.0,
                            j_start::Int = 1, j_end::Int = size(dmatrix, 2),
                            useParent::Bool = false)
    dvec_ngwp, BS_ngwp = ghwt_c2f_bestbasis(dmatrix, GP_star; cfspec = cfspec,
                            flatten = flatten, j_start = j_start, j_end = j_end,
                            useParent = useParent)
    BS_ngwp.description = "NGWP Best Basis"
    return dvec_ngwp, BS_ngwp
end


text:
Dict(
     "" => "Plural-Forms: nplurals=2; plural = n > 1",

     "hello"    => "c'a va t'il - fr",
     raw"These are $(1) houses" => [ raw"C'est $(1) maison",
                                    raw"Ce sont beaucoup($(1)) de maisons",]
    )


text:
#
# Datums
# note that geodesy itself doesn't understand datums, it just stops you doing silly things with them

abstract AbstractDatum

# start a get_datum function to retrieve the datum where possible
get_datum{T <: AbstractDatum}(::Type{T}) = T

#
# Unknown
#
immutable UnknownDatum <: AbstractDatum end

get_datum(::Type{AbstractDatum}) = UnknownDatum                     # make the abstract datum return this
get_datum{T <: GeodesyType}(::Union{Type{T}, T}) = UnknownDatum     # default answer when called on the appropriate thing

# display ???'s for unkown ellipse for compactness
ellipsoid(::Type{UnknownDatum}) = error("The ellipse is unknown")
# show(io::IO, ::Type{UnknownDatum}) = print(io, "???") # this is killing the code generation in type_methods.jl

#
# Known
#
abstract KnownDatum <: AbstractDatum

# build alias's for known datums
"""
When using a "static" datum fixed points on the Earth's surface move over time because of effects like continental drift (the datum is static but the Earth is not)
"""
abstract StaticDatum <: KnownDatum

# And make some datums based on them
immutable WGS84 <: StaticDatum; end
ellipsoid(::Type{WGS84}) = ellipsoid(WGS84_ELLIPSE)
show(io::IO, ::Type{WGS84}) = print(io, "WGS84")

immutable NAD27 <: StaticDatum; end # is this actually dynamic?
ellipsoid(NAD27) = ellipsoid(CLARKE66_ELLIPSE)
show(io::IO, ::Type{NAD27}) = print(io, "NAD27")

immutable ED50 <: StaticDatum; end
ellipsoid(::Type{ED50}) = ellipsoid(HAYFORD_ELLIPSE)
show(io::IO, ::Type{ED50}) = print(io, "ED50")

immutable OSGB36 <: StaticDatum; end # is this actually dynamic?
ellipsoid(::Type{OSGB36}) = ellipsoid(AIRY_ELLIPSE)
show(io::IO, ::Type{OSGB36}) = print(io, "OSGB36")

# The below are dynamic datums (a fixed point on the earth's surface doesn't move with continental drift)
# (aka the datum is dynamic allowing points to be static)
"""
When using a "dynamic" datum fixed points on the Earth's surface have constant coordinates (static) over time, because the datum changes (dynamic) over time to account for effects like continental drift.
Any accurate transformation between static and dynamic datums should include a time
"""
abstract DynDatum <: KnownDatum

# Australia
immutable GDA94  <:  DynDatum; end
ellipsoid(::Type{GDA94}) = ellipsoid(GRS80_ELLIPSE)
ref_date(::Type{GDA94}) = DateTime(1994)
show(io::IO, ::Type{GDA94}) = print(io, "GDA94")

# New Zealand
immutable NZGD2000  <:  DynDatum; end
ellipsoid(::Type{NZGD2000}) = ellipsoid(NZGD2000_ELLIPSE)
ref_date(::Type{NZGD2000}) = DateTime(2000)
show(io::IO, ::Type{NZGD2000}) = print(io, "NZGD2000")

# Europia
immutable ETRS89  <:  DynDatum; end # N.B. adopting the same notation abuse as everyone else, this is actually ETRF89
ellipsoid(::Type{ETRS89}) = ellipsoid(GRS80_ELLIPSE)
ref_date(::Type{ETRS89}) = DateTime(1989)
show(io::IO, ::Type{ETRS89}) = print(io, "ETRS89")

# Americania
immutable NAD83  <:  DynDatum; end
ellipsoid(::Type{NAD83}) = ellipsoid(GRS80_ELLIPSE)
show(io::IO, ::Type{NAD83}) = print(io, "NAD83")
ref_date(::Type{NAD83}) = DateTime(1983)

#
# type for custom geoids
#
abstract AbstractGeoid
get_geoid{T <: AbstractGeoid}(::Type{T}) = T

"""
Function to set / get the geoid directory
"""
function set_geoid_dir(path::AbstractString)
    Geodesy.geodesy_properties.geoid_dir = path
end
get_geoid_dir() = Geodesy.geodesy_properties.geoid_dir

# dont know the Geoid
"""
Default unknown geoid
"""
immutable UnknownGeoid <: AbstractGeoid end
geoid_file(::Type{UnknownGeoid}) = error("Unsure which Geoid to use for an unknown Geoid")

# Known Geoids
abstract KnownGeoid <: AbstractGeoid

# Australia
immutable AusGeoid09 <: KnownGeoid end
geoid_file(::Type{AusGeoid09}) = "ausgeoid09.gtx"

# New Zealand
immutable NZGeoid2016 <: KnownGeoid end
geoid_file(::Type{NZGeoid2016}) = "nzgeoid2016.gtx"

# UK
immutable OSGM02 <: KnownGeoid end
geoid_file(::Type{OSGM02}) = "osgm02etrs89.gtx"

# Google Earth
immutable EGM96 <: KnownGeoid end
geoid_file(::Type{EGM96}) = "egm96-5.gtx"

# Ellipsoidal (a.k.a no geoid)
immutable NoGeoid <: KnownGeoid end
geoid_file(::Type{NoGeoid}) = ""

# function to get a list of all datums
function get_datums(datum::DataType=AbstractDatum, super::ASCIIString="", out=Vector{ASCIIString}(0))
    sub_datums = subtypes(datum)
    if (length(sub_datums) > 0)
        for sd in sub_datums
            out = get_datums(sd, datum==AbstractDatum ? "" : " (" * string(datum) * ")", out)
        end
    else
        push!(out, replace(string(datum) *  super  , "Geodesy.", ""))
    end
    return out
end

# function to try and guess the geoid dir
find_geoid_dir() = get(ENV, "GEOID_DIR", "")

#####################################################################
# dev tool, find stuff in the Proj4 dicts
# e.g find_match(Proj4.epsg, [r"proj=longlat", r"datum=WGS84"])
#####################################################################
function find_match{T,U}(p4_dict::Dict{T,U}, exprs)

    if !isa(exprs, Vector)
        exprs = [exprs]
    end

    matched_key = Vector{T}(0)
    matched_str = Vector{U}(0)
    for key in keys(p4_dict)
        str = p4_dict[key]
        hasmatch = true
        for expr in exprs
            hasmatch &= length(matchall(expr, str)) > 0
        end
        if hasmatch
            push!(matched_key, key)
            push!(matched_str, str)
        end
    end
    return [matched_key   matched_str]
end


text:
# !! auto-generated code, do not edit !!
module sensor_msgs
using RobotOSData.Messages
using RobotOSData.StdMsgs
using RobotOSData.CommonMsgs
struct MagneticField <: Readable
    header::Header
    magnetic_field::geometry_msgs.Vector3
    magnetic_field_covariance::SVector{9, Float64}
end
struct NavSatStatus <: Readable
    status::Int8
    service::UInt16
end
struct LaserScan <: Readable
    header::Header
    angle_min::Float32
    angle_max::Float32
    angle_increment::Float32
    time_increment::Float32
    scan_time::Float32
    range_min::Float32
    range_max::Float32
    ranges::Vector{Float32}
    intensities::Vector{Float32}
end
struct BatteryState <: Readable
    header::Header
    voltage::Float32
    current::Float32
    charge::Float32
    capacity::Float32
    design_capacity::Float32
    percentage::Float32
    power_supply_status::UInt8
    power_supply_health::UInt8
    power_supply_technology::UInt8
    present::Bool
    cell_voltage::Vector{Float32}
    location::String
    serial_number::String
end
struct JoyFeedback <: Readable
    type_::UInt8
    id::UInt8
    intensity::Float32
end
struct JoyFeedbackArray <: Readable
    array::Vector{JoyFeedback}
end
struct Temperature <: Readable
    header::Header
    temperature::Float64
    variance::Float64
end
struct PointField <: Readable
    name::String
    offset::UInt32
    datatype::UInt8
    count::UInt32
end
struct PointCloud2 <: Readable
    header::Header
    height::UInt32
    width::UInt32
    fields::Vector{PointField}
    is_bigendian::Bool
    point_step::UInt32
    row_step::UInt32
    data::Vector{UInt8}
    is_dense::Bool
end
struct Illuminance <: Readable
    header::Header
    illuminance::Float64
    variance::Float64
end
struct LaserEcho <: Readable
    echoes::Vector{Float32}
end
struct MultiEchoLaserScan <: Readable
    header::Header
    angle_min::Float32
    angle_max::Float32
    angle_increment::Float32
    time_increment::Float32
    scan_time::Float32
    range_min::Float32
    range_max::Float32
    ranges::Vector{LaserEcho}
    intensities::Vector{LaserEcho}
end
struct ChannelFloat32 <: Readable
    name::String
    values::Vector{Float32}
end
struct PointCloud <: Readable
    header::Header
    points::Vector{geometry_msgs.Point32}
    channels::Vector{ChannelFloat32}
end
struct NavSatFix <: Readable
    header::Header
    status::NavSatStatus
    latitude::Float64
    longitude::Float64
    altitude::Float64
    position_covariance::SVector{9, Float64}
    position_covariance_type::UInt8
end
struct TimeReference <: Readable
    header::Header
    time_ref::ROSTime
    source::String
end
struct Imu <: Readable
    header::Header
    orientation::geometry_msgs.Quaternion
    orientation_covariance::SVector{9, Float64}
    angular_velocity::geometry_msgs.Vector3
    angular_velocity_covariance::SVector{9, Float64}
    linear_acceleration::geometry_msgs.Vector3
    linear_acceleration_covariance::SVector{9, Float64}
end
struct RegionOfInterest <: Readable
    x_offset::UInt32
    y_offset::UInt32
    height::UInt32
    width::UInt32
    do_rectify::Bool
end
struct CompressedImage <: Readable
    header::Header
    format::String
    data::Vector{UInt8}
end
struct Image <: Readable
    header::Header
    height::UInt32
    width::UInt32
    encoding::String
    is_bigendian::UInt8
    step::UInt32
    data::Vector{UInt8}
end
struct Joy <: Readable
    header::Header
    axes::Vector{Float32}
    buttons::Vector{Int32}
end
struct Range <: Readable
    header::Header
    radiation_type::UInt8
    field_of_view::Float32
    min_range::Float32
    max_range::Float32
    range::Float32
end
struct FluidPressure <: Readable
    header::Header
    fluid_pressure::Float64
    variance::Float64
end
struct JointState <: Readable
    header::Header
    name::Vector{String}
    position::Vector{Float64}
    velocity::Vector{Float64}
    effort::Vector{Float64}
end
struct MultiDOFJointState <: Readable
    header::Header
    joint_names::Vector{String}
    transforms::Vector{geometry_msgs.Transform}
    twist::Vector{geometry_msgs.Twist}
    wrench::Vector{geometry_msgs.Wrench}
end
struct RelativeHumidity <: Readable
    header::Header
    relative_humidity::Float64
    variance::Float64
end
struct CameraInfo <: Readable
    header::Header
    height::UInt32
    width::UInt32
    distortion_model::String
    D::Vector{Float64}
    K::SVector{9, Float64}
    R::SVector{9, Float64}
    P::SVector{12, Float64}
    binning_x::UInt32
    binning_y::UInt32
    roi::RegionOfInterest
end
end


text:
@doc raw"""
Axisymmetric tokamak equilibrium in (x,y,z) coordinates with covariant
components of the vector potential given by
```math
A (x,y,z) = \frac{1}{2} \frac{B_0}{q_0} \, \bigg( \frac{q_0 R_0 x z - r^2 y}{R^2} , \, \frac{q_0 R_0 y z + r^2 x}{R^2} , \, - q_0 R_0 \, \ln \bigg( \frac{R}{R_0} \bigg) \bigg)^T ,
```
resulting in the magnetic field with covariant components
```math
B (x,y,z) = \frac{B_0}{q_0} \, \bigg( - \frac{q_0 R_0 y + x z}{R^2} , \, \frac{q_0 R_0 x - y z}{R^2} , \, \frac{R - R_0}{R} \bigg)^T ,
```
where $R = \sqrt{ x^2 + y^2 }$ and $r = \sqrt{ (R - R_0)^2 + z^2 }$.

Parameters:
* `R₀`: position of magnetic axis
* `B₀`: B-field at magnetic axis
* `q₀`: safety factor at magnetic axis
"""
module AxisymmetricTokamakCartesian

    using RecipesBase

    import ..ElectromagneticFields
    import ..ElectromagneticFields: CartesianEquilibrium, code
    import ..AnalyticCartesianField: X, Y, Z

    export  AxisymmetricTokamakCartesianEquilibrium

    const DEFAULT_R₀ = 1
    const DEFAULT_B₀ = 1
    const DEFAULT_q₀ = 2

    const ITER_R₀ = 6.2
    const ITER_B₀ = 5.3
    const ITER_q₀ = √2

    struct AxisymmetricTokamakCartesianEquilibrium{T <: Number} <: CartesianEquilibrium
        name::String
        R₀::T
        B₀::T
        q₀::T

        function AxisymmetricTokamakCartesianEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number
            new("AxisymmetricTokamakCartesianEquilibrium", R₀, B₀, q₀)
        end
    end

    AxisymmetricTokamakCartesianEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakCartesianEquilibrium{T}(R₀, B₀, q₀)

    function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀)
        AxisymmetricTokamakCartesianEquilibrium(R₀, B₀, q₀)
    end

    function ITER()
        AxisymmetricTokamakCartesianEquilibrium(ITER_R₀, ITER_B₀, ITER_q₀)
    end

    macro code(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀)
        code(init(R₀, B₀, q₀); escape=true)
    end

    macro code_iter()
        code(ITER(); escape=true)
    end

    function Base.show(io::IO, equ::AxisymmetricTokamakCartesianEquilibrium)
        print(io, "Axisymmetric Tokamak Equilibrium in (x,y,z) Coordinates with\n")
        print(io, "  R₀ = ", equ.R₀, "\n")
        print(io, "  B₀ = ", equ.B₀, "\n")
        print(io, "  q₀ = ", equ.q₀)
    end

    R²(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = X(x,equ)^2 + Y(x,equ)^2
    r²(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2
    R(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = sqrt(R²(x,equ))
    r(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = sqrt(r²(x,equ))
    θ(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀)
    ϕ(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = atan(Y(x,equ), X(x,equ))

    ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = + equ.B₀ * (equ.R₀ * X(x,equ) * Z(x,equ) - r²(x,equ) * Y(x,equ) / equ.q₀ ) / R²(x,equ) / 2
    ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = + equ.B₀ * (equ.R₀ * Y(x,equ) * Z(x,equ) + r²(x,equ) * X(x,equ) / equ.q₀ ) / R²(x,equ) / 2
    ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakCartesianEquilibrium) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2

    ElectromagneticFields.get_functions(::AxisymmetricTokamakCartesianEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, R²=R², r²=r²)

    @recipe function f(equ::AxisymmetricTokamakCartesianEquilibrium;
                       nx = 100, ny = 120, levels = 50, size = (400,400),
                       xlims = (  0.5 * equ.R₀,   1.5 * equ.R₀),
                       ylims = (- 0.5 * equ.R₀, + 0.5 * equ.R₀))

        xgrid = LinRange(xlims[1], xlims[2], nx)
        zgrid = LinRange(ylims[1], ylims[2], ny)
        pot   = [ElectromagneticFields.A₂([xgrid[i], 0.0, zgrid[j]], equ) for i in eachindex(xgrid), j in eachindex(zgrid)]

        seriestype   := :contour
        aspect_ratio := :equal
        size   := size
        xlims  := xlims
        ylims  := ylims
        levels := levels
        legend := :none

        (xgrid, zgrid, pot')
    end

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*-------------------------------------------------------------------
#* Author: Rafael de Pelegrini Soares
#*         Maur좩cio Carvalho Maciel
#* $Id$
#*--------------------------------------------------------------------
type mixer_n
	mixer_n()=begin
		NComp=outers.NComp
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Inlet Streams",
				:Lower=>1,
				:Default=>2
			)),
			fill(stream(Dict{Symbol,Any}(
				:Brief=>"Inlet streams",
				:PosX=>0,
				:PosY=>0.5,
				:Symbol=>"_{inMix}"
			)),(Ninlet)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.5059,
				:Symbol=>"_{out}"
			)),
			[
				:(Outlet.F = sum(Inlet.F)),
				:(Outlet.F*Outlet.z([1:NComp]) = sum(Inlet.F*Inlet.z([1:NComp]))),
				:(Outlet.F*Outlet.h = sum(Inlet.F*Inlet.h)),
				:(Outlet.P = min(Inlet.P)),
			],
			[
				"Flow","Composition","Energy Balance","Pressure",
			],
			[:NComp,:Ninlet,],
			[:Inlet,:Outlet,]
		)
	end
	NComp::DanaInteger
	Ninlet::DanaInteger
	Inlet::Array{stream}
	Outlet::streamPH
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export mixer_n
function setEquationFlow(in::mixer_n)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
end
function atributes(in::mixer_n,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/mixer"
	fields[:Brief]="Model of a generic mixer (NOT Handled By the GUI) "
	fields[:Info]="== Assumptions ==
* static
* adiabatic

== Specify ==
* the inlet streams"
	drive!(fields,_)
	return fields
end
mixer_n(_::Dict{Symbol,Any})=begin
	newModel=mixer_n()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
a = Iterable(i for i ∈ [1,2,3])
b = Iterable([4,5,6])

# MonoidAlternative
# =================

@test collect(a ⊕ b) == [collect(a); collect(b)]

# FunctorApplicativeMonad
# =======================

@test map(x -> x*x, a) |> collect == [1, 4, 9]

product_ab = [5, 6, 7, 6, 7, 8, 7, 8, 9]
@test mapn(a, b) do x, y
  x + y
end |> collect == product_ab

h = @syntax_flatmap begin
  x = a
  y = b
  @pure x + y
end
@test h |> collect == product_ab

# Everything Iterable works with Iterable
h2 = @syntax_flatmap begin
  x = a
  y = isodd(x) ? Option(x*x) : Option()
  @pure x + y
end
@test collect(h2) == [2, 12]

# However note that the following does not work
h3() = @syntax_flatmap begin
  x = a
  y = isodd(x) ? Option(x*x) : Option()
  z = b
  @pure x, y, z
end
@test collect(h3()) == [
  (x, x*x, z) for x in a for z in b if isodd(x)
]

# FlipTypes
# =========

it = Iterable(Option(i) for i ∈ [1, 4, 7])
@test map(collect, flip_types(it)) == Option([1, 4, 7])


text:
using RedPitayaDAQServer
using PyPlot

# obtain the URL of the RedPitaya
include("config.jl")

rp = RedPitayaCluster([URLs[1]])

dec = 8
modulus = 4800
base_frequency = 125000000
periods_per_step = 10
samples_per_period = div(modulus, dec)#*periods_per_step
periods_per_frame = 100 # about 0.5 s frame length
frame_period = dec*samples_per_period*periods_per_frame / base_frequency

slow_dac_steps_per_frame = div(periods_per_frame, periods_per_step)

@info samples_per_period, periods_per_frame, slow_dac_steps_per_frame
decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)
passPDMToFastDAC(master(rp), true)
slowDACStepsPerFrame(master(rp), slow_dac_steps_per_frame)

numSlowDACChan(master(rp), 2)
lutA = collect(range(0,0.3,length=slow_dac_steps_per_frame))
lutB = collect(ones(slow_dac_steps_per_frame))

lutEnableDACA = ones(Bool, slow_dac_steps_per_frame)
lutEnableDACA[1:2:end] .= false
lutEnableDACB = ones(Bool, slow_dac_steps_per_frame)
lutEnableDACB[2:2:end] .= false

modeDAC(rp, "STANDARD")
frequencyDAC(rp,1,1, base_frequency / modulus)

freq = frequencyDAC(rp,1,1)
println(" frequency = $(freq)")
jumpSharpnessDAC(rp, 1, 0.01) # controls the sharpness of the jump for the square
signalTypeDAC(rp, 1 , "SQUARE")
phaseDAC(rp, 1, 1, 0.0 ) # Phase has to be given in between 0 and 1
amplitudeDACNext(rp, 1, 1, 0.1)

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "INTERNAL")
masterTrigger(rp, false)
startADC(rp)
masterTrigger(rp, true)

sleep(0.1)

enableDACLUT(master(rp), collect( cat(lutEnableDACA,lutEnableDACB,dims=2)' ) )

lut = cat(lutA,lutB*0.1,dims=2)'
setSlowDACLUT(master(rp), collect(lut))

currFr = enableSlowDAC(master(rp), true, 1, frame_period, 0.5)
uCurrentPeriod = readData(rp, currFr, 1)

sleep(0.5)

lut = cat(lutA,lutB*0.4,dims=2)'
setSlowDACLUT(master(rp), collect(lut))

currFr = enableSlowDAC(master(rp), true, 1, frame_period, 0.5)
uCurrentPeriod2 = readData(rp, currFr, 1)

figure(1)
clf()
subplot(2,1,1)
plot(vec(uCurrentPeriod[:,1,:,:]))
plot(vec(uCurrentPeriod[:,2,:,:]))
legend(("Rx1", "Rx2"))
subplot(2,1,2)
plot(vec(uCurrentPeriod2[:,1,:,:]))
plot(vec(uCurrentPeriod2[:,2,:,:]))
legend(("Rx1", "Rx2"))

savefig("images/slowDACMultiChannelWithOffset.png")

stopADC(rp)

max_stars_repo_name: Fypsilonn/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer


text:
using Statistics, LinearAlgebra, PDMats
import Distributions: IsoNormal, DiagNormal, FullNormal, logpdf
import PDMats: ScalMat, PDiagMat, PDMat

export SphericalGaussian, DiagonalGaussian, FullGaussian,
initVariances!, initMixtures!,lpdf,updateVariances!

abstract type AbstractGaussian <: AbstractMixture end

mutable struct SphericalGaussian{T <:Number} <: AbstractGaussian
        μ  ::Union{Array{T,1},Nothing}
        σ² ::Union{T,Nothing}
        #SphericalGaussian(;μ::Union{Array{T,1},Nothing},σ²::Union{T,Nothing}) where {T} = SphericalGaussian(μ,σ²)
        """SphericalGaussian(μ,σ²) - Spherical Gaussian mixture with mean μ and (single) variance σ²"""
        SphericalGaussian(μ::Union{Array{T,1},Nothing},σ²::Union{T,Nothing}=nothing) where {T} = new{T}(μ,σ²)
        SphericalGaussian(type::Type{T}=Float64) where {T} = new{T}(nothing, nothing)
end

mutable struct DiagonalGaussian{T <:Number} <: AbstractGaussian
    μ::Union{Array{T,1},Nothing}
    σ²::Union{Array{T,1},Nothing}
    """DiagonalGaussian(μ,σ²) - Gaussian mixture with mean μ and variances σ² (and fixed zero covariances)"""
    DiagonalGaussian(μ::Union{Array{T,1},Nothing},σ²::Union{Array{T,1},Nothing}=nothing) where {T} = new{T}(μ,σ²)
    DiagonalGaussian(::Type{T}=Float64) where {T} = new{T}(nothing, nothing)
end

mutable struct FullGaussian{T <:Number} <: AbstractGaussian
    μ::Union{Array{T,1},Nothing}
    σ²::Union{Array{T,2},Nothing}
    """FullGaussian(μ,σ²) - Gaussian mixture with mean μ and variance/covariance matrix σ²"""
    FullGaussian(μ::Union{Array{T,1},Nothing},σ²::Union{Array{T,2},Nothing}=nothing) where {T} = new{T}(μ,σ²)
    FullGaussian(::Type{T}=Float64) where {T} = new{T}(nothing, nothing)
end

function initVariances!(mixtures::Array{T,1}, X; minVariance=0.25, minCovariance=0.0,rng = Random.GLOBAL_RNG) where {T <: SphericalGaussian}
    (N,D)         = size(X)
    K             = length(mixtures)
    varX_byD      = fill(0.0,D)
    for d in 1:D
      varX_byD[d] = var(skipmissing(X[:,d]))
    end
    varX = max(minVariance,mean(varX_byD)/K^2)

    for (i,m) in enumerate(mixtures)
        if isnothing(m.σ²)
            m.σ² = varX
        end
    end
end

function initVariances!(mixtures::Array{T,1}, X; minVariance=0.25, minCovariance=0.0,rng = Random.GLOBAL_RNG) where {T <: DiagonalGaussian}
    (N,D)         = size(X)
    K             = length(mixtures)
    varX_byD      = fill(0.0,D)
    for d in 1:D
      varX_byD[d] = max(minVariance, var(skipmissing(X[:,d])))
    end

    for (i,m) in enumerate(mixtures)
        if isnothing(m.σ²)
            m.σ² = varX_byD
        end
    end

end

function initVariances!(mixtures::Array{T,1}, X; minVariance=0.25, minCovariance=0.0,rng = Random.GLOBAL_RNG) where {T <: FullGaussian}
    (N,D) = size(X)
    K = length(mixtures)
    varX_byD = fill(0.0,D)

    for d in 1:D
      varX_byD[d] = max(minVariance, var(skipmissing(X[:,d])))
    end

    for (i,m) in enumerate(mixtures)
        if isnothing(m.σ²)
            m.σ² = fill(0.0,D,D)
            for d1 in 1:D
                for d2 in 1:D
                    if d1 == d2
                        m.σ²[d1,d2] = varX_byD[d1]
                    else
                        m.σ²[d1,d2] = minCovariance
                    end
                end
            end
        end
    end
end

"""
    initMixtures!(mixtures::Array{T,1}, X; minVariance=0.25, minCovariance=0.0, initStrategy="grid",rng=Random.GLOBAL_RNG)

 The parameter `initStrategy` can be `grid`, `kmeans` or `given`:
 - `grid`: Uniformly cover the space observed by the data
 - `kmeans`: Use the kmeans algorithm. If the data contains missing values, a first run of `predictMissing` is done under init=`grid` to impute the missing values just to allow the kmeans algorithm. Then the em algorithm is used with the output of kmean as init values.
 - `given`: Leave the provided set of initial mixtures

"""
function initMixtures!(mixtures::Array{T,1}, X; minVariance=0.25, minCovariance=0.0, initStrategy="grid",rng = Random.GLOBAL_RNG) where {T <: AbstractGaussian}
    # debug..
    #X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing 2; 3.3 38; missing -2.3; 5.2 -2.4]
    #mixtures = [SphericalGaussian() for i in 1:3]
    # ---
    if initStrategy == "given"
        return
    end

    (N,D) = size(X)
    K     = length(mixtures)

    # count nothing mean mixtures
    nMM = 0
    for (i,m) in enumerate(mixtures)
        if isnothing(m.μ)
            nMM += 1
        end
    end

    if initStrategy == "grid"

        minX = fill(-Inf,D)
        maxX = fill(Inf,D)

        for d in 1:D
           minX[d]  = minimum(skipmissing(X[:,d]))
           maxX[d]  = maximum(skipmissing(X[:,d]))
        end

        rangedμ = zeros(nMM,D)
        for d in 1:D
            rangedμ[:,d] = collect(range(minX[d] + (maxX[d]-minX[d])/(nMM*2) , stop=maxX[d] - (maxX[d]-minX[d])/(nMM*2)  , length=nMM))
            # ex: rangedμ[:,d] = collect(range(minX[d], stop=maxX[d], length=nMM))
        end

        j = 1
        for m in mixtures
           if isnothing(m.μ)
               m.μ = rangedμ[j,:]
               j +=1
           end
        end

    elseif initStrategy == "kmeans"
        if !any(ismissing.(X)) # there are no missing
            kmμ = kmeans(X,K,rng=rng)[2]
            for (k,m) in enumerate(mixtures)
               if isnothing(m.μ)
                   m.μ = kmμ[k,:]
               end
            end
        else # missings are present
            # First pass of predictMissing using initStrategy=grid
            emOut1 = predictMissing(X,K;mixtures=mixtures,verbosity=NONE,minVariance=minVariance,minCovariance=minCovariance,initStrategy="grid",rng=rng,maxIter=10)
            kmμ = kmeans(emOut1.X̂,K,rng=rng)[2]
            for (k,m) in enumerate(mixtures)
               if isnothing(m.μ)
                   m.μ = kmμ[k,:]
               end
            end
        end
    else
        @error "initStrategy $initStrategy not supported by this mixture type"
    end

    initVariances!(mixtures,X,minVariance=minVariance, minCovariance=minCovariance,rng=rng)

end

"""lpdf(m::SphericalGaussian,x,mask) - Log PDF of the mixture given the observation `x`"""
function lpdf(m::SphericalGaussian,x,mask)
    μ  = m.μ[mask]
    σ² = m.σ²
    #d = IsoNormal(μ,ScalMat(length(μ),σ²))
    #return logpdf(d,x)
    return (- (length(x)/2) * log(2π*σ²)  -  norm(x-μ)^2/(2σ²))
end

"""lpdf(m::DiagonalGaussian,x,mask) - Log PDF of the mixture given the observation `x`"""
function lpdf(m::DiagonalGaussian,x,mask)
    μ  = m.μ[mask]
    σ² = m.σ²[mask]
    d = DiagNormal(μ,PDiagMat(σ²))
    return logpdf(d,x)
end

"""lpdf(m::FullGaussian,x,mask) - Log PDF of the mixture given the observation `x`"""
function lpdf(m::FullGaussian,x,mask)
    μ   = m.μ[mask]
    nmd = length(μ)
    σ²  = reshape(m.σ²[mask*mask'],(nmd,nmd))
    σ² = σ² + max(0, -2minimum(eigvals(σ²))) * I # Improve numerical stability https://stackoverflow.com/q/57559589/1586860 (-2 * minimum...) https://stackoverflow.com/a/35612398/1586860

    #=
    try
      d      = FullNormal(μ,PDMat(σ²))
      return logpdf(d,x)
    catch
        println(σ²)
        println(mask)
        println(μ)
        println(x)
        println(σ²^(-1))
        error("Failed PDMat")
    end
    =#
    diff = x .- μ
    return -(nmd/2)*log(2pi)-(1/2)log(det(σ²))-(1/2)*diff'*σ²^(-1)*diff

end

npar(mixtures::Array{T,1}) where {T <: SphericalGaussian} = length(mixtures) * length(mixtures[1].μ) + length(mixtures) # K * D + K
npar(mixtures::Array{T,1}) where {T <: DiagonalGaussian}  = length(mixtures) * length(mixtures[1].μ) + length(mixtures) * length(mixtures[1].μ) # K * D + K * D
npar(mixtures::Array{T,1}) where {T <: FullGaussian} = begin K = length(mixtures); D = length(mixtures[1].μ); K * D + K * (D^2+D)/2 end

function updateVariances!(mixtures::Array{T,1}, X, pₙₖ; minVariance=0.25, minCovariance = 0.0) where {T <: SphericalGaussian}

    # debug stuff..
    #X = [1 10 20; 1.2 12 missing; 3.1 21 41; 2.9 18 39; 1.5 15 25]
    #m1 = SphericalGaussian(μ=[1.0,15,21],σ²=5.0)
    #m2 = SphericalGaussian(μ=[3.0,20,30],σ²=10.0)
    #mixtures= [m1,m2]
    #pₙₖ = [0.9 0.1; 0.8 0.2; 0.1 0.9; 0.1 0.9; 0.4 0.6]
    #Xmask = [true true true; true true false; true true true; true true true; true true true]
    #minVariance=0.25
    # ---

    (N,D) = size(X)
    K = length(mixtures)
    Xmask     =  .! ismissing.(X)
    XdimCount = sum(Xmask, dims=2)

    #    #σ² = [sum([pⱼₓ[n,j] * norm(X[n,:]-μ[j,:])^2 for n in 1:N]) for j in 1:K ] ./ (nⱼ .* D)
    for k in 1:K
        nom = 0.0
        den = dot(XdimCount,pₙₖ[:,k])
        m = mixtures[k]
        for n in 1:N
            if any(Xmask[n,:])
                nom += pₙₖ[n,k] * norm(X[n,Xmask[n,:]]-m.μ[Xmask[n,:]])^2
            end
        end
        if(den> 0 && (nom/den) > minVariance)
            m.σ² = nom/den
        else
            m.σ² = minVariance
        end
    end

end

function updateVariances!(mixtures::Array{T,1}, X, pₙₖ; minVariance=0.25, minCovariance = 0.0) where {T <: DiagonalGaussian}
    # debug stuff..
    #X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4]
    #m1 = DiagonalGaussian([1.0,10.0],[5.0,5.0])
    #m2 = DiagonalGaussian([4.0,40.0],[10.0,10.0])
    #m3 = DiagonalGaussian([4.0,-2.0],[5.0,5.0])
    #mixtures= [m1,m2,m3]
    #pₙₖ = [0.9 0.1 0; 0.7 0.1 0.1; 0.8 0.2 0; 0.7 0.3 0; 0.1 0.9 0; 0.4 0.4 0.2; 0.1 0.9 0; 0.2 0.1 0.7 ; 0 0.1 0.9]
    #minVariance=0.25
    # ---

    (N,D) = size(X)
    K = length(mixtures)
    Xmask     =  .! ismissing.(X)
    #XdimCount = sum(Xmask, dims=2)

    #    #σ² = [sum([pⱼₓ[n,j] * norm(X[n,:]-μ[j,:])^2 for n in 1:N]) for j in 1:K ] ./ (nⱼ .* D)
    for k in 1:K
        m = mixtures[k]
        for d in 1:D
            nom = 0.0
            den = 0.0
            for n in 1:N
                if Xmask[n,d]
                    nom += pₙₖ[n,k] * (X[n,d]-m.μ[d])^2
                    den += pₙₖ[n,k]
                end
            end
            if(den > 0 )
                m.σ²[d] = max(nom/den,minVariance)
            else
                m.σ²[d] = minVariance
            end
        end
    end

end

function updateVariances!(mixtures::Array{T,1}, X, pₙₖ; minVariance=0.25, minCovariance = 0.0) where {T <: FullGaussian}

    # debug stuff..
    #X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4]
    #m1 = FullGaussian([1.0,10.0],[5.0 1; 1.0 5.0])
    #m2 = FullGaussian([4.0,40.0],[10.0 1.0; 1.0 10.0])
    #m3 = FullGaussian([4.0,-2.0],[5.0 1; 1.0 5.0])
    #mixtures= [m1,m2,m3]
    #pₙₖ = [0.9 0.1 0; 0.7 0.1 0.1; 0.8 0.2 0; 0.7 0.3 0; 0.1 0.9 0; 0.4 0.4 0.2; 0.1 0.9 0; 0.2 0.1 0.7 ; 0 0.1 0.9]
    #minVariance=0.25
    # ---

    (N,D) = size(X)
    K = length(mixtures)

    # NDDMAsk is true only if both (N,D1) and (N,D2) are nonmissing values
    NDDMask = fill(false,N,D,D)
    for n in 1:N
        for d1 in 1:D
            for d2 in 1:D
                if !ismissing(X[n,d1]) && !ismissing(X[n,d2])
                    NDDMask[n,d1,d2] = true
                end
            end
        end
    end

    #    #σ² = [sum([pⱼₓ[n,j] * norm(X[n,:]-μ[j,:])^2 for n in 1:N]) for j in 1:K ] ./ (nⱼ .* D)
    for k in 1:K
        m = mixtures[k]
        for d2 in 1:D # out var matrix col
            for d1 in 1:D # out var matrix row
                if d1 >= d2 # lower half of triang
                    nom = 0.0
                    den = 0.0
                    for n in 1:N
                        if NDDMask[n,d1,d2]
                            nom += pₙₖ[n,k] * (X[n,d1]-m.μ[d1])*(X[n,d2]-m.μ[d2])
                            den += pₙₖ[n,k]
                        end
                    end
                    if(den > 0 )
                        if d1 == d2
                            m.σ²[d1,d2] = max(nom/den,minVariance)
                        else
                            m.σ²[d1,d2] = max(nom/den,minCovariance)
                        end
                    else
                        if d1 == d2
                           m.σ²[d1,d2] = minVariance
                        else
                          #m.σ²[d1,d2] = minVariance-0.01 # to avoid singularity in all variances equal to minVariance
                          m.σ²[d1,d2] = minCovariance
                        end
                    end
                else # upper half of the matrix
                    m.σ²[d1,d2] = m.σ²[d2,d1]
                end
            end
        end
    end

end

"""
updateParameters!(mixtures::Array{T,1}, X, pₙₖ; minVariance=0.25, minCovariance)

Find and set the parameters that maximise the likelihood (m-step in the EM algorithm)

"""
#https://github.com/davidavdav/GaussianMixtures.jl/blob/master/src/train.jl
function updateParameters!(mixtures::Array{T,1}, X, pₙₖ; minVariance=0.25, minCovariance = 0.0) where {T <: AbstractGaussian}
    # debug stuff..
    #X = [1 10 20; 1.2 12 missing; 3.1 21 41; 2.9 18 39; 1.5 15 25]
    #m1 = SphericalGaussian(μ=[1.0,15,21],σ²=5.0)
    #m2 = SphericalGaussian(μ=[3.0,20,30],σ²=10.0)
    #mixtures= [m1,m2]
    #pₙₖ = [0.9 0.1; 0.8 0.2; 0.1 0.9; 0.1 0.9; 0.4 0.6]
    #Xmask = [true true true; true true false; true true true; true true true; true true true]

    (N,D) = size(X)
    K = length(mixtures)
    Xmask     =  .! ismissing.(X)

    #nₖ = sum(pₙₖ,dims=1)'
    #n  = sum(nₖ)
    #pₖ = nₖ ./ n

    nkd = fill(0.0,K,D)
    #nkd = [sum(pₙₖ[Xmask[:,d],k]) for k in 1:K, d in 1:D] # number of point associated to a given mixture for a specific dimension

    # updating μ...
    for k in 1:K
        m = mixtures[k]
        for d in 1:D
            nkd[k,d] = sum(pₙₖ[Xmask[:,d],k])
            if nkd[k,d] > 1
                m.μ[d] = sum(pₙₖ[Xmask[:,d],k] .* X[Xmask[:,d],d])/nkd[k,d]
            end
        end
    end

    updateVariances!(mixtures, X, pₙₖ; minVariance=minVariance, minCovariance=minCovariance)
end


text:
using Interpolations, DataStructures, Lasso, DataFrames, Flux
using LinearAlgebra: cond
# include("../src/sindyc.jl")
# include("../src/dmdc.jl");
# include("../utils/combinatorics_utils.jl");
# include("../utils/regression_utils.jl")
# include("../utils/sparse_regression_functions.jl")

#####
##### Fitting function; effectively the default constructor
#####
# """
# Naive implementation of Sparse Identification of Nonlinear DYnamics
# with control (SINDYc). To actually use sparsity, the Lasso.jl package
# is required
#
# There are several nonlinear library terms that can be implemented via passing
# a list of strings and arguments via 'library', but a custom function can also be
# passed using 'custom_func'. Currently implemented library terms are:
#     ["cross_terms", order::Int]
#         Here, 'order' is how high of an order to do
#
# model = sindyc(X, X_grad=nothing, U=nothing, ts=nothing;
#                 library=Dict(),
#                 use_lasso=false,
#                 hard_threshold=nothing,
#                 quantile_threshold=0.1,
#                 var_names = ["x", "y", "z"])
# """
# function sindyc(X, X_grad=nothing, U=nothing, ts=nothing;
#                 library=Dict(),
#                 use_lasso=false,
#                 hard_threshold=nothing,
#                 quantile_threshold=0.1,
#                 num_terms=nothing,
#                 var_names=nothing)
#                 #TODO: allow custom library functions
#     @warn("DEPRECATED SYNTAX")
#     if X_grad == nothing
#         X_grad = numerical_derivative(X)
#     end
#     if var_names == nothing
#         # TODO: default names for more variables
#         default_names = ["x", "y", "z", "z2", "z3"]
#         var_names = default_names[1:size(X,1)]
#     end
#     if ts == nothing
#         ts = range(0, 1, length=size(X,2))
#         if size(U,1) > 1
#             @warn("Multidimensional control signals may not work")
#         end
#     end
#
#     # Get a model with the augmented data
#     if length(library) == 0
#         error("Must add library terms.")
#     end
#     n, m = size(X)
#     library = convert_string2function(library)
#     X_augmented = calc_augmented_data(X, library)
#     condition_number = cond(X_augmented);
#     if condition_number > 1e5
#         @warn("Very large condition number detected ($condition_number); SINDy algorithm may be unstable")
#     end
#     if !use_lasso
#         # i.e. is a dense model
#         A, B = dmdc(X_augmented, X_grad, U)
#         A = A[1:n, :]
#         if U !== nothing
#             B = B[1:n, :]
#         end
#     else
#         # UPDATE: use my own sequential least squares threshold
#         if U == nothing
#             A = sparse_regression(X_augmented, X_grad,
#                                 hard_threshold=hard_threshold,
#                                 quantile_threshold=quantile_threshold,
#                                 num_terms=num_terms)
#         else
#             Ω = [X_augmented; U]
#             AB = sparse_regression(Ω, X_grad,
#                                 hard_threshold=hard_threshold,
#                                 quantile_threshold=quantile_threshold,
#                                 num_terms=num_terms)
#             A = AB[:, 1:size(X_augmented,1)]
#             B = AB[:, size(X_augmented,1)+1:end]
#         end
#
#         # Do lasso with cross validation; loop over variables
#         # A = ones(n, size(X_augmented,1))
#         # df = build_dataframe(X_augmented, X_grad, library)
#         # regressors = names(df)[n+1:end]
#         # for i in 1:n
#             #TODO: the formula here may need to be updated
#             # this_predictor = names(df)[i] # Derivatives
#             # this_f = this_predictor ~ sum(StatsModels.terms.(regressors))
#             # lasso_model = fit(LassoModel, this_f, df;
#             #                 select=MinCVmse(Kfold(3,2)))
#         # end
#     end
#
#     if U == nothing
#         model = sindycModel(ts, A, zeros(n,1), zeros(1, m), (t)->zeros(1),
#                             library, var_names)
#     else
#         # TODO: make U_func work with multiple channels
#         # U_func(t) = CubicSplineInterpolation(ts, vec(U))(t)
#         U_func = generate_map_rows_function(
#                 (d)->CubicSplineInterpolation(ts,d), U)
#         model = sindycModel(ts, A, B, U, U_func,
#                             library, var_names )
#     end
#
#     return model
# end

# Overloading to work with a sindycModel directly
# function sparse_regression(m::sindycModel, X, X_grad)
#   TODO

"""
Naive implementation of Sparse Identification of Nonlinear DYnamics
with control (SINDYc).

New optimizer-object based interface

There are several nonlinear library terms that can be implemented via passing
a list of strings and arguments via 'library', but a custom function can also be
passed using 'custom_func'. Currently implemented library terms are:
    ["cross_terms", order::Int]
        Here, 'order' is how high of an order to do

model = sindyc(X, X_grad=nothing, U=nothing, ts=nothing;
                library=Dict(),
                optimizer=denseSolver(),
                derivative_function=numerical_derivative,
                var_names=nothing)
"""
function sindyc(X, X_grad=nothing, U=nothing, ts=nothing;
                library=Dict(),
                optimizer=denseSolver(),
                derivative_function=numerical_derivative,
                var_names=nothing)
                #TODO: allow custom library functions
    if X_grad == nothing; X_grad = derivative_function(X); end
    if var_names == nothing
        # TODO: default names for more variables
        default_names = ["x", "y", "z", "x2", "y2", "z2"]
        var_names = default_names[1:size(X,1)]
    end
    if ts == nothing
        ts = range(0, 1, length=size(X,2))
        if size(U,1) > 1
            @warn("Multidimensional control signals may not work")
        end
    end

    # Get a model with the augmented data
    if length(library) == 0
        error("Must add library terms.")
    end
    n, m = size(X)
    library = convert_string2function(library)
    X_augmented = calc_augmented_data(X, library)
    condition_number = cond(X_augmented);
    if condition_number > 1e5
        @warn("Very large condition number detected ($condition_number); SINDy algorithm may be unstable")
    end
    # Use my own sequential least squares threshold objects
    if U == nothing
        A = sparse_regression(optimizer, X_augmented, X_grad)
        model = sindycModel(ts, A, zeros(n,1), zeros(1, m), (t)->zeros(1),
                            library, var_names, optimizer)
    else
        Ω = [X_augmented; U]
        AB = sparse_regression(optimizer, Ω, X_grad)
        A = AB[:, 1:size(X_augmented,1)]
        B = AB[:, size(X_augmented,1)+1:end]
        # TODO: make U_func work with multiple channels
        # U_func(t) = CubicSplineInterpolation(ts, vec(U))(t)
        U_func = generate_map_rows_function(
                (d)->CubicSplineInterpolation(ts,d), U)
        model = sindycModel(ts, A, B, U, U_func,
                            library, var_names, optimizer)
    end

    return model
end

"""
Simple implementation of Sparse Identification of Nonlinear DYnamics
(SINDy).

New optimizer-object based interface

There are several nonlinear library terms that can be implemented via passing
a list of strings and arguments via 'library', but a custom function can also be
passed using 'custom_func'. Currently implemented library terms are:
    ["cross_terms", order::Int]
        Here, 'order' is how high of an order to do

model = sindyc(X, X_grad=nothing, ts=nothing;
                library=Dict(),
                optimizer=denseSolver(),
                derivative_function=numerical_derivative,
                var_names=nothing)
"""
function sindy(X, X_grad=nothing, ts=nothing;
                library=Dict(),
                optimizer=denseSolver(),
                derivative_function=numerical_derivative,
                var_names=nothing)
                #TODO: allow custom library functions
    if X_grad == nothing; X_grad = derivative_function(X); end
    if var_names == nothing
        # TODO: default names for more variables
        default_names = ["x", "y", "z", "x2", "y2", "z2"]
        var_names = default_names[1:size(X,1)]
    end
    if ts == nothing
        ts = range(0, 1, length=size(X,2))
    end

    # Get a model with the augmented data
    if length(library) == 0
        error("Must add library terms.")
    end
    n, m = size(X)
    library = convert_string2function(library)
    X_augmented = calc_augmented_data(X, library)
    condition_number = cond(X_augmented);
    if condition_number > 1e5
        @warn("Very large condition number detected ($condition_number); SINDy algorithm may be unstable")
    end
    # Use my own sequential least squares threshold objects
    A = sparse_regression(optimizer, X_augmented, X_grad)
    model = sindyModel(ts, A, library, var_names, optimizer)

    return model
end

###
### Helper functions for making many similar models
###

"""
Makes a model like the given template, with a different training set
"""
function sindy_retrain(m::sindyModel, X::Matrix, X_grad::Matrix)
    sindy(X, X_grad, m.ts;
        library=m.library,
        optimizer=m.optimizer,
        var_names=m.variable_names)
end

"""
Makes a model like the given template, with a different training set
"""
function sindyc_retrain(m::sindycModel, X::Matrix, X_grad::Matrix)
    if size(X,2) !== size(m.U,2)
        error("DataError: Must pass new control signal if data is different")
    end
    sindyc(X, X_grad, m.U, m.ts;
        library=m.library,
        optimizer=m.optimizer,
        var_names=m.variable_names)
end

function sindyc_retrain(m::sindycModel, X::Matrix, X_grad::Matrix, U::VecOrMat)
    if size(X,2) !== size(m.U,2)
        error("DataError: Control signal ($(size(m.U,2))) should be same length as data ($(size(X,2)))")
    end
    sindyc(X, X_grad, m.U, m.ts;
        library=m.library,
        optimizer=m.optimizer,
        var_names=m.variable_names)
end

#####
##### Helper functions for pretty printing
#####
"""
print_equations(model::sindycModel;
                     var_names=nothing, tol=1e-4, digits=1)

Prints equations of a SINDy model using saved or passed variable names
    Note: if 0.0 is displayed, this means the coefficient was nonzero to the
    tolerance specified in 'tol' but below the rounding threshold
"""
function print_equations(model::DynamicalSystemModel;
                         var_names=nothing, tol=1e-4, digits=1)
    n, m = size(model.A)
    if var_names == nothing
        # default_names = ["x"; "y"; "z"]
        # var_names = default_names[1:n]
        var_names = model.variable_names;
    end
    term_names = build_term_names(model, var_names)
    for i_var in 1:n
        this_var = var_names[i_var]
        print("d$this_var/dt = ")
        first_term = true
        for i_term in 1:m
            val = model.A[i_var, i_term]
            if abs(val) > tol
                this_term = term_names[i_term]
                val = round(val, digits=digits)
                if first_term
                    print("$val$this_term")
                    first_term = false
                else
                    print(" + $val$this_term")
                end
            end
        end
        println("\b\b")
    end
    # println()
end

function get_nonzero_terms(model::DynamicalSystemModel; linear_indices=false)
    # A = reshape(model.A, length(model.A))
    if linear_indices
        lin = LinearIndices(model.A)
        return lin[findall(abs.(model.A).>0.0)]
    else
        return findall(abs.(model.A).>0.0)
    end
end

function get_nonzero_term_names(model::DynamicalSystemModel)
    rhs_names = build_term_names(model)
    term_ind = get_nonzero_terms(model)
    lhs_names = model.variable_names
    final_term_names = ones(String, size(term_ind))
    for (i, mat_ind) in enumerate(term_ind)
        lhs_var, rhs_var = Tuple(mat_ind)
        lhs = lhs_names[lhs_var]
        rhs = rhs_names[rhs_var]
        final_term_names[i] = "d$(lhs)dt_$rhs"
    end
    return final_term_names
end

#####
##### Helper functions for nonlinear terms
#####
"""
Calculates higher powers of individual rows of X, but not cross terms
"""
function calc_power_terms(X, order)
    n = size(X, 1)
    all_powers = nchoosek(order, n)
    popfirst!(all_powers) # i.e. remove the identity
    power_terms = nothing

    for p_list in all_powers
        this_power_term = zeros(size(X))
        for (i,p) in enumerate(p_list)
            this_power_term[i,:] = X[i,:] .^ p
        end
        if power_terms == nothing
            power_terms = copy(this_power_term)
        else
            power_terms = vcat(power_terms, this_power_term)
        end
    end

    return power_terms
end

"""
Calculates higher powers of individual rows of X, AND cross terms
    NOT cumulative, but a vector can be passed for 'order'
"""
function calc_cross_terms(X, order)
    if order isa Vector
        # Call recursively
        cross_terms = calc_cross_terms(X, order[1])
        for i in 2:length(order)
            this_cross_terms = calc_cross_terms(X, order[i])
            cross_terms = vcat(cross_terms, this_cross_terms)
        end
        return cross_terms
    end
    if ndims(X) == 1
        X = reshape(X, (length(X), 1))
    end
    n, m = size(X)
    # TODO: This is NOT combinations, right?
    all_powers = calc_permutations(n, order)
    cross_terms = nothing

    for p_list in all_powers
        this_cross_term = ones(eltype(X), 1, m)
        for p in p_list
            # Keep the X slice 2 dimensional
            this_cross_term .*= X[p:p,:]
        end
        if cross_terms == nothing
            cross_terms = copy(this_cross_term)
        else
            cross_terms = vcat(cross_terms, this_cross_term)
        end
    end

    return cross_terms
end

"""
Returns a constant; has the right function signature
"""
function calc_constant(X, tmp)
    return ones(eltype(X), 1, size(X,2))
end

FUNCTION_DICT = Dict("cross_terms"=>calc_cross_terms,
                     "power_terms"=>calc_power_terms,
                     "constant"=>calc_constant)

#####
##### Helper functions for variable names
#####
 """Convert library strings to functions"""
 function convert_string2function(library)
     func_library = OrderedDict{Function, Any}()
     for f in keys(library)
         if f isa String
             f_key = FUNCTION_DICT[f]
         else
             f_key = f
         end
         func_library[f_key] = library[f] # Key is now the function
     end
     return func_library
 end

 function build_term_names(library::Union{Dict, OrderedDict}, var_names)
     term_names = copy(var_names)
     # Now we pass through library functions
     library_vec = collect(library)
     for (f, args) in library_vec
         these_term_names = vec(f(var_names, args))
         term_names = vcat(term_names, these_term_names)
     end

     return term_names
 end
 build_term_names(model::DynamicalSystemModel, var_names) =
     build_term_names(model.library, var_names)
 build_term_names(model::DynamicalSystemModel) =
     build_term_names(model.library, model.variable_names)

 #####
 ##### Helper functions for calculation
 #####
"""Uses saved library functions to calculate the input space"""
function calc_augmented_data(X, library)
    # Apply functions and save
    X_augmented = copy(X)
    for (f, args) in pairs(library)
        X_augmented = vcat(X_augmented, f(X, args))
    end

    n, m = size(X_augmented)
    if m == 1
        # ode solvers want a vector
        X_augmented = reshape(X_augmented, (n))
    end

    return X_augmented
end

# Flux version!
#   Vertical concatenation doesn't easily work
function calc_augmented_data(X::TrackedArray, library)
    # Apply functions and save
    X_augmented = copy(X)
    for (f, args) in pairs(library)
        X_augmented = vcat(X_augmented,
                Flux.data.(f(X, args)))
    end

    n, m = size(X_augmented)
    if m == 1
        # ode solvers want a vector
        X_augmented = reshape(X_augmented, (n))
    end

    return X_augmented
end

"""
Build dataframe from the augmented data
"""
function build_dataframe(X_augmented, dXdt, library)
    var_names = ["x", "y", "z"]
    df = DataFrame()
    # First add derivatives
    for i in 1:size(dXdt,1)
        name = var_names[i]
        name_sym = Symbol("d$name/dt")
        df[!, name_sym] = vec(dXdt[i,:])
    end
    # Now other terms
    term_names = build_term_names(library, var_names)
    for i in 1:length(term_names)
        name = term_names[i]
        if length(name) == 0
            name = "const"
        end
        name_sym = Symbol(name)
        df[!, name_sym] = vec(X_augmented[i,:])
    end

    return df
end

export calc_augmented_data, convert_string2function, calc_constant,
    FUNCTION_DICT, calc_cross_terms, calc_power_terms,
    build_dataframe, sindyc, sindy, print_equations,
    get_nonzero_terms, get_nonzero_term_names, build_term_names,
    sindyc_retrain, sindy_retrain


text:
module EMLcompressor
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("compressor/centrifugal_compressor.jl")
end


text:
module QuadratureRules

using LinearAlgebra

export gauleg

include("gauleg.jl")

end


text:
module TestDoctest

using Documenter
using Test
using ContextManagers

function test()
    if VERSION < v"1.6"
        @test_broken false
        return
    end
    doctest(ContextManagers; manual = true)
end

end  # module


text:
# Tracing charged particle in a static analytic dipole magnetic field.
# This example shows how to trace protons of a certain energy in a analytic
# Earth-like magnetic dipole field. Pyplot is used for interactively changing
# views in 3D.
#
# Hongyang Zhou, hyzhou@umich.edu

using TestParticle
using TestParticle: getB_dipole, getE_dipole, sph2cart, Rₑ
using OrdinaryDiffEq
using PyPlot

Ek = 5e7 # [eV]

m = TestParticle.mᵢ
q = TestParticle.qᵢ
c = TestParticle.c
Rₑ = TestParticle.Rₑ

# initial velocity, [m/s]
v₀ = sph2cart(c*sqrt(1-1/(1+Ek*q/(m*c^2))^2), 0.0, π/4)
# initial position, [m]
r₀ = sph2cart(2.5*Rₑ, 0.0, π/2)
stateinit = [r₀..., v₀...]
# obtain field
param = prepare(getE_dipole, getB_dipole)
tspan = (0.0, 2.0)

prob = ODEProblem(trace_analytic!, stateinit, tspan, param)

sol = solve(prob, Tsit5(); save_idxs=[1,2,3])

x = getindex.(sol.u,1) / Rₑ
y = getindex.(sol.u,2) / Rₑ
z = getindex.(sol.u,3) / Rₑ

## Visualization

using3D()
fig = plt.figure()
ax = fig.gca(projection="3d")

plot(x,y,z)

for ϕ in range(0, stop=2*π, length=10)
   plot(fieldline(ϕ)..., color="r", alpha=0.3)
end


text:
"""
Lightweight 3D attitude representation module.

Defines the abstract type `Abstract3DRotation`, and implements four concrete subtypes:
- Unit quaternion (`RQuat`)
- Rotation matrix (`RMatrix`)
- Axis-angle (`RAxAng`)
- Euler angles (`REuler`)

All of them support the essential attitude operations used in this package:
- Equality (`==`) and approximate equality (`≈`)
- Inversion (`r'`)
- Composition (`r1 ∘ r2`)
- Vector coordinate transformation (`r * v == r(v)`)
- Conversion to and from other representations (`r1 = R1(r2)`)

`RQuat` operates as a primary representation. It implements all of the
operations natively, and therefore it can be used for promotion whenever a
specific operation cannot be directly performed on its original operand(s).

This allows for a common interface; any `Abstract3DRotation` subtype need only implement
conversions to and from `RQuat` to be fully compatible with all other subtypes
and operations. However, for efficiency reasons, some direct methods are
provided for the implemented subtypes.
"""
module Attitude

using StaticArrays
using StructArrays
using LinearAlgebra
using Flight.Quaternions
using Flight.Plotting
#using ..Quaternions: UnitQuat, Quat #also works (but relies on folder hierarchy)

export Abstract3DRotation, RQuat, RAxAng, REuler, RMatrix, Rx, Ry, Rz, dt

const half_π = π/2

"""
Compute the skew-symmetric matrix corresponding to a 3-element vector.
"""
function skew(v::AbstractVector)

    @SMatrix [0 -v[3] v[2];
            v[3] 0 -v[1];
            -v[2] v[1] 0]

end

######################### Abstract3DRotation ###########################

"""
Generic 3D rotation descriptor
"""
abstract type Abstract3DRotation end

"Function call notation for coordinate transformations"
(r::Abstract3DRotation)(v::AbstractVector{<:Real}) = r * v

############################# RQuat ###############################

"Unit quaternion representation"
struct RQuat <: Abstract3DRotation
    _u::UnitQuat #no normalization is performed for an UnitQuat input
end

#generic AbstractVectors are normalized by default
function RQuat(v::AbstractVector{<:Real}; normalization::Bool = true)
    RQuat(UnitQuat(SVector{4,Float64}(v); normalization))
end
RQuat() = RQuat(UnitQuat(1.0))
RQuat(r::Abstract3DRotation) = convert(RQuat, r)

Base.getindex(r::RQuat, i) = getindex(r._u, i)
Base.length(r::RQuat) = length(r._u)
Base.size(r::RQuat) = size(r._u)
# Base.iterate(r::RQuat, state = 1) = iterate(getfield(getfield(r._u, :_q), :_sv)[:], state)

# direct iteration unsupported because it allocates. doing x.=r[:] instead of
#x.=r gets the underlying SVector first by forwarding getindex, then iterates on
#the StateVector itself. this avoids allocations.

Base.show(io::IO, r::RQuat) = print(io, "$(typeof(r))($(r[:]))")

#Strict equality (accounts for the unit quaternions' double cover of SO(3))
Base.:(==)(r1::RQuat, r2::RQuat) = (r1._u == r2._u || r1._u == -r2._u)
# Approximate equality (accounts for the unit quaternions' double cover of SO(3))
Base.:(≈)(r1::RQuat, r2::RQuat) = (r1._u ≈ r2._u || r1._u ≈ -r2._u)
# Composition
Base.:∘(r1::RQuat, r2::RQuat) = RQuat(r1._u * r2._u)
# Inversion
Base.adjoint(r::RQuat) = RQuat(r._u', normalization = false)

# Coordinate transformation (function call notation also supported)
function Base.:*(r_ab::RQuat, v_b_in::AbstractVector{<:Real})::SVector{3,Float64}
    v_b = SVector{3, Float64}(v_b_in) #casting to SVector yields a 3x speed gain
    q = r_ab._u; q_re = q.real; q_im = q.imag
    v_a = v_b + 2q_im × (q_re * v_b + q_im × v_b)
    return v_a
end

LinearAlgebra.norm(r::RQuat) = norm(r._u)
LinearAlgebra.normalize(r::RQuat) = RQuat(normalize(r._u))

"""
    dt(r_ab::RQuat, ω_ab_b::AbstractVector{<:Real})

Unit quaternion time derivative

If `r_ab` represents the rotation from axes εa to axes εb, and `ω_ab_b` is the
angular velocity of εb with respect to εa projected in εb, then `ṙ_ab ==
dt(r_ab, ω_ab_b)`
"""
dt(r_ab::RQuat, ω_ab_b::AbstractVector{<:Real}) = 0.5 * (r_ab._u * FreeQuat(imag=ω_ab_b))

##### RQuat fallbacks #####

#require each Abstract3DRotation subtype to implement conversions to and from RQuat
Base.convert(::Type{RQuat}, r::R) where {R<:Abstract3DRotation} = error("Implement $R to RQuat conversion")
Base.convert(::Type{R}, r::RQuat) where {R<:Abstract3DRotation} = error("Implement RQuat to $R conversion")

#this enables a two-step conversion via RQuat, but requires defining trivial
#conversions to avoid infinite recursion
Base.convert(::Type{R}, r::Abstract3DRotation) where {R<:Abstract3DRotation} = convert(R, RQuat(r))
Base.convert(::Type{R}, r::R) where {R<:Abstract3DRotation} = r
Base.convert(::Type{RQuat}, r::RQuat) = r

#unless the representation defines its own methods, fall back to RQuat;
#promote rules don't work here, because promoting two inputs of the same subtype
#leaves them unchanged, and here we need them converted to RQuat
Base.adjoint(r::R) where {R<:Abstract3DRotation} = convert(R, RQuat(r)')
Base.:*(r::Abstract3DRotation, v::AbstractVector{<:Real}) = RQuat(r) * v
Base.:(≈)(r1::Abstract3DRotation, r2::Abstract3DRotation) = RQuat(r1) ≈ RQuat(r2)
Base.:∘(r1::Abstract3DRotation, r2::Abstract3DRotation) = RQuat(r1) ∘ RQuat(r2)

#absolute equality is not defined in general, because comparing two different
#subtypes requires promotion to RQuat, whose outcome is affected by floating
#point precision

######################### RMatrix #######################

"Rotation matrix representation"
struct RMatrix <: Abstract3DRotation
    _mat::SMatrix{3, 3, Float64, 9}
    function RMatrix(input::AbstractArray{<:Real, 2}; normalization::Bool = true)
        return normalization ? new(qr(input).Q) : new(input)
    end
end

#normalize by default
RMatrix() = RMatrix(SMatrix{3,3,Float64}(I), normalization = false)
RMatrix(r::Abstract3DRotation) = convert(RMatrix, r)

Base.size(r::RMatrix) = size(getfield(r,:_mat))
Base.getindex(r::RMatrix, i...) = getindex(getfield(r,:_mat), i...)

Base.convert(::Type{RMatrix}, r::RMatrix) = r

function Base.convert(::Type{RMatrix}, r::RQuat)

    q = normalize(r._u) #cheap
    q_sq = q[:] .* q[:]
    dq12 = 2*q[1]*q[2]; dq13 = 2*q[1]*q[3]; dq14 = 2*q[1]*q[4];
    dq23 = 2*q[2]*q[3]; dq24 = 2*q[2]*q[4]; dq34 = 2*q[3]*q[4];

    M = @SMatrix [
    1 - 2*(q_sq[3]+q_sq[4])  dq23 - dq14                dq24 + dq13;
    dq23 + dq14              1 - 2*(q_sq[2]+q_sq[4])    dq34 - dq12;
    dq24 - dq13              dq34 + dq12                1 - 2*(q_sq[2]+q_sq[3])
    ]

    RMatrix(M, normalization = false)

end

function Base.convert(::Type{RQuat}, r::RMatrix)
    R = r._mat
    tr_R = tr(R)
    i_max = findmax(vcat(tr_R, diag(R)))[2]

    #these 4-element vectors are all proportional to the the desired unit
    #quaternion's components. to minimize numerical errors, we choose the one
    #with the largest norm (indicated by i_max), and construct the quaternion by
    #normalizing it
    if i_max == 1
        v = SVector{4,Float64}(
            1 + tr_R,
            R[3,2] - R[2,3],
            R[1,3] - R[3,1],
            R[2,1] - R[1,2]
        )
    elseif i_max == 2
        v = SVector{4,Float64}(
            R[3,2] - R[2,3],
            1 + 2R[1,1] - tr_R,
            R[1,2] + R[2,1],
            R[3,1] + R[1,3]
        )
    elseif i_max == 3
        v = SVector{4,Float64}(
            R[1,3] - R[3,1],
            R[1,2] + R[2,1],
            1 + 2R[2,2] - tr_R,
            R[2,3] + R[3,2]
        )
    else #i_max == 4
        v = SVector{4, Float64}(
            R[2,1] - R[1,2],
            R[3,1] + R[1,3],
            R[2,3] + R[3,2],
            1 + 2R[3,3] - tr_R
        )
    end

    return RQuat(v, normalization = true)

end

#RMatrix also implements all operators natively
Base.:(==)(r1::RMatrix, r2::RMatrix) = (r1._mat == r2._mat)
Base.:(≈)(r1::RMatrix, r2::RMatrix) = r1._mat ≈ r2._mat

# Coordinate transformation
Base.:*(r::RMatrix, v::AbstractVector{<:Real}) = r._mat * v

# Composition
Base.:∘(r1::RMatrix, r2::RMatrix) = RMatrix(r1._mat * r2._mat, normalization = false)
Base.:*(r1::RMatrix, r2::RMatrix) = r1 ∘ r2 #allow overload of * for composition

# Inversion
Base.adjoint(r::RMatrix) = RMatrix(r._mat', normalization = false)

#let the constructor handle it
LinearAlgebra.normalize(r::RMatrix) = RMatrix(r._mat, normalization = true)

LinearAlgebra.det(r::RMatrix) = det(r._mat)

"""
    dt(r_ab::RMatrix, ω_ab_b::AbstractVector{<:Real})

Rotation matrix time derivative

If `R_ab` represents the rotation from axes εa to axes εb, and `ω_ab_b` is the
angular velocity of εb with respect to εa projected in εb, then `Ṙ_ab ==
dt(R_ab, ω_ab_b)`
"""
function dt(R_ab::RMatrix, ω_ab_b::AbstractVector{<:Real})
    Ω_ab_b = skew(ω_ab_b)
    return R_ab._mat * Ω_ab_b
end

############################# RAxAng ###############################

"Axis-angle representation"
struct RAxAng <: Abstract3DRotation
    axis::SVector{3, Float64}
    angle::Float64
    function RAxAng(axis::AbstractVector{<:Real}, angle::Real; normalization::Bool = true)
        axis = SVector{3,Float64}(axis) #normalization will be faster for an SVector
        return normalization ? new(normalize(axis), angle) : new(axis, angle)
    end
end

RAxAng(r::Abstract3DRotation) = convert(RAxAng, r)
RAxAng(::Nothing, ::Real; normalization::Bool = false) = RAxAng()
RAxAng(input::Tuple{Union{Nothing, AbstractVector{<:Real}}, Real}) = RAxAng(input...)

#shorthand constructors for rotations around X, Y, Z
Rz(ψ::Real = 0.0) = RAxAng(SVector{3}(0.0, 0.0, 1.0), ψ, normalization = false)
Ry(θ::Real = 0.0) = RAxAng(SVector{3}(0.0, 1.0, 0.0), θ, normalization = false)
Rx(φ::Real = 0.0) = RAxAng(SVector{3}(1.0, 0.0, 0.0), φ, normalization = false)

RAxAng() = Rx(0)

#conversions to and from RQuat
function Base.convert(::Type{RAxAng}, r::RQuat)
    q_re = r._u.real
    q_im = r._u.imag
    norm_im = norm(q_im)
    μ = 2atan(norm_im, q_re)
    u = (norm_im > 0 ? q_im / norm_im : nothing)
    return RAxAng(u, μ, normalization = false)
end

function Base.convert(::Type{RQuat}, r::RAxAng)
    u = r.axis
    μ = r.angle
    RQuat(UnitQuat(real = cos(0.5μ), imag = u*sin(0.5μ), normalization = false))
end

# Inversion (the only native operation we can easily define for axis-angle)
Base.adjoint(r::RAxAng) = RAxAng(r.axis, -r.angle)

############################# REuler #############################

"Euler angle representation (convention is ZYX)"
struct REuler <: Abstract3DRotation
    ψ::Float64 #heading
    θ::Float64 #inclination
    φ::Float64 #bank
    function REuler(ψ, θ, φ)
        @assert abs(ψ<=π) "Heading must be within [-π, π]"
        @assert abs(θ<=half_π) "Inclination must be within [-π/2, π/2]"
        @assert abs(φ<=π) "Bank must be within [-π, π]"
        new(ψ, θ, φ)
    end
end

REuler(r::Abstract3DRotation) = convert(REuler, r)
REuler(input::Tuple{Real, Real, Real}) = REuler(input...)
REuler(; ψ = 0, θ = 0, φ = 0) = REuler(ψ, θ, φ)

#conversions to and from RQuat
function Base.convert(::Type{REuler}, r::RQuat)
        q = r._u
        q_sq = q[:] .* q[:]

        ψ = atan( 2*(q[1]*q[4] + q[2]*q[3]), 1 - 2*(q_sq[3] + q_sq[4]))
        θ = asin( clamp(2*(q[1]*q[3] - q[2]*q[4]), -1, 1) )
        φ = atan( 2*(q[1]*q[2] + q[3]*q[4]), 1 - 2*(q_sq[2] + q_sq[3]))

        return REuler(ψ, θ, φ)
end

function Base.convert(::Type{RQuat}, r::REuler)
    Rz(r.ψ) ∘ Ry(r.θ) ∘ Rx(r.φ)
end

########################## Plotting #################################

#unless a more specialized method is defined, a TimeHistory{Rotations} is
#converted to REuler for plotting
@recipe function plot_rotation(th::TimeHistory{<:AbstractVector{<:Abstract3DRotation}};
                                rot_ref = "", rot_target = "")

    v_euler = Vector{REuler}(undef, length(th.data))
    for i in 1:length(v_euler)
        v_euler[i] = REuler(th.data[i])
    end
    sa = StructArray(v_euler)
    data = hcat(sa.ψ, sa.θ, sa.φ)./π #plot as π factors

    label --> ["Heading" "Inclination" "Bank"]
    yguide --> hcat(L"$\psi_{%$rot_ref %$rot_target} \ (\pi \ rad)$",
                    L"$\theta_{%$rot_ref %$rot_target} \ (\pi \ rad)$",
                    L"$\phi_{%$rot_ref %$rot_target} \ (\pi \ rad)$")
    th_split --> :h #custom TimeHistory attribute

    return TimeHistory(th.t, data)

end

end


text:
"""The default search location for Pseudopotential data files"""
DFTK_DATADIR = joinpath(dirname(pathof(DFTK)), "..", "data")


text:
# Copyright (c) 2022 Dongdong Kong. All rights reserved.
# This work is licensed under the terms of the MIT license.
# For a copy, see <https://opensource.org/licenses/MIT>.

module CMIP

using DataFrames
# import nctools: nc_calendar | not work
using nctools
import nctools.Ipaper: dates_miss, dates_nmiss

include("CMIPFiles_info.jl")

# include("unit_convert.jl")
include("heat_index.jl")

export Tem_F2C, Tem_C2F, heat_index

end


text:
using EMSimJulia
using Test

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


text:
const log2π = log(2*π)
const sqrt2π = √(2*π)
const sqrtπ = √π
# Euler–Mascheroni constant
const γ = 0.577215665

# Ramanujan's approximation to the Γ function:
Γ(x::Real) = sqrtπ*((x-1)/exp(1))^(x-1) * (8*(x-1)^3 + 4*(x-1)^2 + (x-1) + 1/30)^(1/6)


text:
using RDatasets, DataFrames, GLM, Plots, Measures; pyplot()

df = dataset("datasets", "anscombe")

model1 = lm(@formula(Y1 ~ X1), df)
model2 = lm(@formula(Y2 ~ X2), df)
model3 = lm(@formula(Y3 ~ X3), df)
model4 = lm(@formula(Y4 ~ X4), df)

println("Model 1. Coefficients: ", coef(model1),"\t R squared: ",r2(model1))
println("Model 2. Coefficients: ", coef(model2),"\t R squared: ",r2(model2))
println("Model 3. Coefficients: ", coef(model3),"\t R squared: ",r2(model3))
println("Model 4. Coefficients: ", coef(model4),"\t R squared: ",r2(model4))

yHat(model, X) = coef(model)' * [ 1 , X ]
xlims = [0, 20]

p1 = scatter(df.X1, df.Y1, c=:blue, msw=0, ms=8)
p1 = plot!(xlims, [yHat(model1, x) for x in xlims], c=:red, xlims=(xlims))

p2 = scatter(df.X2, df.Y2, c=:blue, msw=0, ms=8)
p2 = plot!(xlims, [yHat(model2, x) for x in xlims], c=:red, xlims=(xlims))

p3 = scatter(df.X3, df.Y3, c=:blue, msw=0, ms=8)
p3 = plot!(xlims, [yHat(model3, x) for x in xlims], c=:red, xlims=(xlims))

p4 = scatter(df.X4, df.Y4, c=:blue, msw=0, ms=8)
p4 = plot!(xlims, [yHat(model4, x) for x in xlims], c=:red, msw=0, xlims=(xlims))

plot(p1, p2, p3, p4, layout = (2,2), xlims=(0,20), ylims=(0,14),
	legend=:none, xlabel = "x", ylabel = "y",
    	size=(1200, 800), margin = 10mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module Foo
    #include("using_json.jl")
    include("parse.jl")
end


text:
#####
##### Navier-Stokes and tracer advection equations
#####

""" Store previous value of the source term and calculate current source term. """
function calculate_interior_tendency_contributions!(G, arch, grid, advection, coriolis, buoyancy, surface_waves, closure,
                                                    U, C, pHY′, K, F, clock)

    workgroup, worksize = work_layout(grid, :xyz)

    calculate_Gu_kernel! = calculate_Gu!(device(arch), workgroup, worksize)
    calculate_Gv_kernel! = calculate_Gv!(device(arch), workgroup, worksize)
    calculate_Gw_kernel! = calculate_Gw!(device(arch), workgroup, worksize)
    calculate_Gc_kernel! = calculate_Gc!(device(arch), workgroup, worksize)

    barrier = Event(device(arch))

    Gu_event = calculate_Gu_kernel!(G.u, grid, advection, coriolis, surface_waves, closure,
                                    U, C, K, F, pHY′, clock, dependencies=barrier)

    Gv_event = calculate_Gv_kernel!(G.v, grid, advection, coriolis, surface_waves, closure,
                                    U, C, K, F, pHY′, clock, dependencies=barrier)

    Gw_event = calculate_Gw_kernel!(G.w, grid, advection, coriolis, surface_waves, closure,
                                    U, C, K, F, clock, dependencies=barrier)

    events = [Gu_event, Gv_event, Gw_event]

    for tracer_index in 1:length(C)
        @inbounds Gc = G[tracer_index+3]
        @inbounds Fc = F[tracer_index+3]
        @inbounds  c = C[tracer_index]

        Gc_event = calculate_Gc_kernel!(Gc, grid, c, Val(tracer_index), advection, closure,
                                        buoyancy, U, C, K, Fc, clock, dependencies=barrier)

        push!(events, Gc_event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

#####
##### Tendency calculators for u, v, w-velocity
#####

""" Calculate the right-hand-side of the u-velocity equation. """
@kernel function calculate_Gu!(Gu, grid, advection, coriolis, surface_waves, closure, U, C, K, F, pHY′, clock)
    i, j, k = @index(Global, NTuple)

    @inbounds Gu[i, j, k] = u_velocity_tendency(i, j, k, grid, advection, coriolis, surface_waves,
                                                closure, U, C, K, F, pHY′, clock)
end

""" Calculate the right-hand-side of the v-velocity equation. """
@kernel function calculate_Gv!(Gv, grid, advection, coriolis, surface_waves, closure, U, C, K, F, pHY′, clock)
    i, j, k = @index(Global, NTuple)

    @inbounds Gv[i, j, k] = v_velocity_tendency(i, j, k, grid, advection, coriolis, surface_waves,
                                                closure, U, C, K, F, pHY′, clock)
end

""" Calculate the right-hand-side of the w-velocity equation. """
@kernel function calculate_Gw!(Gw, grid, advection, coriolis, surface_waves, closure, U, C, K, F, clock)
    i, j, k = @index(Global, NTuple)

    @inbounds Gw[i, j, k] = w_velocity_tendency(i, j, k, grid, advection, coriolis, surface_waves,
                                                closure, U, C, K, F, clock)
end

#####
##### Tracer(s)
#####

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
@kernel function calculate_Gc!(Gc, grid, c, tracer_index, advection, closure, buoyancy, U, C, K, Fc, clock)
    i, j, k = @index(Global, NTuple)

    @inbounds Gc[i, j, k] = tracer_tendency(i, j, k, grid, c, tracer_index,
                                            advection, closure, buoyancy, U, C, K, Fc, clock)
end

#####
##### Boundary contributions to tendencies due to user-prescribed fluxes
#####

""" Apply boundary conditions by adding flux divergences to the right-hand-side. """
function calculate_boundary_tendency_contributions!(Gⁿ, arch, U, C, clock, state)

    barrier = Event(device(arch))

    events = []

    # Velocity fields
    for i in 1:3
        x_bcs_event = apply_x_bcs!(Gⁿ[i], U[i], arch, barrier, clock, state)
        y_bcs_event = apply_y_bcs!(Gⁿ[i], U[i], arch, barrier, clock, state)
        z_bcs_event = apply_z_bcs!(Gⁿ[i], U[i], arch, barrier, clock, state)

        push!(events, x_bcs_event, y_bcs_event, z_bcs_event)
    end

    # Tracer fields
    for i in 4:length(Gⁿ)
        x_bcs_event = apply_x_bcs!(Gⁿ[i], C[i-3], arch, barrier, clock, state)
        y_bcs_event = apply_y_bcs!(Gⁿ[i], C[i-3], arch, barrier, clock, state)
        z_bcs_event = apply_z_bcs!(Gⁿ[i], C[i-3], arch, barrier, clock, state)

        push!(events, x_bcs_event, y_bcs_event, z_bcs_event)
    end

    events = filter(e -> typeof(e) <: Event, events)

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

#####
##### Vertical integrals
#####

"""
Update the hydrostatic pressure perturbation pHY′. This is done by integrating
the `buoyancy_perturbation` downwards:

    `pHY′ = ∫ buoyancy_perturbation dz` from `z=0` down to `z=-Lz`
"""
@kernel function update_hydrostatic_pressure!(pHY′, grid, buoyancy, C)
    i, j = @index(Global, NTuple)

    @inbounds pHY′[i, j, grid.Nz] = - ℑzᵃᵃᶠ(i, j, grid.Nz+1, grid, buoyancy_perturbation, buoyancy, C) * ΔzF(i, j, grid.Nz+1, grid)

    @unroll for k in grid.Nz-1 : -1 : 1
        @inbounds pHY′[i, j, k] =
            pHY′[i, j, k+1] - ℑzᵃᵃᶠ(i, j, k+1, grid, buoyancy_perturbation, buoyancy, C) * ΔzF(i, j, k+1, grid)
    end
end

#####
##### Source term storage
#####

""" Store source terms for `u`, `v`, and `w`. """
@kernel function store_velocity_tendencies!(G⁻, grid::AbstractGrid{FT}, G⁰) where FT
    i, j, k = @index(Global, NTuple)
    @inbounds G⁻.u[i, j, k] = G⁰.u[i, j, k]
    @inbounds G⁻.v[i, j, k] = G⁰.v[i, j, k]
    @inbounds G⁻.w[i, j, k] = G⁰.w[i, j, k]
end

""" Store previous source terms for a tracer before updating them. """
@kernel function store_tracer_tendency!(Gc⁻, grid::AbstractGrid{FT}, Gc⁰) where FT
    i, j, k = @index(Global, NTuple)
    @inbounds Gc⁻[i, j, k] = Gc⁰[i, j, k]
end

""" Store previous source terms before updating them. """
function store_tendencies!(G⁻, arch, grid, G⁰)

    barrier = Event(device(arch))

    workgroup, worksize = work_layout(grid, :xyz)

    store_velocity_tendencies_kernel! = store_velocity_tendencies!(device(arch), workgroup, worksize)
    store_tracer_tendency_kernel! = store_tracer_tendency!(device(arch), workgroup, worksize)

    velocities_event = store_velocity_tendencies_kernel!(G⁻, grid, G⁰, dependencies=barrier)

    events = [velocities_event]

    # Tracer fields
    for i in 4:length(G⁻)
        @inbounds Gc⁻ = G⁻[i]
        @inbounds Gc⁰ = G⁰[i]
        tracer_event = store_tracer_tendency_kernel!(Gc⁻, grid, Gc⁰, dependencies=barrier)
        push!(events, tracer_event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end


text:
module TestTasks

# using Revise
using Test
using MLJ
using Random
using Random: MersenneTwister

X=(x=10:10:44, y=1:4, z=collect("abcd"))
task = SupervisedTask(data=X, target=:y, is_probabilistic=true)

@testset "Shuffling" begin
    task0=deepcopy(task)
    rng = MersenneTwister(1234)
    shuffle!(rng, task0)
    @test task.X != task0.X
    @test task.y != task0.y
    @test MLJ.selectrows(task.X, task0.y) == task0.X

    Random.seed!(1234)
    rng = MersenneTwister(1234)
    task1=shuffle(task)
    shuffle!(rng, task1)
    @test task.X != task1.X
    @test task.y != task1.y
    @test MLJ.selectrows(task.X, task1.y) == task1.X

    Random.seed!(1234)
    task0_a = deepcopy(task)
    shuffle!(task0_a)
    @test task0.X == task0_a.X
    @test task0.y == task0_a.y
    @test MLJ.selectrows(task.X, task0_a.y) == task0_a.X

    Random.seed!(1234)
    rng = MersenneTwister(1234)
    task1_a = shuffle(rng, shuffle(task))
    @test task1.X == task1_a.X
    @test task1.y == task1_a.y
    @test MLJ.selectrows(task.X, task1_a.y) == task1_a.X
end

@testset "Indexing" begin
    task2 = task[:]
    @test count(fieldnames(typeof(task))) do fld
        getproperty(task2, fld) != getproperty(task, fld)
    end == 0
    @test task[2:3].X.z == ['b', 'c']
    @test task[2:3].y == [2, 3]
end

@testset "Constructors" begin
    df = (x=10:10:44, y=1:4, z=collect("abcd"), w=[1.0, 3.0, missing])
    types = Dict(:x => Continuous, :z => Multiclass, :w => Count)

    # following 3 lines were commented because of interfering
    # depreciation warning and replaced by the subsequent two lines:
    # task = @test_logs((:warn, r"Missing values encountered"),
    task =supervised(data=df, types=types, target=:y,
               ignore=:y, is_probabilistic=false)

    @test scitype_union(task.X.x) <: Continuous
    @test scitype_union(task.X.w) === Union{Count, Missing}
    @test scitype_union(task.y) <: Count
    @test_logs((:info, r"\nis_probabilistic = true"),
               supervised(task.X, task.y, is_probabilistic=true))
end

end # module
true


text:
module TrussIO

using ..TrussTopOptProblems: TrussFEACrossSec, TrussFEAMaterial

using StaticArrays
export  parse_truss_json,
        parse_support_load_json

include("parse_json.jl")

end