text:
using Test
using StaticArrays
using Proj4
import PROJ_jll

@testset "Error handling" begin
    @test Proj4.proj_errno_string(0) == nothing
    @test Proj4.proj_errno_string(1) == "Operation not permitted"
    @test Proj4.proj_errno_string(2) == "No such file or directory"

    # this throws an internal error that is handled by our log_func
    @test_throws Proj4.PROJError Proj4.proj_create("+proj=bobbyjoe")
    # ensure we have reset the error state
    @test Proj4.proj_context_errno() == 0

end

@testset "Database" begin
    crs = Proj4.proj_create_from_database("EPSG", "4326", Proj4.PJ_CATEGORY_CRS, false)
    Proj4.proj_errno(crs)
    Proj4.proj_get_id_code(crs, 0)
end

@testset "Create" begin
    pj_latlon = Proj4.proj_create("EPSG:4326")

    esriwkt = """GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]"""
    pj2 = Proj4.proj_create(esriwkt)

    @test startswith(Proj4.proj_as_wkt(pj2, Proj4.PJ_WKT1_ESRI, C_NULL), "GEOGCS")
    @test startswith(Proj4.proj_as_wkt(pj2, Proj4.PJ_WKT2_2018, C_NULL), "GEOGCRS")
end

@testset "Transformation between CRS" begin
    # taken from http://osgeo-org.1560.x6.nabble.com/PROJ-PROJ-6-0-0-proj-create-operation-factory-context-behavior-td5403470.html
    src = Proj4.proj_create("IGNF:REUN47GAUSSL")   # area : 55.17,-21.42,55.92,-20.76
    tgt = Proj4.proj_create("IGNF:RGAF09UTM20")    # area : -63.2,14.25,-60.73,18.2
    factory = Proj4.proj_create_operation_factory_context(C_NULL)
    @test factory != C_NULL
    @test Proj4.proj_context_errno() == 0
    results = Proj4.proj_create_operations(src, tgt, factory)
    @test results != C_NULL
    n = Proj4.proj_list_get_count(results)
    for i = 1:n
        operation = Proj4.proj_list_get(results, i - 1)
        hasballpark = Bool(Proj4.proj_coordoperation_has_ballpark_transformation(operation))
        print("""Operation $i
            Name: $(Proj4.proj_get_name(operation))
            Has Ballpark: $hasballpark
        """)
        info = Proj4.proj_pj_info(operation)
        print("""Info $i:
            ID:$(unsafe_string(info.id))
            Desc:$(unsafe_string(info.description))
            Def:$(unsafe_string(info.definition))
            Inv:$(Bool(info.has_inverse))
            Acc:$(info.accuracy)
        """)
        Proj4.proj_destroy(operation)
    end
    Proj4.proj_list_destroy(results)

    Proj4.proj_operation_factory_context_destroy(factory)
    Proj4.proj_destroy(src)
    Proj4.proj_destroy(tgt)
end

@testset "Quickstart" begin
    # https://proj.org/development/quickstart.html

    pj = Proj4.proj_create_crs_to_crs(
        "EPSG:4326",  # source
        "+proj=utm +zone=32 +datum=WGS84",  # target, also EPSG:32632
    )

    # This will ensure that the order of coordinates for the input CRS
    # will be longitude, latitude, whereas EPSG:4326 mandates latitude,
    # longitude
    pj_for_gis = Proj4.proj_normalize_for_visualization(pj)
    Proj4.proj_destroy(pj)
    pj = pj_for_gis

    # a coordinate union representing Copenhagen: 55d N, 12d E
    # Given that we have used proj_normalize_for_visualization(), the order of
    # coordinates is longitude, latitude, and values are expressed in degrees.
    a = Proj4.proj_coord(12, 55)
    @test a isa AbstractVector
    @test a isa SVector{4,Float64}
    @test a isa Proj4.Coord
    @test eltype(a) == Float64
    @test length(a) == 4
    @test sum(a) == Inf
    @test isbits(a)
    @test a[1] === 12.0
    @test a[2] === 55.0
    @test a[3] === 0.0
    @test a[4] === Inf

    # transform to UTM zone 32
    b = Proj4.proj_trans(pj, Proj4.PJ_FWD, a)
    @test b[1] ≈ 691875.632
    @test b[2] ≈ 6098907.825
    @test b[3] === 0.0
    @test b[4] === Inf

    # inverse transform, back to geographical
    b = Proj4.proj_trans(pj, Proj4.PJ_INV, b)
    @test b[1] ≈ 12.0
    @test b[2] ≈ 55.0
    @test b[3] === 0.0
    @test b[4] === Inf

    # Clean up
    pj = Proj4.proj_destroy(pj)
    pj = Proj4.proj_destroy(pj)
    # Julia crashes if proj_destroy is called twice on the same non nullptr,
    # since the contents at that address are already wiped. To prevent this,
    # change the input binding directly after, e.g. pj = proj_destroy(pj)  # -> C_NULL
end

@testset "inverse transformation" begin
    trans = Proj4.Transformation("EPSG:4326", "+proj=utm +zone=32 +datum=WGS84")

    # for custom / proj strings, no description can be looked up in the database
    @test repr(trans) == """
    Transformation
        source: WGS 84
        target: unknown"""

    trans⁻¹ = inv(trans)
    trans¹ = inv(trans⁻¹)

    # inverse twice should get the same transform back
    wkt_type = Proj4.PJ_WKT2_2019
    @test Proj4.proj_as_wkt(trans⁻¹.pj, wkt_type) != Proj4.proj_as_wkt(trans.pj, wkt_type)
    @test Proj4.proj_as_wkt(trans¹.pj, wkt_type) == Proj4.proj_as_wkt(trans.pj, wkt_type)
end

@testset "single transformation" begin
    source_crs = Proj4.proj_create("EPSG:4326")
    target_crs = Proj4.proj_create("EPSG:28992")
    @test source_crs isa Ptr{Nothing}
    @test target_crs isa Ptr{Nothing}
    trans = Proj4.Transformation(source_crs, target_crs)

    # Check that altitude and time inputs are correctly forwarded from the transformation
    @test trans(SA_F64[52.16, 5.39, 5, 2020]) ≈ SA[155191.3538124342, 463537.1362732911, 5.0, 2020.0]

    a = Proj4.proj_coord(52.16, 5.39)
    b = Proj4.proj_trans(trans.pj, Proj4.PJ_FWD, a)
    @test a != b
    @test b ≈ Proj4.proj_coord(155191.3538124342, 463537.1362732911)

    # with always_xy=True, we need to use lon/lat, and still get x/y out
    trans = Proj4.Transformation(source_crs, target_crs, always_xy = true)
    a = Proj4.proj_coord(5.39, 52.16)
    b = Proj4.proj_trans(trans.pj, Proj4.PJ_FWD, a)
    @test a != b
    @test b ≈ Proj4.proj_coord(155191.3538124342, 463537.1362732911)
end

@testset "dense 4D coord vector transformation" begin
    source_crs = Proj4.proj_create("EPSG:4326")
    target_crs = Proj4.proj_create("EPSG:28992")
    trans = Proj4.Transformation(source_crs, target_crs, always_xy = true)
    # This array is mutated in place. Note that this array needs to have 4D elements,
    # with 2D elements it will only do every other one
    A = [Proj4.proj_coord(5.39, 52.16) for _ = 1:5]
    err = Proj4.proj_trans_array(trans.pj, Proj4.PJ_FWD, length(A), A)
    @test err == 0
    B = [Proj4.proj_coord(155191.3538124342, 463537.1362732911) for _ = 1:5]
    @test A ≈ B

    # since A is not in target_crs, we will get an error on a second call
    err = Proj4.proj_trans_array(trans.pj, Proj4.PJ_FWD, length(A), A)
    errno = Proj4.proj_errno(trans.pj)
    @test err == errno == -14
    @test Proj4.proj_errno_string(errno) == "latitude or longitude exceeded limits"
    # reset error
    Proj4.proj_errno_reset(trans.pj)
    @test Proj4.proj_errno(trans.pj) == 0
    # test triggering finalizer manually
    finalize(trans)
    @test trans.pj === C_NULL
end

@testset "generic array transformation" begin
    source_crs = Proj4.proj_create("EPSG:4326")
    target_crs = Proj4.proj_create("EPSG:28992")
    trans = Proj4.Transformation(source_crs, target_crs, always_xy = true)

    # inplace transformation of vector of 2D coordinates
    # using https://proj.org/development/reference/functions.html#c.proj_trans_generic
    A = [SA[5.39, 52.16] for _ = 1:5]
    st = sizeof(first(A))
    ptr = pointer(A)
    n = length(A)
    # 8 is sizeof(Float64), so ptr + 8 points to first latitude element
    n_done = Proj4.proj_trans_generic(
        trans.pj,
        Proj4.PJ_FWD,
        ptr,
        st,
        n,
        ptr + 8,
        st,
        n,
        C_NULL,
        C_NULL,
        C_NULL,
        C_NULL,
        C_NULL,
        C_NULL,
    )
    @test n_done == n
    B = [SA[155191.3538124342, 463537.1362732911] for _ = 1:5]
    @test A ≈ B
end

@testset "compose" begin
    trans1 = Proj4.Transformation("EPSG:4326", "EPSG:28992", always_xy = true)
    trans2 = Proj4.Transformation("EPSG:32632", "EPSG:2027", always_xy = true)
    trans = trans1 ∘ trans2  # same as compose(trans1, trans2)
    source_crs = Proj4.proj_get_source_crs(trans.pj)
    target_crs = Proj4.proj_get_target_crs(trans.pj)

    # trans1 source is the new source
    source_crs1 = Proj4.proj_get_source_crs(trans1.pj)
    source_wkt = Proj4.proj_as_wkt(source_crs, Proj4.PJ_WKT2_2019, C_NULL)
    source_wkt1 = Proj4.proj_as_wkt(source_crs1, Proj4.PJ_WKT2_2019, C_NULL)
    @test source_wkt == source_wkt1

    # trans2 source is the new target
    target_crs1 = Proj4.proj_get_target_crs(trans2.pj)
    target_wkt = Proj4.proj_as_wkt(target_crs, Proj4.PJ_WKT2_2019, C_NULL)
    target_wkt1 = Proj4.proj_as_wkt(target_crs1, Proj4.PJ_WKT2_2019, C_NULL)
    @test target_wkt == target_wkt1

    # which we can also see from show
    @test repr(trans) == """
    Transformation
        source: WGS 84 (with axis order normalized for visualization)
        target: NAD27(76) / UTM zone 15N"""
end

@testset "in and output types" begin
    trans = Proj4.Transformation("EPSG:4326", "EPSG:28992", always_xy = true)
    trans(Proj4.proj_coord(5.39, 52.16))
    b = trans(SA[5.39, 52.16, 0.0, 0.0])

    # StaticVector, Vector or Tuple
    x, y = 5.39, 52.16
    for inpoint in [SA[x, y], [x, y], (x, y)]
        p = if inpoint isa Vector
            Test.@inferred SArray{S, Float64, 1, L} where {S<:Tuple, L} trans(inpoint)
        else
            Test.@inferred trans(inpoint)
        end
        @test p ≈ [155191.3538124342, 463537.1362732911]
        @test p isa SVector{2, Float64}
    end

    # StaticVectors like GeometryBasics.Point are returned, also with Float32
    p = Test.@inferred trans(Point(5.39f0, 52.16f0))
    @test p ≈ [155191.3538124342, 463537.1362732911]
    @test p isa Point{2, Float32}

    # Integer input will go to the ::Any method and become Float64
    @test trans(SA[1,2]) isa SVector{2, Float64}
    @test trans([5,52]) ≈ [128410.08537081012, 445806.50883314764]

    # SVector{3, Float64}
    p = Test.@inferred trans(SA[5.39, 52.16, 2.0])
    @test p ≈ [155191.35381147722, 463537.13624246384, 2.0]
    @test p isa SVector{3, Float64}

    # SVector{4, Float64}
    p = Test.@inferred trans(SA[5.39, 52.16, 0.0, Inf])
    @test p ≈ [155191.3538124342, 463537.1362732911, 0.0, Inf]
    @test p isa SVector{4, Float64}
end

@testset "network" begin
    as_before = Proj4.network_enabled()

    # turn off network, no z transformation
    @test !Proj4.enable_network!(false)
    @test !Proj4.network_enabled()
    trans_z = Proj4.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy = true)
    @test trans_z((151, -33, 5))[3] == 5
    # turn on network, z transformation
    @test Proj4.enable_network!(true)
    @test Proj4.network_enabled()
    trans_z = Proj4.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy = true)
    z = trans_z((151, -33, 5))[3]
    if isinf(z)
        # TODO on CI this hits julia 1.3 all OS, julia 1.6 and nightly Ubuntu only
        @warn "networking not configured correctly"
    else
        @test trans_z((151, -33, 5))[3] ≈ 5.28067830334755
    end

    # 0 args turns it on as well
    Proj4.enable_network!(false)
    Proj4.enable_network!()
    @test Proj4.network_enabled()

    # restore setting as outside the testset
    Proj4.enable_network!(as_before)
end


text:
@testset "absolute.jl" begin
    @testset "Absolute()-Suite1" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(-5, 5, "x", trailer)
        y = SeaPearl.IntVar(-10, 10, "y", trailer)

        constraint = SeaPearl.Absolute(x, y, trailer)

        @test constraint in x.onDomainChange
        @test constraint in y.onDomainChange
        @test y.domain.min.value == -10
        @test y.domain.max.value == 10
        @test x.domain.min.value == -5
        @test x.domain.max.value == 5
        @test constraint.active.value

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test y.domain.max.value == 5
        @test length(x.domain) == 11

        SeaPearl.removeAbove!(x.domain, -2)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
        @test y.domain.min.value == 2

        SeaPearl.removeBelow!(x.domain, -4)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
        @test y.domain.max.value == 4

    end

    @testset "Absolute()-Suite2" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(-5, 0, "x", trailer)
        y = SeaPearl.IntVar(4, 4, "y", trailer)

        constraint = SeaPearl.Absolute(x, y, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.isbound(y)
        @test SeaPearl.isbound(x)
        @test SeaPearl.assignedValue(x) == -4
    end

    @testset "Absolute()-Suite3" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(7, 7, "x", trailer)
        y = SeaPearl.IntVar(-1000, 13, "y", trailer)

        constraint = SeaPearl.Absolute(x, y, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.isbound(x)
        @test SeaPearl.isbound(y)
        @test SeaPearl.assignedValue(x) == 7
    end

    @testset "Absolute()-Suite4" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntVar(7, 7, "a", trailer)
        b = SeaPearl.IntVar(-7, -7, "b", trailer)
        c = SeaPearl.IntVar(6, 13, "c", trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        x = SeaPearl.IntVar(-5, 5, "x", trailer)
        constraint1 = SeaPearl.Absolute(x, a, trailer)
        @test !SeaPearl.propagate!(constraint1, toPropagate, prunedDomains)

        x = SeaPearl.IntVar(-5, 5, "x", trailer)
        constraint2 = SeaPearl.Absolute(b, x, trailer)
        @test !SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)

        x = SeaPearl.IntVar(-5, 5, "x", trailer)
        constraint3 = SeaPearl.Absolute(x, c, trailer)
        @test !SeaPearl.propagate!(constraint3, toPropagate, prunedDomains)
    end

end

ext: jl


text:
using DSP, Distributions

t = linspace(0,5pi,512);
x1 = sin(t);
x2 = 0.4*cos(4t);
x3 = 0.2*cos(8t);
x = x1 + x2 + x3;

responsetype = Lowpass(0.2)
prototype = Elliptic(4, 0.5, 30)
tf = convert(TFFilter, digitalfilter(responsetype, prototype))
numerator_coefs   = coefb(tf)
denominator_coefs = coefa(tf)

responsetype = Bandpass(10, 40; fs=1000)
prototype = Butterworth(4)
xb = filt(digitalfilter(responsetype, prototype), x);

using Winston
plot(xb)

max_stars_repo_path: Module 3/Chapter06/main-dsp.jl

max_issues_repo_path: Module 3/Chapter06/main-dsp.jl

max_forks_repo_path: Module 3/Chapter06/main-dsp.jl


text:
module Pyroheliometer
using Dates
using DataFrames, CSV
using FITSIO
using Statistics

function make_pyrohelio_file_dataframe(pyro_path::AbstractString)
	pyro_files = readdir(pyro_path)
	filter!(f->match(r"neid_ljpyrohelio_chv0_(\d+).tel",f)!=nothing , pyro_files )
	df_pyrohelio_files = DataFrame(:filename=>joinpath.(pyro_path,pyro_files), :start_date=>map(f->Date(match(r"neid_ljpyrohelio_chv0_(\d+).tel$",f)[1],DateFormat("yyyymmdd")),pyro_files))
end

begin
	function pick_pyrohelio_file(pyrohelio_files::DataFrame, fn::AbstractString)
		m = match(r"neidL2_(\d+T\d+).fits",fn)
		@assert m != nothing
		t_start = DateTime(m[1],DateFormat("yyyymmddTHHMMSS"))
		pick_pyrohelio_file(pyrohelio_files,t_start)
	end

	function pick_pyrohelio_file(pyrohelio_files::DataFrame, t::DateTime)
		if !issorted(pyrohelio_files, :start_date)
			sort!(pyrohelio_files, :start_date)
		end
		idx = searchsortedlast(pyrohelio_files.start_date, t)
		if !(1 <= idx <= size(pyrohelio_files,1))
			throw("pick_pyrohelio_file: idx = " * string(idx))
		end
		return pyrohelio_files.filename[idx]
	end
end

# Cache so don't reread CSV files every time
pyrohelio_cache = Dict{String,DataFrame}()

function clear_pyrohelio_cache!()
	global pyrohelio_cache = Dict{String,DataFrame}()
end

begin
	function read_pyrohelio_file(fn::AbstractString; force::Bool = false)
                global pyrohelio_cache
		if !force && haskey(pyrohelio_cache,fn)
			df_pyrohelio_data = pyrohelio_cache[fn]
		else
			df_pyrohelio_data = CSV.read(fn, header=[:time_str,:Voltage,:FluxDensity], types=[String,Float64,Float64], DataFrame)

			datefmt = DateFormat("yyyy-mm-ddTHH:MM:SS.sss")
			df_pyrohelio_data.Time = map(tstr->DateTime(view(tstr,1:23),datefmt),collect(df_pyrohelio_data.time_str))
			#delete!(pyrohelio_data_from_prev_file,[:time_str])
			pyrohelio_cache[fn] = df_pyrohelio_data
		end
		df_pyrohelio_data
	end

	function read_pyrohelio_file(pyrohelio_files::DataFrame, t::DateTime)
		fn = pick_pyrohelio_file(pyrohelio_files, t)
		#read_pyrohelio_file(joinpath(pyro_path,fn))
		read_pyrohelio_file(fn)
	end
	function read_pyrohelio_file(pyrohelio_files::DataFrame, fn::AbstractString)
		m = match(r"neidL\d_(\d+T\d+).fits",fn)
		@assert m != nothing
		t_start = DateTime(m[1],DateFormat("yyyymmddTHHMMSS"))
		fn = pick_pyrohelio_file(pyrohelio_files, t_start)
		#read_pyrohelio_file(joinpath(pyro_path,fn))
		read_pyrohelio_file(fn)
	end
end

begin
	function get_pyrohelio_data(pyrohelio_files::DataFrame, t_start::DateTime, exptime::Real)
		t_stop = t_start + Dates.Second(exptime)

		df = read_pyrohelio_file(pyrohelio_files,t_start)
		idx_start = searchsortedfirst(df.Time,t_start)
		if !(1 <= idx_start <= size(df,1))
			throw("get_pyrohelio_data: idx_start = " * string(idx_start))
		end
		idx_stop = searchsortedlast(df.Time,t_stop)
		if !(1 <= idx_stop <= size(df,1))
			throw("get_pyrohelio_data: idx_stop = " * string(idx_stop))
		end
		#delta_idx = searchsortedlast(view(pyrohelio_data_from_prev_file.time,idx_start:size(pyrohelio_data_from_prev_file,1)),t_stop)
		#idx_stop = idx_start+delta_idx-1
		df[idx_start:idx_stop,:]
	end

	function get_pyrohelio_data(pyrohelio_files::DataFrame, fn::AbstractString, exptime::Real)
		m = match(r"neidL\d_(\d+T\d+).fits",fn)
		@assert m != nothing
		t_start = DateTime(m[1],DateFormat("yyyymmddTHHMMSS"))
		#=
                fits = FITS(fn)
		@assert length(fits) >= 1
		hdr = read_header(fits[1])
		exptime = hdr["EXPTIME"]
		close(fits)
                =#
		exptime_sec = round(Int64,exptime)
		get_pyrohelio_data(pyrohelio_files,t_start, exptime_sec)
	end

end

function get_pyrohelio_mean_Δt(df::DataFrame; time_start::DateTime = first(df.Time) )
	@assert hasproperty(df,:Time)
	@assert hasproperty(df,:FluxDensity)
	dt = map(dt->dt.value/1000,df.Time .-time_start)
	sum(dt.*df.FluxDensity)/sum(df.FluxDensity)
end

begin
	function get_pyrohelio_summary(pyrohelio_files::DataFrame, fn::AbstractString, exptime::Real)
		m = match(r"neidL\d_(\d+T\d+).fits",fn)
		@assert m != nothing
		time_start = DateTime(m[1],DateFormat("yyyymmddTHHMMSS"))
		#=
        	fits = FITS(fn)
			@assert length(fits) >= 1
			hdr = read_header(fits[1])
			exptime = hdr["EXPTIME"]
			close(fits)
        	=#
		exptime_sec = round(Int64,exptime)
		#println(fn,": exptime = ", exptime_sec)
		try
			df = get_pyrohelio_data(pyrohelio_files, time_start, exptime_sec)
			#println(fn, ": pyrheliometer_flux = ", size(df.FluxDensity,1))
			#time_start = first(df.time)  # Comment out so use FITS file start as reference time
			mean_Δt = get_pyrohelio_mean_Δt(df, time_start=time_start)
			mean_flux = mean(df.FluxDensity)
			rms_flux = sqrt(var(df.FluxDensity,mean=mean_flux,corrected=false))
			(min_flux, max_flux) = extrema(df.FluxDensity)
			(;filename=basename(fn), time_start, exptime, mean_Δt, mean_pyroflux=mean_flux, rms_pyroflux=rms_flux, min_pyroflux=min_flux, max_pyroflux=max_flux)
		catch exepction
			(;filename=basename(fn), time_start, exptime, mean_Δt=missing, mean_pyroflux=missing, rms_pyroflux=missing,  min_pyroflux=missing, max_pyroflux=missing)
		end
	end

	function get_pyrohelio_summary(fn::String) #, exptime::Real)
		m = match(r"neidL\d_(\d+T\d+)\.fits",fn)
		@assert m != nothing
		t_start = DateTime(m[1],DateFormat("yyyymmddTHHMMSS"))
		try
			f = FITS(fn)
			@assert length(f) >= 1
			hdr = read_header(f[1])
			exptime = hdr["EXPTIME"]
			hdr_pyrhelio = read_header(f["SOLAR"])
			stel = hdr_pyrhelio["STEL"]
			staz = hdr_pyrhelio["STAZ"]
			#println(fn,": exptime = ", exptime)

			t = read(f["SOLAR"],"Time")
        	v = read(f["SOLAR"],"Voltage")
        	fd = read(f["SOLAR"],"FluxDensity")
        	df_tmp = DataFrame(:Time=> DateTime.(t), :Voltage=>v, :FluxDensity=>fd)
			#println(fn, ": pyrheliometer_flux = ", size(df_tmp.FluxDensity,1))
			#df_tmp.Time = DateTime.(df_tmp.Time)
        	mean_Δt = get_pyrohelio_mean_Δt(df_tmp)
			mean_pyroflux = mean(df_tmp.FluxDensity)
        	rms_pyroflux = sqrt(var(df_tmp.FluxDensity, corrected=false))
        	(min_pyroflux,max_pyroflux) = extrema(df_tmp.FluxDensity)
			close(f)
			return (;filename=basename(fn), time_start=t_start, exptime, mean_Δt, mean_pyroflux, rms_pyroflux, min_pyroflux, max_pyroflux, stel, staz)
		catch exepction
			return (;filename=basename(fn), time_start=t_start, exptime=missing, mean_Δt=missing, mean_pyroflux=missing, rms_pyroflux=missing,  min_pyroflux=missing, max_pyroflux=missing, stel=missing, staz=missing)
		end
	end

end

export make_pyrohelio_file_dataframe, pick_pyrohelio_file, read_pyrohelio_file
export get_pyrohelio_data, get_pyrohelio_summary
end # module


text:
using FreeRay
using Test, Plots

#### Geometry
# Point source 1000 Hz at 80 m
source = Source(
    frequency = 100,
    height = -5
)

# Multiple recievers on grid
receiver = Receiver(
    depth_point = 1001,
    range_point = 10001,
    depth = Vec2(-1000f0,0f0),
    range = Vec2(0f0,10f0)
)

# Terrain geometry
terrain = Terrain(
    interp_type = "C",
    profile = (Vec2([0f0,5f0,10f0],[0f0,0f0,0f0]))
)

### Boundary condition
Zc = 12.81 + 11.62im
Theta, Rmag, Rphase = R_coeff(Zc;len=100)

brc = Vec3(Theta,Rmag,Rphase)
trc = Vec3([0f0,45f0,90f0],[0f0,0f0,0f0],[0f0,0f0,0f0])

reflection = Reflection_Coeff(
    top_coeff = trc,
    bottom_coeff = brc
)

### Sound speed profile
#sspl = Vec2([-1000f0,-500f0,0f0],[243f0,293,343]) # case 3
sspl = Vec2([-1000f0,-500f0,0f0],[343f0,343f0,343f0]) # case 1
#sspl = Vec2([-1000f0,-500f0,0f0],[443f0,393f0,343f0]) # case 2
#sspl = Vec2([-1000f0,-300f0,-100f0,0f0],[333f0,333f0,353f0,343f0]) #case 4
ssp = SSP(sound_speed_profile= sspl)

### Analysis
opt = Analysis(
    filename = "Case0_Bellhop_f10",
    analyse = "RG",
    option1 = "CFW",
    option2 = "F*",
    num_ray = 161,
    alpha = Vec2(-80.0f0,80.0f0),
    box = Vec2(10f0,1000.0f0),
    step= 0
)

Environment(opt,source, receiver,ssp,terrain,reflection)

# Run Bellhop
fn= opt.filename
filename = "$fn"
run_bellhop = `bellhop $filename`
@time run(run_bellhop)

p1 = PlotRay("$filename.ray",
        xlabs = "Range, m",
        ylabs = "Height, m")

#savefig(p1,"ray1.png")

#p2 = PlotShd("$filename.shd";
#        xlabs = "Range, m",
#        ylabs = "Transmission loss, dB",
#        cblabs = "dB",
#        climb = (40,80))

#savefig(p2,"transc4.png")

#PlotTlr("$filename.shd", -7;
#        xlabs = "Range, m",
#        ylabs = "Transmission loss, dB",
#        xlim = (0,10000),
#        ylim = (20,120))

rmfile(filename)

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

#read_shd("Case0_Bellhop_f10.shd")


text:
module test_Plot2

import ..ModiaMath
using Test
t = range(0.0, stop=10.0, length=100)

series = Dict{AbstractString,Any}()
series["time"] = t
series["phi"]  = sin.(t)
series["w"]    = cos.(t)
series["r"]    = hcat(0.4 * cos.(t), 0.5 * sin.(t), 0.3*cos.(t))

time = ModiaMath.RealScalar(:time, unit="s")
phi  = ModiaMath.RealScalar(:phi, unit="rad")
w    = ModiaMath.RealScalar(:w, unit="rad/s")
r    = ModiaMath.RealSVector3(:r, unit="m")
var  = Dict{Symbol,Any}(:time => time, :phi => phi, :w => w, :r => r)

result = ModiaMath.ResultWithVariables(series, var, "ModiaMath/test/Result/test_Plot2.jl")

ModiaMath.plot(result, :phi    , prefix="sim 1: ", heading="Sine(time)")
ModiaMath.plot(result, "r"     , prefix="sim 1: ", figure=2)
ModiaMath.plot(result, "r[2]"  , figure=3)

println("\n... Next plot should give a warning:")
ModiaMath.plot(result, "r[2:3]", figure=4)

println("\n... Next plot should give a warning:")
ModiaMath.plot(result, "phi", xAxis="r", figure=4)

println("\n... figure=4 is closed")
ModiaMath.closeFigure(4)

# Print result variables
println("\n... result variables = ", result)

# Add next simulation run to plot
result.series["phi"] = 1.2*sin.(t)
result.series["r"]   = hcat(0.2 * cos.(t), 0.3 * sin.(t), 0.5*cos.(t))
ModiaMath.plot(result, :phi, heading="Sine(time)", figure=1, prefix="sim 2: ", reuse=true)
ModiaMath.plot(result, :r  ,                       figure=2, prefix="sim 2: ", reuse=true)

end


text:
using BayesianLinearRegression
using Test
using LinearAlgebra

const BLR = BayesianLinearRegression

function runtests()
    X = randn(100,30)
    β = randn(30)
    y = X * β .+ randn(100)

    m = BayesianLinReg(X,y)
    fit!(m)

    @testset "log-determinant of Hessian" begin
        @test BLR.logdetH(m) ≈ logdet(hessian(m))
    end

    @testset "Inverse Hessian" begin
        @test m.Hinv * hessian(m) ≈ I
    end

    @testset "Sum of Squared Residuals" begin
        @test BLR.ssr(m) ≈ BLR.normSquared(y - predict(m,X; uncertainty=false))
    end

    # This comes from the equating γ = γ between the α and β updates
    @testset "N invariant" begin
        α = priorPrecision(m)
        β = noisePrecision(m)
        w² = BLR.normSquared(m.weights)
        ssr = BLR.ssr(m)
        @test α * w² + β * ssr ≈ m.N
    end
end

runtests()


text:
# Analysis

using Clustering
using Statistics
using StatsBase

"""
x,y,photons=getxy(states::Vector{Any})

returns a vector of all x,y,photons over all states
"""
function getxy(states::Vector{Any})
    # count number of emitters
    nemitters = 0;
    for nn = 1:length(states)
        nemitters += states[nn].n
    end

    x = Vector{Float32}(undef, nemitters)
    y = Vector{Float32}(undef, nemitters)
    photons = Vector{Float32}(undef, nemitters)

    cnt = 0;
    for ss = 1:length(states)
        for nn = 1:states[ss].n
            cnt += 1
            x[cnt] = states[ss].x[nn]
            y[cnt] = states[ss].y[nn]
            photons[cnt] = states[ss].photons[nn]
        end
    end
    return x,y,photons
end

function getn(states::Vector{Any})
    maxemitters = 0
    nstates=length(states)

    traj_n=Vector{Int32}((undef),nstates)
    for nn = 1:nstates
        traj_n[nn]=states[nn].n
        maxemitters=max(maxemitters,states[nn].n)
    end

    posterior_n=zeros(Float32, maxemitters+1)
    for nn = 1:nstates
        posterior_n[states[nn].n+1]+=states[nn].n
    end
    posterior_n=posterior_n/sum(posterior_n)

    cntmax = maximum(posterior_n)
    coords = findall(x -> x == cntmax, posterior_n)
    map_n=Int(coords[1]-1)

    return map_n,posterior_n,traj_n
end

"""
    getnstates(states::Vector{Any})

return states with n number of emitters
"""
function getnstates(states::Vector{Any},n::Int)
    nstates=length(states)
    idx=zeros(Bool,nstates)

    for ii=1:length(states)
        if states[ii].n==n
            idx[ii]=true
        end
    end

    return states[idx]
end

"""
states_mapn=getmapn(states::Vector{Any})

returns a vector of states that where states.n == map_n
"""
function getmapnstates(states::Vector{Any})
    map_n,posterior_n,traj_n=getn(states)

    nstates=length(states)
    idx=zeros(Bool,nstates)

    for ii=1:length(states)
        if states[ii].n==map_n
            idx[ii]=true
        end
    end

    return states[idx],map_n
end

"""
    clusterstates(states::Vector{Any},n::Int32)

uses kmeans clustering group x,y locations into n clusters
and returns a State_Results structure.

"""
function clusterstates(states::Vector{Any},n::Int)
    x,y,photons=getxy(states)
    X=cat(dims=1,x',y')
    R = kmeans(X, n; maxiter=200, display=:iter)
    @assert nclusters(R) == n # verify the number of clusters
    a = assignments(R) # get the assignments of points to clusters
    c = counts(R) # get the cluster sizes
    M = R.centers # get the cluster centers

    state_mapn=StateFlatBg_Results(n)

    for ii=1:n
        x_n=X[1,a.==ii]
        y_n=X[2,a.==ii]
        p_n=photons[a.==ii]
        state_mapn.x[ii]=mean(x_n)
        state_mapn.y[ii]=mean(y_n)
        state_mapn.photons[ii]=mean(p_n)
        state_mapn.σ_x[ii]=std(x_n)
        state_mapn.σ_y[ii]=std(y_n)
        state_mapn.σ_photons[ii]=std(p_n)
    end
    return state_mapn
end

function getmapn(states::Vector{Any})
    mapnstates,n=getmapnstates(states)
    return clusterstates(mapnstates,n)
end

"""
    getposterior(states::Vector{Any}, sz::Int32, zoom::Int32)

Get a 2D posterior probability image of source locations

# Arguments
- `states::Vector{Any}`: Array of states found by RJMCMC`
- `sz::Int32' : Linear size of analyzed ROI (pixels)
- `zoom::Int32' : Pixel subsampling

"""
function getposterior(states::Vector{Any}, sz::Int, zoom::Int)
    x,y=getxy(states)
    xbins = range(0.5, sz+0.5;step=1 / zoom)
    ybins = range(0.5, sz+0.5;step=1 / zoom)
    h = fit(Histogram,(x,y),(xbins,ybins))
    return h.weights
end


text:
__precompile__()

module Simulation_Env_module

push!(LOAD_PATH, "./")
push!(LOAD_PATH, "../")
push!(LOAD_PATH, "../../VICAS/src/")

using DeepRL
using POMDPs
using PyPlot
using CAS_QMDP
using Distributions

include("common_const.jl")

export Simple_CASIM_Env, Dummy_Problem
export step!, sample_action
export convert_ac_state, get_qvals, lowfi_values, normalize_qvals, normalize_obs!, actions
export IND_rho, IND_theta, IND_phi, IND_v_own, IND_v_int, IND_deviation, IND_dist_to_dest, IND_prev_dist_to_dest
export STATE_DIM, AUG_STATE_DIM, TERM_VAR, LOW_FI_POLICY, LOW_FI_POLICY_FILENAME, DEST_RANGE
export V_MAX, NMAC_RANGE, SENSING_RANGE, SPAWN_RATE_SET, CONST_NUM_AC_SET, RELAX_FACTOR
export POSSIBLE_MAX_REW_VICAS, POSSIBLE_MIN_REW_VICAS
export terminal_state

"""
	keeps the angles in [0, 2*pi]
"""
function norm_angle(angle::Float64)
	return ((angle % (2 * pi)) + 2 * pi) % (2 * pi)
end

"""
	state of an Aircraft object, reflecting the observation of the environemnt
	s = [ρ, θ, ϕ, v_own, v_int]
"""
mutable struct Aircraft_State
	rho::Float64
	theta::Float64
	phi::Float64
	v_own::Float64
	v_int::Float64
end

"""
	convert Aircraft_State to Float64 array
	s = [ρ, θ, ϕ, v_own, v_int]
"""
function convert_ac_state(s::Aircraft_State)
	return Float64[s.rho, s.theta, s.phi, s.v_own, s.v_int]
end

"""
	Define the dynamics of the Aircraft, self property
"""
mutable struct Dynamics
	x::Float64
	y::Float64
	dt::Float64
	v::Float64
	heading::Float64
	turn_rate::Float64
end

"""
	agent object
"""
mutable struct Aircraft
	id::Base.Random.UUID
	is_ego::Bool
	dist_to_ego::Float64
	section::Int64
	x_init::Float64
	y_init::Float64
	x_dest::Float64
	y_dest::Float64
	advisory::Symbol
	dyn::Dynamics
	state::Aircraft_State
end # module

"""
	initiate the ego agent, i.e., the learning agent
"""
function init_ego_agent()
	id = Base.Random.uuid1()
	is_ego = true
	section = 0
	dist_to_ego = 0.

	x_init = 0.
	y_init = 0.
	th = 2 * pi * rand()
	x_dest = cos(th) * DEST_RANGE
	y_dest = sin(th) * DEST_RANGE

	# dynamics
	x = x_init
	y = y_init
	v = (V_MAX - V_MIN) * rand() + V_MIN
	heading = norm_angle(atan2(y_dest, x_dest))
	turn_rate = 0.

	advisory = :COC
	dyn = Dynamics(x, y, DT, v, heading, turn_rate)
	state = Aircraft_State(0., 0., 0., 0., 0.)
	ego_agent = Aircraft(id, is_ego, dist_to_ego, section, x_init, y_init,
							x_dest, y_dest, advisory, dyn, state)
	return ego_agent
end

"""
	Generate init and dest coords for env agents
	This has to be truly representitive of the CASIM simulation
"""
function get_init_and_dest_coords(ego_agent::Aircraft)
	# Use encounter model:
	intruder_dist = Categorical(vec(BEARING_DIST[:, 2]))
	intruder_angle = BEARING_DIST[:, 1][rand(intruder_dist)]
	if intruder_angle > π
		intruder_angle -= 2 * π
	end # intruder_angle ∈ [-pi, pi], centered at 0
	init_angle = norm_angle(ego_agent.dyn.heading + intruder_angle)
	dest_angle = norm_angle(init_angle + deg2rad(rand() * 270 + 45)) # 45 - 315 uniform?
	x_init = cos(init_angle) * SENSING_RANGE * RELAX_FACTOR
	y_init = sin(init_angle) * SENSING_RANGE * RELAX_FACTOR
	x_dest = cos(dest_angle) * SENSING_RANGE * RELAX_FACTOR
	y_dest = sin(dest_angle) * SENSING_RANGE * RELAX_FACTOR

	# Use Gaussian:
	# RELAX_FACTOR = 1.5
	# init_angle = norm_angle(ego_agent.dyn.heading + rand(Normal(0, deg2rad(70))))
	# dest_angle = norm_angle(init_angle + rand(Normal(pi, deg2rad(90)))) # TODO
	# x_init = cos(init_angle) * SENSING_RANGE * RELAX_FACTOR
	# y_init = sin(init_angle) * SENSING_RANGE * RELAX_FACTOR
	# x_dest = cos(dest_angle) * SENSING_RANGE * RELAX_FACTOR
	# y_dest = sin(dest_angle) * SENSING_RANGE * RELAX_FACTOR

	# init_angle = norm_angle(ego_agent.dyn.heading + rand(Normal(0, deg2rad(90))))
	# dest_angle = norm_angle(init_angle + rand(Normal(pi, deg2rad(100)))) # TODO

	# # @ outdated @
	# x_init = cos(init_angle) * SENSING_RANGE * RELAX_FACTOR * (1 - rand() * 0.4)
	# y_init = sin(init_angle) * SENSING_RANGE * RELAX_FACTOR * (1 - rand() * 0.4)
	# x_dest = cos(dest_angle) * SENSING_RANGE * RELAX_FACTOR * (1 - rand() * 0.4)
	# y_dest = sin(dest_angle) * SENSING_RANGE * RELAX_FACTOR * (1 - rand() * 0.4)

	return x_init, y_init, x_dest, y_dest
end

"""
	initiate single env agent
"""
function init_one_env_agent(ego_agent::Aircraft, num_sections::Int64)
	id = Base.Random.uuid1()
	is_ego = false
	section = 0

	x_init, y_init, x_dest, y_dest = get_init_and_dest_coords(ego_agent)

	# dynamics
	x = x_init
	y = y_init
	v = (V_MAX - V_MIN) * rand() + V_MIN
	# adding noise to the initial heading of env_agents
	heading = norm_angle(atan2(y_dest - y_init, x_dest - x_init) + rand(Normal(0, deg2rad(20))))
	turn_rate = 0.
	dist_to_ego = sqrt(x^2 + y^2)

	advisory = :COC
	dyn = Dynamics(x, y, DT, v, heading, turn_rate)
	state = Aircraft_State(0., 0., 0., 0., 0.)
	return Aircraft(id, is_ego, dist_to_ego, section, x_init, y_init,
						x_dest, y_dest, advisory, dyn, state)
end

"""
	initiate an array of env agents
"""
function init_env_agents(ego_agent::Aircraft, init_num_env_agents::Int64, num_sections::Int64)
	env_agents = Vector{Aircraft}(init_num_env_agents)
	for i in 1:init_num_env_agents
		env_agents[i] = init_one_env_agent(ego_agent, num_sections)
	end
	return env_agents
end

"""
	define an interface helper
"""
struct Dummy_Problem
	discount_factor::Float64
end

"""
	define the simulation environemnt object
"""
mutable struct Simple_CASIM_Env <: AbstractEnvironment
	scheme::AbstractString
	t::Int64
	sim_horizon::Int64
	num_sections::Int64
	init_num_env_agents::Int64
	const_num_env_agents::Int64
	ego_agent::Aircraft
	ego_agent_dist_to_dest::Float64
	env_agents::Vector{Aircraft}
	section_ac_arr::Dict{Int64, Vector{Aircraft}}
	full_state::Vector{Float64}
	action_space::Vector{Symbol}
	advisory_to_ind_dict::Dict{Symbol, Int64}
	advisory_to_action_dict::Dict{Symbol, Float64}
	discount_factor::Float64
	low_fi_policy::Matrix{Float64}
	reward_weights::Dict{AbstractString, Float64}
	stochastic_env_policy::Bool
	stochastic_ego_policy::Bool
	# POMDPs
	rng::AbstractRNG
	problem::Dummy_Problem
	# visualization option
	render::Bool
	# deep correction
	correction::Bool
	correction_weight::Float64
	nn_policy::Bool
	random_populate_replay_buffer::Bool

	function Simple_CASIM_Env(;
		scheme::AbstractString="_TEST_",
		sim_horizon::Int64=5000,
		num_sections::Int64=NUM_SECTIONS,
		init_num_env_agents::Int64=6,
		const_num_env_agents::Int64=rand(CONST_NUM_AC_SET),
		discount_factor::Float64=0.95,
		low_fi_policy_filename::AbstractString=LOW_FI_POLICY_FILENAME,
		rng::MersenneTwister=MersenneTwister(1),
		render::Bool=false,
		correction::Bool=false,
		stochastic_env_policy::Bool=false,
		stochastic_ego_policy::Bool=false,
		nn_policy::Bool=false,
		correction_weight::Float64=0.5,
		random_populate_replay_buffer::Bool=false)

		this = new()
		this.scheme = scheme
		this.t = 1
		this.sim_horizon = sim_horizon
		this.num_sections = num_sections
		this.init_num_env_agents = init_num_env_agents
		this.const_num_env_agents = const_num_env_agents
		this.ego_agent = init_ego_agent()
		this.ego_agent_dist_to_dest = DEST_RANGE
		this.env_agents = init_env_agents(this.ego_agent, init_num_env_agents, num_sections)
		this.section_ac_arr = Dict{Int64, Vector{Aircraft}}()
		this.full_state = zeros(num_sections * STATE_DIM + AUG_STATE_DIM)
		this.action_space = Symbol[:SR, :WR, :KEEP, :WL, :SL, :COC]
		this.advisory_to_ind_dict = Dict(:SR => 1, :WR => 2, :KEEP => 3, :WL => 4, :SL => 5, :COC => 6)
		this.advisory_to_action_dict = Dict(:SR => deg2rad(-10.0), :WR => deg2rad(-5.0), :KEEP => 0.0,
											:WL => deg2rad(5.0), :SL => deg2rad(10.0), :COC => -1.0)
		this.discount_factor = discount_factor
		this.low_fi_policy = load_alphas(low_fi_policy_filename)
		# POMDPs
		this.rng = rng
		this.problem = Dummy_Problem(discount_factor)
		this.render = render
		if render == true
			show()
		end
		this.correction = correction
		this.correction_weight = correction_weight
		this.reward_weights = Dict("penalty_action" => PEN_ACTION,
								   "penalty_closeness" => PEN_CLOSENESS,
								   "penalty_nmac" => PEN_NMAC,
								   "penalty_conflict" => PEN_CONFLICT,
								   "penalty_deviation" => PEN_DEVIATION,
								   "penalty_digression" => PEN_DIGRESSION,
								   "reward_destination" => REW_DESTINATION)
		this.stochastic_env_policy = stochastic_env_policy
		this.stochastic_ego_policy = stochastic_ego_policy
		this.nn_policy = nn_policy
		this.random_populate_replay_buffer = random_populate_replay_buffer
		return this
	end
end # struct

"""
	visualize the environemnt
"""
function render_env(env::Simple_CASIM_Env)
	All_AC = vcat(env.env_agents, env.ego_agent)
	for ac in All_AC
		if ac.advisory == :COC
			scatter(ac.dyn.x, ac.dyn.y, marker="o", color="blue", s=12)
		else
			scatter(ac.dyn.x, ac.dyn.y, marker="o", color="red", s=12)
		end

		arrow_len = AIRSPACE_DIM / 6
		arrow(
			ac.dyn.x, ac.dyn.y,
			cos(ac.dyn.heading) * arrow_len * 0.8,
			sin(ac.dyn.heading) * arrow_len * 0.8,
			head_width=AIRSPACE_DIM / 20,
			width=1,
			head_length=AIRSPACE_DIM / 20,
			overhang=0.5,
			head_starts_at_zero="true",
			facecolor="black",
			length_includes_head="true")

		scatter(ac.x_dest, ac.y_dest, marker=",", color="magenta", s=12)

		plot([ac.dyn.x; ac.x_dest], [ac.dyn.y; ac.y_dest],
			linestyle="--", color="black", linewidth=0.5)

		# text(ac.dyn.x, ac.dyn.y, string(ac.section))
	end

	# alert circle
	th = linspace(-pi, pi, 30)
	if env.ego_agent.advisory == :COC
		plot(env.ego_agent.dyn.x + SENSING_RANGE * cos.(th),
			env.ego_agent.dyn.y + SENSING_RANGE * sin.(th),
			linestyle="--", color="green", linewidth=0.4)
	else
		plot(env.ego_agent.dyn.x + SENSING_RANGE * cos.(th),
			env.ego_agent.dyn.y + SENSING_RANGE * sin.(th),
			linestyle="--", color="red", linewidth=0.4)
	end

	# NMAC circle
	plot(env.ego_agent.dyn.x + NMAC_RANGE * cos.(th),
			env.ego_agent.dyn.y + NMAC_RANGE * sin.(th),
			linestyle="--", color="red", linewidth=0.4)

	# plot sections and highlight the closest intruder in each section
	# section_angle = 2 * pi / env.num_sections
	# for i in 0 : env.num_sections / 2 - 1
	# 	th_1 = i * section_angle + env.ego_agent.dyn.heading
	# 	th_2 = th_1 + pi
	# 	plot(SENSING_RANGE * [cos(th_1), cos(th_2)], SENSING_RANGE * [sin(th_1), sin(th_2)],
	# 		linestyle="--", color="red", linewidth=0.4)
	# end

	# println(env.full_state)
	for i in 1 : env.num_sections
		rho = env.full_state[(i - 1) * STATE_DIM + IND_rho]
		if rho != TERM_VAR
			theta = env.full_state[(i - 1) * STATE_DIM + IND_theta]
			scatter(rho * cos(theta), rho * sin(theta),
				marker="o", color="green", s=40, linewidths=0.6)
		end
	end

	xlabel("X (m)")
	ylabel("Y (m)")
	title(string(env.ego_agent.advisory) * string(env.const_num_env_agents / 100) * ", " * string(env.t) * ": " * string(length(env.env_agents)))
	axis("equal")
	ax = gca()
	ax[:set_xlim]([-AIRSPACE_DIM * 1.5, AIRSPACE_DIM * 1.5])
	ax[:set_ylim]([-AIRSPACE_DIM * 1.5, AIRSPACE_DIM * 1.5])
	pause(0.1)
	clf()
end

"""
	overload :< operator for the Aircraft object
"""
function Base.:<(a::Aircraft, b::Aircraft)
	if a.dyn.x^2 + a.dyn.y^2 < b.dyn.x^2 + b.dyn.y^2
		return true
	else
		return false
	end
end

"""
	overload :> operator for the Aircraft object
"""
function Base.:>(a::Aircraft, b::Aircraft)
	if a.dyn.x^2 + a.dyn.y^2 > b.dyn.x^2 + b.dyn.y^2
		return true
	else
		return false
	end
end

"""
	sort the Aircraft in an array, wrt distance to the ego agent
	using quick sort, returns ascending order in terms of the distance to the ego agent
"""
function sort_agents!(agents_arr::Vector{Aircraft}, lo::Int64, hi::Int64)
    i, j = lo, hi
    while i < hi
        pivot = agents_arr[(lo + hi) >>> 1]
        while i <= j
            while agents_arr[i] < pivot; i += 1; end
            while agents_arr[j] > pivot; j -= 1; end
            if i <= j
                agents_arr[i], agents_arr[j] = agents_arr[j], agents_arr[i]
                i, j = i + 1, j - 1
            end
        end
        if lo < j; sort_agents!(agents_arr, lo, j); end
        lo, j = i, hi
    end
end

"""
	sort the Aircraft in an array, wrt the maximum q_vals wrt the eo agent
	using quick sort, returns ascending order in terms of the maximum q_vals
"""
function sort_agents_q_vals!(env::Simple_CASIM_Env, agents_arr::Vector{Aircraft}, lo::Int64, hi::Int64)
	i, j = lo, hi
    while i < hi
        pivot = agents_arr[(lo + hi) >>> 1]
        pivot_max_q = maximum(get_qvals(env.low_fi_policy, get_state(env.ego_agent, pivot)))
        while i <= j
        	agent_i_max_q = maximum(get_qvals(env.low_fi_policy, get_state(env.ego_agent, agents_arr[i])))
        	agent_j_max_q = maximum(get_qvals(env.low_fi_policy, get_state(env.ego_agent, agents_arr[j])))
            while agent_i_max_q < pivot_max_q
            	i += 1
            	agent_i_max_q = maximum(get_qvals(env.low_fi_policy, get_state(env.ego_agent, agents_arr[i])))
            end
            while agent_j_max_q > pivot_max_q
            	j -= 1
            	agent_j_max_q = maximum(get_qvals(env.low_fi_policy, get_state(env.ego_agent, agents_arr[j])))
            end
            if i <= j
                agents_arr[i], agents_arr[j] = agents_arr[j], agents_arr[i]
                i, j = i + 1, j - 1
            end
        end
        if lo < j; sort_agents_q_vals!(env, agents_arr, lo, j); end
        lo, j = i, hi
    end
end

"""
	update the state of an Aircraft object with Nearest_Intruder_Tracker
	OR with the most dangerous intruder traker
	OR with the Multi_Intruder_Tracker
"""
function get_state!(ac::Aircraft, env::Simple_CASIM_Env)
	# search for the closest intruder from ownship
	if length(env.env_agents) > 0

		# # The closest intruder
		# closest_intruder = nothing
		# closest_dist = Inf
		# All_AC = vcat(env.env_agents, env.ego_agent)
		# for other_ac in All_AC
		# 	if other_ac.id != ac.id
		# 		dist = sqrt((ac.dyn.x - other_ac.dyn.x)^2 + (ac.dyn.y - other_ac.dyn.y)^2)
		# 		if dist < closest_dist
		# 			closest_dist = dist
		# 			closest_intruder = other_ac
		# 		end
		# 	end
		# end

		# # The most dangerous intruder in terms of low max q_vals
		closest_intruder = nothing
		min_max_q = Inf
		All_AC = vcat(env.env_agents, env.ego_agent)
		for other_ac in All_AC
			if other_ac.id != ac.id
				max_q = maximum(get_qvals(env.low_fi_policy, get_state(ac, other_ac)))
				if max_q < min_max_q
					min_max_q = max_q
					closest_intruder = other_ac
				end
			end
		end
		closest_dist = sqrt((ac.dyn.x - closest_intruder.dyn.x)^2 + (ac.dyn.y - closest_intruder.dyn.y)^2)

		# s = [ρ, θ, ϕ, v_own, v_int]
		intruder_pos_angle = norm_angle(atan2(closest_intruder.dyn.y - ac.dyn.y,
											  closest_intruder.dyn.x - ac.dyn.x))
		theta = norm_angle(intruder_pos_angle - ac.dyn.heading)
		phi = norm_angle(closest_intruder.dyn.heading - ac.dyn.heading)

		ac.state = Aircraft_State(closest_dist, theta, phi, ac.dyn.v, closest_intruder.dyn.v)
	else # no env_agents
		# following VICAS, set the state to terminal state
		ac.state = Aircraft_State(TERM_VAR, TERM_VAR, TERM_VAR, TERM_VAR, TERM_VAR)
	end
end

"""
	helper: get state for ego_agent
"""
function get_state(ego::Aircraft, int::Aircraft)
	rho = sqrt(int.dyn.x^2 + int.dyn.y^2)
	theta = norm_angle(atan2(int.dyn.y, int.dyn.x))
	phi = norm_angle(int.dyn.heading - ego.dyn.heading)
	return [rho, theta, phi, ego.dyn.v, int.dyn.v]
end

"""
	get section id of an ac
"""
function get_section_id!(ac::Aircraft, env::Simple_CASIM_Env)
	# ac_pos_angle = norm_angle(atan2(ac.dyn.y, ac.dyn.x))
	# Relative to the heading of the ego_agent
	ac_pos_angle = norm_angle(norm_angle(atan2(ac.dyn.y, ac.dyn.x)) - norm_angle(env.ego_agent.dyn.heading))
	ac.section = floor(ac_pos_angle / (2 * pi / env.num_sections)) + 1
end

"""
	update the states of each agent in the env as well as the full_state
"""
function update_states!(env::Simple_CASIM_Env)
	# Two other different ways of forming full_state:
	get_state!(env.ego_agent, env) # update intruder observation: closest or most dangerous

	#
	# intruders = Vector{Aircraft}(0)
	# for ac in env.env_agents
	# 	get_state!(ac, env) # update state for env_agents using closest
	# 	if check_within_sensing(ac, env) # collect intruders for ego_agent
	# 		push!(intruders, ac)
	# 	end
	# end

	# hi = length(intruders)
	# lo = (hi == 0) ? 0 : 1

	# if hi > 0
	# 	# 1. Using multiple closest intruders
	# 	# Ascending order in terms of distance from ego_agent
	# 	sort_agents!(intruders, lo, hi)
	# 	# 2. Using multiple dangerous intruders
	# 	# Ascending order in terms of maximum q_vals wrt ego_agent
	# 	# The minimum maximum q_vals => the most dangerous
	# 	# sort_agents_q_vals!(env, intruders, lo, hi)

	# 	for i in 1 : env.num_sections
	# 		if i <= length(intruders)
	# 			env.full_state[STATE_DIM * (i - 1) + 1 : STATE_DIM * i] .= get_state(env.ego_agent, intruders[i])
	# 		else
	# 			env.full_state[STATE_DIM * (i - 1) + 1 : STATE_DIM * i] .= terminal_state
	# 		end
	# 	end
	# else # no intruder
	# 	env.full_state[1 : STATE_DIM * env.num_sections] .= TERM_VAR
	# end
	#

	# 3. Using closest intruders in each sections
	# clear section_ac_arr

	env.section_ac_arr = Dict{Int64, Vector{Aircraft}}()
	for isection in 1:env.num_sections
		env.section_ac_arr[isection] = Aircraft[]
	end

	# update state for the ego agent
	get_state!(env.ego_agent, env)

	for ac in env.env_agents
		get_state!(ac, env) # update state for env_agents using closest
		if check_within_sensing(ac, env)
			get_section_id!(ac, env)
			push!(env.section_ac_arr[ac.section], ac)
		else
			ac.section = 0
		end
	end

	# sort ac in each section
	for isection in 1:env.num_sections
		hi = length(env.section_ac_arr[isection])
		lo = (hi == 0) ? 0 : 1
		if hi > 0
			sort_agents!(env.section_ac_arr[isection], lo, hi)
			# OR
			# sort_agents_q_vals!(env, env.section_ac_arr[isection], lo, hi)

			# Resulting section_ac_arr is in ascending order,
			# taking ac that is the closest OR the most dangerous (indexed by 1)
			# wrt the ego ac for each section to update the full_state
			env.full_state[STATE_DIM * (isection - 1) + 1 : STATE_DIM * isection] .=
				get_state(env.ego_agent, env.section_ac_arr[isection][1])
		else # no intruder in the section
			# set the corresponding state var to terminal values
			env.full_state[STATE_DIM * (isection - 1) + 1 : STATE_DIM * isection] .=
				[TERM_VAR, TERM_VAR, TERM_VAR, TERM_VAR, TERM_VAR]
		end
	end
	#

	# Append the destination information to the full_state
	# 1. deviation from the destination, measured by rad, ∈[0, π]
	deviation = norm_angle(env.ego_agent.dyn.heading -
							norm_angle(atan2(env.ego_agent.y_dest - env.ego_agent.dyn.y,
						 					env.ego_agent.x_dest - env.ego_agent.dyn.x)))# ∈ [0, 2*π)
	if deviation > π
		deviation = deviation - 2 * π
	end # deviation ∈ [-π, π]

	# 2. Current distance to destination
	dist_to_dest = sqrt((env.ego_agent.dyn.y - env.ego_agent.y_dest)^2 +
		(env.ego_agent.dyn.x - env.ego_agent.x_dest)^2)

	# 3. previous distance to destination
	prev_dist_to_dest = env.ego_agent_dist_to_dest
	## update ego_agent_dist_to_dest
	env.ego_agent_dist_to_dest = dist_to_dest

	# update state entries
	env.full_state[end - AUG_STATE_DIM + 1 : end] .= [deviation, dist_to_dest, prev_dist_to_dest]

	return env.full_state
end

"""
	update advisory for env_agents following low_fi_policy (VICAS)
"""
function get_env_agent_advisory!(ac::Aircraft, env::Simple_CASIM_Env)
	qvals = get_qvals(env.low_fi_policy, convert_ac_state(ac.state))
	if env.stochastic_env_policy == false
		# definitive
		q_ind = indmax(qvals)
	else
		# softmax
		qvals -= maximum(qvals) # preventing extreme values for power
		exp_qvals = exp.(qvals)
		cat_dist = Categorical(exp_qvals / sum(exp_qvals))
		q_ind = rand(cat_dist)
	end
	ac.advisory = env.action_space[q_ind]
end

"""
	 travel towards destination if there is no conflict
"""
function go_to_dest!(ac::Aircraft, env::Simple_CASIM_Env)
	# set threshold for confirming the correct heading towards destination
	threshold = deg2rad(10.)
	heading_vec = [cos(ac.dyn.heading); sin(ac.dyn.heading); 0.]
	desti_dir_angle = norm_angle(atan2(ac.y_dest - ac.dyn.y, ac.x_dest - ac.dyn.x))
	heading_error = norm_angle(abs(desti_dir_angle - ac.dyn.heading))

	if heading_error < threshold
		ac.dyn.turn_rate = 0.
		ac.dyn.heading = desti_dir_angle
	else
		desti_vec = [cos(desti_dir_angle); sin(desti_dir_angle); 0.]
		cross_result = cross(heading_vec, desti_vec)

		if cross_result[end] < 0
			side = :Right
		elseif cross_result[end] > 0
			side = :Left
		else
			side = :OnDir
		end

		# get turn rate
		if side == :Right
			ac.dyn.turn_rate = env.advisory_to_action_dict[:SR]
		elseif side == :Left
			ac.dyn.turn_rate = env.advisory_to_action_dict[:SL]
		else # dest on heading
			ac.dyn.turn_rate = 0. # keep heading
		end
	end
end

"""
	update dynamics for env_agents, env_agents follow VICAS by default
"""
function update_agents_dynamics!(env::Simple_CASIM_Env, advisory::Symbol)
	# get advisory for ego_agent
	env.ego_agent.advisory = advisory
	if advisory != :COC # conflict!
		env.ego_agent.dyn.turn_rate = env.advisory_to_action_dict[advisory]
	else # COC
		go_to_dest!(env.ego_agent, env)
	end

	env.ego_agent.dyn.heading =
		norm_angle(env.ego_agent.dyn.heading + env.ego_agent.dyn.turn_rate * env.ego_agent.dyn.dt)
	ego_x_disp = cos(env.ego_agent.dyn.heading) * env.ego_agent.dyn.v * env.ego_agent.dyn.dt
	ego_y_disp = sin(env.ego_agent.dyn.heading) * env.ego_agent.dyn.v * env.ego_agent.dyn.dt
	# update destination coords
	env.ego_agent.x_dest -= ego_x_disp
	env.ego_agent.y_dest -= ego_y_disp

	for ac in env.env_agents
		# update destination coords
		ac.x_dest -= ego_x_disp
		ac.y_dest -= ego_y_disp
		# get advisories for env_agents
		get_env_agent_advisory!(ac, env)
		if ac.advisory != :COC # conflict!
			ac.dyn.turn_rate = env.advisory_to_action_dict[ac.advisory]
		else # COC
			go_to_dest!(ac, env)
		end
		ac.dyn.heading = norm_angle(ac.dyn.heading + ac.dyn.turn_rate * ac.dyn.dt)
		ac_x_disp = cos(ac.dyn.heading) * ac.dyn.v * ac.dyn.dt
		ac_y_disp = sin(ac.dyn.heading) * ac.dyn.v * ac.dyn.dt
		ac.dyn.x += ac_x_disp - ego_x_disp
		ac.dyn.y += ac_y_disp - ego_y_disp
	end
end

"""
	check if an env_agent goes out of the sensing range
"""
function check_within_sensing(ac::Aircraft, env::Simple_CASIM_Env)
	if sqrt((ac.dyn.x - env.ego_agent.dyn.x)^2 + (ac.dyn.y - env.ego_agent.dyn.y)^2) <= SENSING_RANGE
		return true
	else
		return false
	end
end

function check_within_extended_sensing(ac::Aircraft, env::Simple_CASIM_Env)
	if sqrt((ac.dyn.x - env.ego_agent.dyn.x)^2 + (ac.dyn.y - env.ego_agent.dyn.y)^2) <= SENSING_RANGE * RELAX_FACTOR
		return true
	else
		return false
	end
end

"""
	check if an agent reaches destination
"""
function check_dest(ac::Aircraft)
	if sqrt((ac.dyn.x - ac.x_dest)^2 + (ac.dyn.y - ac.y_dest)^2) <= DEST_CRITERION
		return true
	else
		return false
	end
end

"""
	control (maintain) the traffic density in the airspace as constant
"""
function airspace_control!(env::Simple_CASIM_Env)
	# # constant number of_env agents
	# while length(env.env_agents) < env.const_num_env_agents
	# 	ac = init_one_env_agent(env.ego_agent, env.num_sections)
	# 	push!(env.env_agents, ac)
	# end
	# for i = 1:length(env.env_agents)
	# 	if check_within_sensing(env.env_agents[i], env) == false || check_dest(env.env_agents[i]) == true
	# 		deleteat!(env.env_agents, i)
	# 		ac = init_one_env_agent(env.ego_agent, env.num_sections)
	# 		push!(env.env_agents, ac)
	# 	end
	# end

	i = 1
	while i <= length(env.env_agents)
		if check_within_extended_sensing(env.env_agents[i], env) == false || check_dest(env.env_agents[i]) == true
		# if check_within_sensing(env.env_agents[i], env) == false || check_dest(env.env_agents[i]) == true
			deleteat!(env.env_agents, i)
		end
		i += 1
	end

	# Number of env_agents controlled by Poisson distribution parameterized by spawn rates
	# SPAWN_RATE = env.const_num_env_agents / 100
	# num_new_ac = rand(Poisson(SPAWN_RATE))
	# for i = 1 : num_new_ac
	# 	ac = init_one_env_agent(env.ego_agent, env.num_sections)
	# 	push!(env.env_agents, ac)
	# end

	# Number of env_agents controlled by constant number
	const_num_ac = rand(CONST_NUM_AC_SET)
	while length(env.env_agents) < const_num_ac
		ac = init_one_env_agent(env.ego_agent, env.num_sections)
		push!(env.env_agents, ac)
	end

end

"""
	reward function
	obs is not normalized
"""
function reward_function(env::Simple_CASIM_Env, obs::Vector{Float64}, advisory::Symbol)
	reward = 0.0
	num_nmac = 0

	action = rad2deg(env.advisory_to_action_dict[advisory]) # in [deg]

	# 1. penalize conflict
	if advisory == :COC
		action = 0.0
	else
		reward -= PEN_CONFLICT
	end

	# 2. penalize action
	reward -= PEN_ACTION * action^2

	# 3. penalize distance between intruders
	int_reward = 0.0
	for isection in 1:env.num_sections
		rho = obs[(isection - 1) * STATE_DIM + IND_rho]
		if rho < NMAC_RANGE
			num_nmac += 1
			int_reward -= PEN_NMAC
		end
		int_reward -= PEN_CLOSENESS * exp(- (rho - NMAC_RANGE) / NMAC_RANGE)
	end
	# average the int_reward over sections?
	int_reward /= env.num_sections
	reward += int_reward

	# Additional reward apart from those the same as VICAS:

	# 4. reward reaching destination
	if check_dest(env.ego_agent) == true
		reward += REW_DESTINATION
	end

	# Experiment on the following rewards:

	# 5. penalize pointing far from destination
	reward -= rad2deg(abs(obs[IND_deviation])) * PEN_DEVIATION

	# # 6. penalize digression from destination in terms of distance
	reward -= (obs[IND_dist_to_dest] - obs[IND_prev_dist_to_dest]) * PEN_DIGRESSION
	# OR
	# 6. rewarding getting closer to the destination
	# reward += exp(- obs[IND_dist_to_dest] / DEST_RANGE) * PEN_DIGRESSION

    # scale the step reward (scale reward, do not shift reward!)
    # norm_reward = reward / (POSSIBLE_MAX_REW - POSSIBLE_MIN_REW)
    norm_reward = reward / (POSSIBLE_MAX_REW_VICAS - POSSIBLE_MIN_REW_VICAS)
    # println("norm_rew = ", norm_reward)

    # # non-scaled step reward
    # norm_reward = reward

	return norm_reward, num_nmac
end

"""
	normalize the input to the network
	s = [ρ, θ, ϕ, v_own, v_int]
"""
function normalize_obs!(env::Simple_CASIM_Env, obs::Vector{Float64})
	for isection in 0 : env.num_sections - 1
		if obs[isection * STATE_DIM + IND_rho] == TERM_VAR
			obs[isection * STATE_DIM + 1 : isection * STATE_DIM + STATE_DIM] *= 0 # if terminal, set empty
			# obs[isection * STATE_DIM + 1 : isection * STATE_DIM + STATE_DIM] .= 1 # if terminal, set to 1
		else
			obs[isection * STATE_DIM + IND_rho] /= SENSING_RANGE
			obs[isection * STATE_DIM + IND_theta] /= (2 * pi)
			obs[isection * STATE_DIM + IND_phi] /= (2 * pi)
			# obs[isection * STATE_DIM + IND_v_own] = (obs[isection * STATE_DIM + IND_v_own] - V_MIN) / (V_MAX - V_MIN + 1e-12)
			obs[isection * STATE_DIM + IND_v_own] /= V_MAX
			obs[isection * STATE_DIM + IND_v_int] /= V_MAX
		end
	end
	obs[IND_deviation] /= pi
	obs[IND_dist_to_dest] /= DEST_RANGE
	obs[IND_prev_dist_to_dest] /= DEST_RANGE
end

"""
	for deep correction
"""
function normalize_qvals(qvals::Vector{Float64}; epsilon::Float64=1e-12)
    # return qvals / sum(qvals)
    return qvals / sqrt(max(sum(qvals .^ 2), epsilon))
end

"""
	Utility fusion using max-min
	This function gives min_qvals over multiple intruders
	Action is obtained by taking argmax over actions
"""
function lowfi_values(env::Simple_CASIM_Env, s::Array{Float64, 1})
    # s is the normalized full_state, term_state set to empty (zeros)
    q_matrix = Matrix{Float64}(length(env.action_space), 0) # initiate an empty q_matrix
    rhos = [s[(i - 1) * STATE_DIM + IND_rho] for i in 1 : env.num_sections]
    for i in 1 : env.num_sections
    	if rhos[i] != 0 # section not empty, since s is normalized, empty state is set to zeros
    		section_state = s[STATE_DIM * (i - 1) + 1 : STATE_DIM * i]
    		# de-normalize the state (restore to the original scales)
        	section_state[IND_rho]   *= SENSING_RANGE
        	section_state[IND_theta] *= (2 * pi)
        	section_state[IND_phi]   *= (2 * pi)
        	section_state[IND_v_own] *= V_MAX
        	section_state[IND_v_int] *= V_MAX
        	# only appending qvals to q_matrix if the section_state is non-terminal
        	qvals = get_qvals(env.low_fi_policy, section_state)
        	q_matrix = hcat(q_matrix, qvals)
    	end
    end
    if size(q_matrix, 2) > 0
    	min_qvals = minimum(q_matrix, 2)
    else # q_matrix is empty, then return terminal state values
    	# This only happens when there is no intruder.
    	min_qvals = get_qvals(env.low_fi_policy, terminal_state)
    end

	# q_lo = normalize_qvals(vec(min_qvals)) # deprecated method
	q_lo = vec(min_qvals)

    return reshape(q_lo, 1, :) # Array{Float64, 2}, convert to row vector
end

##################################################################################
## DeepRL.jl interface implementation                                           ##
##################################################################################

function Base.reset(env::Simple_CASIM_Env)
	# env.const_num_env_agents = rand(SPAWN_RATE_SET) # reset SPAWN_RATE for the episode
	env.const_num_env_agents = rand(CONST_NUM_AC_SET)
	env.t = 1
	env.ego_agent = init_ego_agent()
	env.env_agents = init_env_agents(env.ego_agent, env.init_num_env_agents, env.num_sections)
	# env.const_num_env_agents = env.init_num_env_agents # + rand(0:5)
	update_states!(env) # update observation
	obs = env.full_state
	normalize_obs!(env, obs)
	return obs
end

function POMDPs.actions(env::Simple_CASIM_Env)
	return env.action_space
end

function DeepRL.sample_action(env::Simple_CASIM_Env)
    # return rand(env.rng, actions(env))
    return rand(actions(env)) # ignoring rng
end

function POMDPs.n_actions(env::Simple_CASIM_Env)
	return length(env.action_space)
end

function DeepRL.obs_dimensions(env::Simple_CASIM_Env)
	# obs = full_state
	return size(env.full_state)
end

function POMDPs.discount(problem::Dummy_Problem)
	return problem.discount_factor
end

function POMDPs.action_index(problem::Dummy_Problem, action::Symbol)
	# Symbol[:SR, :WR, :KEEP, :WL, :SL, :COC]
	if action == :SR
		return 1
	elseif action == :WR
		return 2
	elseif action == :KEEP
		return 3
	elseif action == :WL
		return 4
	elseif action == :SL
		return 5
	else
		return 6
	end
end

function DeepRL.step!(env::Simple_CASIM_Env, advisory::Symbol)
	update_agents_dynamics!(env, advisory)
	done = check_dest(env.ego_agent)
	obs = update_states!(env) # update observation, non-normalized (for getting reward)
	# println(obs)
	if env.render
    	render_env(env)
    end
	reward, num_nmac = reward_function(env, obs, advisory)
	# reward is scaled to fit the output of lowfi_values and correction network
	env.t += 1
    info = env.t
	if !done
    	airspace_control!(env)
    end
    normalize_obs!(env, obs)
    # the returned obs is normalized, therefore, the states in the replay buffer is normalized
	return obs, reward, num_nmac, done, info
end

end # module

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
module ZygoteFFTs
using Zygote
using AbstractFFTs
using FillArrays

import Zygote:@adjoint

AbstractFFTs.bfft(x::Fill, dims...) = AbstractFFTs.bfft(collect(x), dims...)
AbstractFFTs.rfft(x::Fill, dims...) = AbstractFFTs.rfft(collect(x), dims...)
AbstractFFTs.irfft(x::Fill, d, dims...) = AbstractFFTs.irfft(collect(x), d, dims...)

@adjoint function fft(xs)
    return AbstractFFTs.fft(xs), function(Δ)
        return (AbstractFFTs.bfft(Δ),)
    end
end

# all of the plans normalize their inverse, while we need the unnormalized one.
@adjoint function *(P::AbstractFFTs.Plan, xs)
    return P * xs, function(Δ)
        N = prod(size(xs)[[P.region...]])
        return (nothing, N * (P \ Δ))
    end
end

@adjoint function \(P::AbstractFFTs.Plan, xs)
    return P \ xs, function(Δ)
        N = length(Δ)
        return (nothing, 1/N * (P * Δ))
    end
end

@adjoint function ifft(xs)
    return AbstractFFTs.ifft(xs), function(Δ)
        N = length(xs)
        return (1/N* AbstractFFTs.fft(Δ),)
    end
end

@adjoint function bfft(xs)
    return AbstractFFTs.bfft(xs), function(Δ)
        return (AbstractFFTs.fft(Δ),)
    end
end

# to actually use rfft, one needs to insure that everything that happens in the
# Fourier domain could've been done in the space domain with real numbers. This
# means enforcing conjugate symmetry along all transformed dimensions besides
# the first. Otherwise this is going to result in *very* weird behavior.
@adjoint function rfft(xs::AbstractArray{<:Real})
    return AbstractFFTs.rfft(xs), function(Δ)
        N = length(Δ)
        originalSize = size(xs,1)
        return (AbstractFFTs.brfft(Δ, originalSize),)
    end
end

@adjoint function irfft(xs, d)
    return AbstractFFTs.irfft(xs, d), function(Δ)
        total = length(Δ)
        fullTransform = 1/total * AbstractFFTs.rfft(real.(Δ))
        return (fullTransform, nothing)
    end
end

@adjoint function brfft(xs, d)
    return AbstractFFTs.brfft(xs, d), function(Δ)
        fullTransform = AbstractFFTs.rfft(real.(Δ))
        return (fullTransform, nothing)
    end
end

@adjoint function fftshift(x)
    return fftshift(x), function(Δ)
        return (ifftshift(Δ),)
    end
end

@adjoint function ifftshift(x)
    return ifftshift(x), function(Δ)
        return (fftshift(Δ),)
    end
end

# if we're specifying the dimensions
@adjoint function fft(xs, dims)
    return AbstractFFTs.fft(xs, dims), function(Δ)
        # dims can be int, array or tuple,
        # convert to collection for use as index
        dims = collect(dims)
        return (AbstractFFTs.bfft(Δ, dims), nothing)
    end
end

@adjoint function bfft(xs, dims)
    return AbstractFFTs.ifft(xs, dims), function(Δ)
        dims = collect(dims)
        return (AbstractFFTs.fft(Δ, dims),nothing)
    end
end

@adjoint function ifft(xs, dims)
    return AbstractFFTs.ifft(xs, dims), function(Δ)
        dims = collect(dims)
        N = prod(collect(size(xs))[dims])
        return (1/N * AbstractFFTs.fft(Δ, dims),nothing)
    end
end

@adjoint function rfft(xs, dims)
    return AbstractFFTs.rfft(xs, dims), function(Δ)
        dims = collect(dims)
        N = prod(collect(size(xs))[dims])
        return (N * AbstractFFTs.irfft(Δ, size(xs,dims[1]), dims), nothing)
    end
end

@adjoint function irfft(xs, d, dims)
    return AbstractFFTs.ifft(xs, dims), function(Δ)
        dims = collect(dims)
        N = prod(collect(size(xs))[dims])
        return (1/N * AbstractFFTs.rfft(Δ, dims), nothing, nothing)
    end
end
@adjoint function brfft(xs, d, dims)
    return AbstractFFTs.ifft(xs, dims), function(Δ)
        dims = collect(dims)
        return (AbstractFFTs.rfft(Δ, dims), nothing, nothing)
    end
end

@adjoint function fftshift(x, dims)
    return fftshift(x), function(Δ)
        return (ifftshift(Δ, dims), nothing)
    end
end

@adjoint function ifftshift(x, dims)
    return ifftshift(x), function(Δ)
        return (fftshift(Δ, dims), nothing)
    end
end

end # module


text:
module ImageAnalysis

using FileIO, ColorTypes, FFTW

include("IO.jl")
include("utils.jl")
include("stats.jl")
include("histogram.jl")
include("integralImages.jl")
include("FFTconvolutions.jl")
include("kernels.jl")
include("thresholding.jl")
include("morphology.jl")
include("CCL.jl")
include("convexHullSIMP.jl")
include("entropy.jl")
include("moments.jl")
include("integral_array_utils.jl")

greet() = print("Hello World!")

end # module

max_stars_repo_name: BrucePBon/ImageAnalysis

max_issues_repo_name: BrucePBon/ImageAnalysis

max_forks_repo_name: BrucePBon/ImageAnalysis


text:
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# Last edited: 30 December, 2014
# ===============================

# =====================================================================
# 2014-12-30: - cell variables can be created including bounddary cells
# =====================================================================

# ============================================================

"""
creates a cell variable over the domain m and assigns value phi0 to it
function createCellVariable(m::MeshStructure, phi0::Real)
"""
function createCellVariable(m::MeshStructure, phi0::Real)
# creates a cell variable and assigns value phi0 to it
CellValue(m, phi0*ones(tuple(m.dims.+2...)))
end

function createCellVariable(phi::CellValue)
  # similar to copyCell
  CellValue(phi.domain, Base.copy(phi.value))
end

# ============================================================
function createCellVariable(m::MeshStructure, phi0::Array{T}) where T<:Real
# creates a cell variable and assigns value phi0 to it
if prod(m.dims.+2)==length(phi0)
  CellValue(m, phi0)
elseif prod(m.dims)==length(phi0)
  d=m.dimension
  phival = zeros(tuple(m.dims.+2...))
  BC = createBC(m) # Neumann boundaries
  if d==1 || d==1.5
    phival[2:end-1] = phi0
  elseif d==2 d==2 || d==2.5 || d==2.8
    phival[2:end-1, 2:end-1] = phi0
  elseif d==3 || d==3.2
    phival[2:end-1,2:end-1,2:end-1] = phi0
  end
  phi = CellValue(m, phival)
  cellBoundary!(phi, BC)
else
  error("JFVM: Matrix must be the same size as the domain.")
end
end

# ==============================================================
function createCellVariable(m::MeshStructure, phi0::Array{T}, BC::BoundaryCondition) where T<:Real
# creates a cell variable and assigns value phi0 to it
if prod(m.dims.+2)==length(phi0)
  error("JFVM: Matrix must be the same size as the domain.")
elseif prod(m.dims)==length(phi0)
  d=m.dimension
  phival = zeros(tuple(m.dims.+2...))
  if d==1 || d==1.5
    phival[2:end-1] = phi0
  elseif d==2 || d==2.5 || d==2.8
    phival[2:end-1, 2:end-1] = phi0
  elseif d==3 || d==3.2
    phival[2:end-1,2:end-1,2:end-1] = phi0
  end
  phi = CellValue(m, phival)
  cellBoundary!(phi, BC)
else
  error("JFVM: Matrix must be the same size as the domain.")
end
end

# ============================================================
function createCellVariable(m::MeshStructure, phi0::Real, BC::BoundaryCondition)
phi = CellValue(m, phi0*ones(tuple(m.dims.+2...)))
cellBoundary!(phi, BC)
end

# ============================================================
function createFaceVariable(m::MeshStructure, phi0::Array{T,1}) where T<:Real
# creates a face variable based on the mesh structure
d=m.dimension
  if d==1 || d==1.5
    FaceValue(m,
	      ones(m.dims[1]+1)*phi0[1],
	      [1.0],
	      [1.0])
  elseif d==2 || d==2.5 || d==2.8
    FaceValue(m,
	      ones(m.dims[1]+1, m.dims[2])*phi0[1],
	      ones(m.dims[1], m.dims[2]+1)*phi0[2],
	      [1.0])
  elseif d==3 || d==3.2
    FaceValue(m,
	      ones(m.dims[1]+1, m.dims[2], m.dims[3])*phi0[1],
	      ones(m.dims[1], m.dims[2]+1, m.dims[3])*phi0[2],
	      ones(m.dims[1], m.dims[2], m.dims[3]+1)*phi0[3])
  end
end

function createFaceVariable(m::MeshStructure, phi0::Real)
# creates a face variable based on the mesh structure
d=m.dimension
  if d==1 || d==1.5
    FaceValue(m,
	      ones(m.dims[1]+1)*phi0,
	      [1.0],
	      [1.0])
  elseif d==2 || d==2.5 || d==2.8
    FaceValue(m,
	      ones(m.dims[1]+1, m.dims[2])*phi0,
	      ones(m.dims[1], m.dims[2]+1)*phi0,
	      [1.0])
  elseif d==3 || d==3.2
    FaceValue(m,
	      ones(m.dims[1]+1, m.dims[2], m.dims[3])*phi0,
	      ones(m.dims[1], m.dims[2]+1, m.dims[3])*phi0,
	      ones(m.dims[1], m.dims[2], m.dims[3]+1)*phi0)
  end
end

# ============================================================
function createCellVector(m::MeshStructure, phi0::Array{T,1}) where T<:Real
# creates a cell vector based on the mesh structure
d=m.dimension
  if d==1 || d==1.5
    CellVector(m,
	      ones(m.dims[1])*phi0[1],
	      [1.0],
	      [1.0])
  elseif d==2 || d==2.5 || d==2.8
    CellVector(m,
	      ones(m.dims[1], m.dims[2])*phi0[1],
	      ones(m.dims[1], m.dims[2])*phi0[2],
	      [1.0])
  elseif d==3 || d==3.2
    CellVector(m,
	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0[1],
	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0[2],
	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0[3])
  end
end

function createCellVector(m::MeshStructure, phi0::Real)
# creates a cell vector based on the mesh structure
d=m.dimension
  if d==1 || d==1.5
    CellVector(m,
	      ones(m.dims[1])*phi0,
	      [1.0],
	      [1.0])
  elseif d==2 || d==2.5 || d==2.8
    CellVector(m,
	      ones(m.dims[1], m.dims[2])*phi0,
	      ones(m.dims[1], m.dims[2])*phi0,
	      [1.0])
  elseif d==3 || d==3.2
    CellVector(m,
	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0,
	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0,
	      ones(m.dims[1], m.dims[2], m.dims[3])*phi0)
  end
end

# ================== copy function for cell variables ===================
function copyCell(phi::CellValue)
  CellValue(phi.domain, Base.copy(phi.value))
end


text:
#__precompile__()

module EAGODomainReduction

using MathProgBase
using EAGOSmoothMcCormickGrad
using IntervalArithmetic
using IntervalContractors

# import basic functions to overload
import Base: +, -, /, *, ^, exp, exp2, exp10, log, log2, log10, sign, step, abs,
             sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh,
             atanh, one, zero, sqrt, min, max, convert, promote_rule

# replace with using once IntervalContractor updated
import IntervalContractors: plus_rev, minus_rev, mul_rev, power_rev, sqr_rev,
                            sqrt_rev, abs_rev, sin_rev, cos_rev, tan_rev,
                            asin_rev, log_rev, exp_rev, inv_rev

# exports main FBBT functions
export Generate_Tape, Generate_TapeList, Generate_Fixed_Tape, Generate_Fixed_TapeList,
       FFBT_Refine, SetConstraint!, SetVarBounds!, GetVarBounds, ForwardPass!,
       ReversePass!, DAGContractor!, NodeFinder, Tape, TapeList, getDAG

# exports supplemental reverse contractors (push to IntervalContractor Library)
export tanh_rev, div_revDR, acos_rev, atan_rev, sinh_rev, cosh_rev, tanh_rev,
       asinh_rev, acosh_rev, atanh_rev, min_rev, max_rev, step_rev, sign_rev,
       exp2_rev, exp10_rev, log2_rev, log10_rev, one_rev, zero_rev

# exports OBBT functions
export Variable_DR!, Variable_DR_Imp!,
       STD_Linear_RR!,
       STD_LP_Probe! #LP_contractor, poorLP_contractor

# includes utility subroutines
include("src/utils.jl")

# includes subroutines for optimality-based bound tightening
include("src/OBBT/Variable_Dual.jl")
include("src/OBBT/Linear_RR.jl")
include("src/OBBT/STD_LP_Probe.jl")
include("src/OBBT/poorLP.jl")

# includes subroutines for feasibility-based bound tightening
include("src/FBBT/NodeFinder.jl")
include("src/FBBT/DAGprop.jl")
include("src/FBBT/Utilities.jl")

#=
function __init__()
end
=#

end # module


text:
using Pkg

for p in ("ArgParse", "Knet", "AutoGrad", "Gym")
    if !haskey(Pkg.installed(),p)
        Pkg.add(p)
        if p == "Gym"
            ENV["GYM_ENVS"] = "atari:algorithmic:box2d:classic_control"
            Pkg.build("Gym")
        end
    end
end

"""
julia reinforce_discrete.jl

This example implements the REINFORCE algorithm from
`Simple statistical gradient-following algorithms for
connectionist reinforcement learning.`,  Williams, Ronald J.
Machine learning, 8(3-4):229–256, 1992. This example also
demonstrates the usage of the `@zerograd` function for
stopping the gradient flow.
"""
module REINFORCE_DISCRETE

using Gym, ArgParse, Knet, AutoGrad, Statistics

function predict_linear(w, ob)
    linear = w["w"] * ob .+ w["b"]
    return linear
end

function sample_action(linear)
    linear = Array(linear)
    probs = vec(exp.(linear) ./ sum(exp.(linear); dims=1))
    c_probs = cumsum(probs)
    return findmax(c_probs .> rand())[2]
end

AutoGrad.@zerograd sample_action(linear)

function play(w, ob)
    linear = predict_linear(w, ob)
    action = sample_action(linear)
    return action, linear
end

function play_episode(w, env, o)
    ob = reset!(env)
    rewards = Float32[]
    linears = Any[]
    actions = Array{Int,1}()
    total = 0

    for t=1:env.spec.max_episode_steps
        ob_inp = convert(o["atype"], reshape(ob, size(ob, 1), 1))
        action, linear = play(w, ob_inp)
        push!(linears, linear)
        ob, reward, done, _ = step!(env, action-1)

        total += reward

        push!(rewards, reward)
        push!(actions, action)

        o["render"] && render(env)

        done && break
    end
    return linears, actions, rewards, total
end

function loss(w, env, o; totalr=nothing)
    linears, actions, rewards, total = play_episode(w, env, o)
    totalr[1] = total

    y = zeros(size(w["w"], 1), length(actions))

    for i=1:length(actions)
        y[actions[i], i] = 1.0
    end

    actions = convert(o["atype"], y)

    discounted = discount(rewards; γ=o["gamma"])
    discounted = discounted .- mean(discounted)#mean R as a baseline

    discounted = convert(o["atype"], reshape(discounted, 1, size(actions, 2)))

    logps = sum(logp(hcat(linears...); dims=1) .* actions; dims=1)
    -sum(logps .* discounted) ./ size(actions, 2)
end

lossgradient = grad(loss)

function init_weights(input, output, atype)
    w = Dict()
    w["w"] = xavier(output, input)
    w["b"] = zeros(output, 1)

    for k in keys(w)
        w[k] = convert(atype, w[k])
    end

    return w
end

function discount(rewards; γ=0.9)
    discounted = zeros(Float32, length(rewards), 1)
    discounted[end] = rewards[end]

    for i=(length(rewards)-1):-1:1
        discounted[i] = rewards[i] + γ * discounted[i+1]
    end
    return discounted
end

function train!(w, opts, env, o)
    totalr = [0.0]
    g = lossgradient(w, env, o; totalr=totalr)
    update!(w, g, opts)
    return totalr[1]
end

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="(c) Ozan Arkan Can, 2018. Demonstration of the REINFORCE algorithm on the discrete action space."
    @add_arg_table s begin
        ("--env_id"; default="CartPole-v1"; help="environment name")
        ("--episodes"; arg_type=Int; default=20; help="number of episodes")
        ("--gamma"; arg_type=Float64; default=0.99; help="doscount factor")
        ("--threshold"; arg_type=Int; default=1000; help="stop the episode even it is not terminal after number of steps exceeds the threshold")
        ("--lr"; arg_type=Float64; default=0.01; help="learning rate")
        ("--render"; help = "render the environment"; action = :store_true)
        ("--usegpu"; action=:store_true; help="use GPU or not")
    end

    Knet.seed!(12345)
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end

    o = parse_args(args, s)
    o["atype"] = !o["usegpu"] ? Array{Float32} : KnetArray{Float32}

    env = GymEnv(o["env_id"])
    seed!(env, 12345)

    INPUT = env.observation_space.shape[1]
    OUTPUT = env.action_space.n

    w = init_weights(INPUT, OUTPUT, o["atype"])
    opts = Dict()
    for k in keys(w)
        opts[k] = Rmsprop(lr=o["lr"])
    end

    for i=1:o["episodes"]
        total = train!(w, opts, env, o)
        println("episode $i , total rewards: $total")
    end
end

PROGRAM_FILE=="reinforce_discrete.jl" && main(ARGS)

end


text:
using Base.Test

using Celeste: Model, Transform, SensitiveFloats, MCMC, Synthetic
using StatsBase

include(joinpath(Pkg.dir("Celeste"), "test", "SampleData.jl"))

# helper to create synthetic data
function generate_single_star_data(; lnr=7.7251)
    # generate images (data)
    ea, vp, catalog = SampleData.gen_sample_star_dataset()

    # adjust fluxes to match input lnr
    old_fluxes = catalog[1].star_fluxes
    _, colors = Model.fluxes_to_colors(old_fluxes)
    new_fluxes = Model.colors_to_fluxes(lnr, colors)
    catalog[1].star_fluxes = new_fluxes

    # cache fluxes location, etc
    ce0    = deepcopy(catalog[1])
    images = deepcopy(ea.images)

    # make mod/dat images
    dat_images = Synthetic.gen_blob(images, catalog; expectation=false);

    # make sure ea images data reflects dat images for vb fitting
    for b in 1:5
      ea.images[b].pixels[:] = dat_images[b].pixels[:]
    end

    # save truedf for python plot comparison
    truedf = MCMC.catalog_to_data_frame_row(ce0)
    return truedf, ce0, dat_images, ea, vp
end

function generate_params()
    lnfluxes = MCMC.sample_logfluxes(; is_star=true, lnr=nothing)
    lu       = .01 * randn(2)
    return vcat([lnfluxes, lu]...)
end

###############
#### tests ####
###############

function test_mcmc_catalog_to_data_frame_row()
    # create slightly less bright log re
    truedf, ce0, dat_images, ea, vp = generate_single_star_data(; lnr=6.)
    @test true
end

function test_star_loglike()
    # gen data and initial params
    truedf, ce0, dat_images, ea, vp = generate_single_star_data(; lnr=6.)

    # generate star params
    lnfluxes = MCMC.sample_logfluxes(; is_star=true, lnr=nothing)
    lu       = .01 * randn(2)
    th = vcat([lnfluxes, lu]...)

    # create loglike
    init_pos = deepcopy(vp[1][ids.pos])
    star_loglike, constrain_pos, unconstrain_pos =
        MCMC.make_star_loglike(dat_images, init_pos)

    ll = star_loglike(th)
    @test !isnan(ll)

    # make sure constriain/unconstrain work
    uu = constrain_pos(lu)
    lu2 = unconstrain_pos(uu)
    @test isapprox(lu, lu2)

end

function test_logflux_logprior()
    # test sample (with all args)
    lnfluxes = MCMC.sample_logfluxes(; is_star=true, lnr=nothing)
    lnfluxes = MCMC.sample_logfluxes(; is_star=false, lnr=nothing)
    lnfluxes = MCMC.sample_logfluxes(; is_star=true, lnr=5.)
    lnfluxes = MCMC.sample_logfluxes(; is_star=false, lnr=5.)

    ll = MCMC.logflux_logprior(lnfluxes; is_star=true)
    @test !isnan(ll)

    ll = MCMC.logflux_logprior(lnfluxes; is_star=false)
    @test !isnan(ll)
end

function test_slicesample()
    function lnpdf(th)
        return -1*sum(th.*th)
    end

    th = randn(5)
    chain, lls = MCMC.slicesample_chain(lnpdf, th, 10; print_skip=20)
    @test true
end

println("Running mcmc tests")
test_mcmc_catalog_to_data_frame_row()
test_star_loglike()
test_logflux_logprior()
test_slicesample()


text:
"""
Compute the spectral tendency of the surface pressure logarithm
"""
function surface_pressure_tendency!(Prog::PrognosticVariables{NF}, # Prognostic variables
                                   Diag::PrognosticVariables{NF}, # Diagnostic variables
                                   l2::Int,                       # leapfrog index 2 (time step used for tendencies)
                                   C::Constants{NF}
                                   ) where {NF<:AbstractFloat}

    @unpack pres_surf= Prog
    @unpack pres_surf_tend = Diag.Tendencies

    @unpack u_grid,v_grid,div_grid = Diag.gridvars
    @unpack u_mean,v_mean,d_mean,dumc,px,py = Diag.miscvars
    @unpack dhs = C

    _,_,nlev = size(u_grid)

    #Initialise mean fields.
    fill!(u_mean,0)
    fill!(v_mean,0)
    fill!(d_mean,0) #d-mean is calculated here but actually used in vertical_velocity_tendency!()

    #..and calculate values
    for k in 1:nlev
        u_mean += u_grid[:,:,k]*dhs[k]
        v_mean += v_grid[:,:,k]*dhs[k]
        d_mean += div_grid[:,:,k]*dhs[k]
    end

    #Now use the mean fields
    grad!(pres_surf[:,:,l2], dumc[:,:,2], dumc[:,:,3])
    px = gridded(dumc[:,:,2]*3600, scale=true)
    py = gridded(dumc[:,:,3]*3600, scale=true) #3600 factor from Paxton/Chantry. I think this is to correct for the underflow rescaling earlier

    pres_surf_tend = spectral(-u_mean.*px - v_mean.*py)
    pres_surf_tend[1,1] = Complex{NF}(zero)

end

"""
Compute the spectral tendency of the "vertical" velocity
"""
function vertical_velocity_tendency!(Prog::PrognosticVariables{NF}, # Prognostic variables
                                     Diag::PrognosticVariables{NF}, # Diagnostic variables
                                     C::Constants{NF}
                                     ) where {NF<:AbstractFloat}

    @unpack u_grid,v_grid,div_grid = Diag.gridvars
    @unpack u_mean,v_mean,d_mean,sigma_tend,sigma_m,puv,px,py = Diag.miscvars
    @unpack dhs = C

    _,_,nlev = size(u_grid)

    #Initialise
    sigma_tend[:,:,:] = 0.0
    sigma_m[:,:,:] = 0.0

    for k in 1:nlev
        puv[:,:,k] = (u_grid[:,:,k] - u_mean)*px + (v_grid[:,:,k] - v_mean)*py
    end

    for k in 1:nlev
        sigma_tend[:,:,k+1] = sigma_tend[:,:,k] - dhs[k]*(puv[:,:,k] + div_grid[:,:,k] - d_mean)
        sigma_m[:,:,k+1] = sigma_m[:,:,k] - dhs[k]*puv[:,:,k]
    end

end

"""
Compute the temperature anomaly in grid point space
"""
function temperature_grid_anomaly!(Diag::PrognosticVariables{NF}, # Diagnostic variables
                                   C::Constants{NF}
                                   ) where {NF<:AbstractFloat}

    @unpack temp_grid,temp_grid_anomaly = Diag.gridvars
    @unpack temp_ref = C

    _,_,nlev = size(temp_grid)

    for k in 1:nlev
        temp_grid_anomaly[:,:,k] = temp_grid[:,:,k] .- temp_ref[k]
    end

end

"""
Compute the spectral tendency of the zonal wind
"""
function zonal_wind_tendency!(Diag::PrognosticVariables{NF}, # Diagnostic variables
                              C::Constants{NF}
                              )where {NF<:AbstractFloat}

    @unpack u_tend = Diag.Tendencies
    @unpack u_grid,v_grid,vor_grid,temp_grid_anomaly= Diag.gridvars
    @unpack sigma_tend,px,py,arbitrary_array = Diag.miscvars
    @unpack rgas,dhsr = C

    _,_,nlev = size(u_grid)

    #Update px,py
    px = rgas*px
    py = rgas*py

    #Initialise arbitrary_array
    fill!(arbitrary_array,0)

    for k in 2:nlev
        arbitrary_array[:,:,k] = sigma_tend[:,:,k].*(u_grid[:,:,k] - u_grid[:,:,k-1])
    end

    for k in 1:nlev
        u_tend[:,:,k] = u_tend[:,:,k]
                        + v_grid[:,:,k].*vor_grid[:,:,k]
                        - temp_grid_anomaly[:,:,k]*px
                         - (arbitrary_array[:,:,k+1] + arbitrary_array[:,:,k])*dhsr[k]
    end

end

"""
Compute the spectral tendency of the meridional wind
"""
function meridional_wind_tendency!(Diag::PrognosticVariables{NF}, # Diagnostic variables
                                   C::Constants{NF}
                                  )where {NF<:AbstractFloat}

    @unpack v_tend = Diag.Tendencies
    @unpack vor_grid,u_grid,v_grid,temp_grid_anomaly =Diag.gridvars
    @unpack sigma_tend,arbitrary_array,px,py = Diag.miscvars
    @unpack dhsr= C

    _,_,nlev = size(u_grid)

    for k in 2:nlev
        arbitrary_array[:,:,k] = sigma_tend[:,:,k].*(v_grid[:,:,k] - v_grid[:,:,k-1])
    end

    for k in 1:nlev
        v_tend[:,:,k] = v_tend[:,:,k]
                        -u_grid[:,:,k].*vor_grid[:,:,k]
                        - temp_grid_anomaly[:,:,k]*py
                        - (arbitrary_array[:,:,k+1] + arbitrary_array[:,:,k])*dhsr[k]
    end

end

"""
Compute the spectral temperature tendency
"""
function temperature_tendency!(Diag::PrognosticVariables{NF}, # Diagnostic variables
                               C::Constants{NF}
                               )where {NF<:AbstractFloat}

    @unpack temp_tend = Diag.Tendencies
    @unpack div_grid,temp_grid_anomaly =Diag.gridvars
    @unpack d_mean,arbitrary_array,sigma_tend,sigma_m,puv= Diag.miscvars
    @unpack temp_ref,tref3,dhsr,fsgr,akap = C

    _,_,nlev = size(div_grid)

    for k in 2:nlev
        arbitrary_array[:,:,k] = sigma_tend[:,:,k].*(temp_grid_anomaly[:,:,k] - temp_grid_anomaly[:,:,k-1])
                    + sigma_m[:,:,k].*(temp_ref[k] - temp_ref[k-1])
    end

    for k in 1:nlev
        temp_tend[:,:,k] = temp_tend[:,:,k]
                        + temp_grid_anomaly[:,:,k].*div_grid[:,:,k]
                        - (arbitrary_array[:,:,k+1] + arbitrary_array[:,:,k])*dhsr[k]
                        + fsgr[k]*temp_grid_anomaly[:,:,k].*(sigma_tend[:,:,k+1] + sigma_tend[:,:,k])
                        + tref3[k]*(sigma_m[:,:,k+1] + sigma_m[:,:,k])
                        + akap*(t_grid[:,:,k].*puv[:,:,k] - temp_grid_anomaly[:,:,k].*d_mean)
    end

end

"""
Compute the humidity tendency
"""
function humidity_tendency!(Diag::PrognosticVariables{NF}, # Diagnostic variables
                          C::Constants{NF}
                          )where {NF<:AbstractFloat}

    @unpack div_grid,humid_grid = Diag.gridvars
    @unpack arbitrary_array,sigma_tend= Diag.miscvars
    @unpack dhsr = C

    _,_,nlev = size(div_grid)

    for k in 2:nlev
        arbitrary_array[:,:,k] = sigma_tend[:,:,k].*(humid_grid[:,:,k,itr] - humid_grid[:,:,k-1,itr])
    end

    arbitrary_array[:,:,2:3] .= zero(NF)

    for k in 1:nlev
        humid_tend[:,:,k,itr] = humid_tend + humid_grid[:,:,k,itr].*div_grid[:,:,k]
                - (arbitrary_array[:,:,k+1] + arbitrary_array[:,:,k])*dhsr[k]
        end

end