text:
# Ensure that Logging.jl works when loaded with `import` rather than `using`
# since `import` does not bring constants like `DEBUG` into scope:
module InnerModule

import Logging
@Logging.configure(level=DEBUG)

using Compat
using Base.Test

function test_non_global_interpolation(y)
    @info("y = $y")
end

test_non_global_interpolation(5)

function test_log_macro_common(flags)
    for (macroname, isshown) in flags
        ex = Expr(:macrocall, Symbol(macroname), "test message")
        if isshown
            @test ex |> macroexpand != :nothing
        else
            @test ex |> macroexpand == :nothing
        end
    end
end

test_log_macro_common([(:@debug, true), (:@info, true), (:@warn, true),
    (:@err, true), (:@critical, true)])

@Logging.configure(level=INFO)
test_log_macro_common([(:@debug, false), (:@info, true), (:@warn, true),
    (:@err, true), (:@critical, true)])

@Logging.configure(level=WARNING)
test_log_macro_common([(:@debug, false), (:@info, false), (:@warn, true),
    (:@err, true), (:@critical, true)])

@Logging.configure(level=ERROR)
test_log_macro_common([(:@debug, false), (:@info, false), (:@warn, false),
    (:@err, true), (:@critical, true)])

@Logging.configure(level=CRITICAL)
test_log_macro_common([(:@debug, false), (:@info, false), (:@warn, false),
    (:@err, false), (:@critical, true)])

@Logging.configure(level=OFF)
test_log_macro_common([(:@debug, false), (:@info, false), (:@warn, false),
    (:@err, false), (:@critical, false)])

end


text:
using ClimateDataIO
using Test
using Dates
# using DataFrames

# Last Edit: 29.07.2021

# Define 7zip location constant
if Sys.iswindows() & (VERSION >= VersionNumber("1.3.0"))
	const exe7z = joinpath(splitdir(Sys.BINDIR)[1],"libexec","7z.exe")
elseif Sys.iswindows() & (VERSION < VersionNumber("1.3.0"))
	const exe7z = joinpath(splitdir(Sys.BINDIR)[1],"7z.exe")
elseif Sys.isunix() & (VERSION >= VersionNumber("1.3.0"))
	const exe7z = joinpath(splitdir(Sys.BINDIR)[1],"libexec","7z")
elseif Sys.isunix() & (VERSION < VersionNumber("1.3.0"))
	const exe7z = joinpath(splitdir(Sys.BINDIR)[1],"7z")
end

# TODO: Add tests for GHG_LOAD and GHG_READ Biomet loading scenarios

@testset "============================== SLT TESTS  ==============================" begin
	# SLTLOAD: Check a known set of data
	println("\n====  SLTLOAD  ====")
	src = splitdir(@__FILE__)[1]
	mindate = DateTime(2016,11,23,14,4)
	maxdate = DateTime(2016,11,23,18)
	Data = ClimateDataIO.slt_load(src,mindate,maxdate,verbose=true)
	tmean, Dmean, Dstd, Dmin, Dmax = ClimateDataIO.slt_load(src,mindate,maxdate,verbose=true,average=true)
	@test size(Data) == (141377,13) || "SLTLOAD: Output size incorrect, should be (141377,13)"
	@test Data.Time[1] == DateTime(2016,11,23,14,4) || "SLTLOAD: First timestamp should be 2016-11-23T14:04:00"
	@test Data.Time[end] == DateTime(2016,11,23,17,59,59,900) || "SLTLOAD: Last timestamp should be 2016-11-23T17:59:59.9"

	tmean, Dmean, Dstd, Dmin, Dmax = ClimateDataIO.slt_load(src,mindate,maxdate,verbose=true,average=true)
	@test size(Dmean) == (8,12) || "SLTLOAD: Output size incorrect, should be (8,12)"
	@test size(Dstd) == (8,12) || "SLTLOAD: Output size incorrect, should be (8,12)"
	@test size(Dmin) == (8,12) || "SLTLOAD: Output size incorrect, should be (8,12)"
	@test size(Dmax) == (8,12) || "SLTLOAD: Output size incorrect, should be (8,12)"
	@test length(tmean) == 8 || "SLTLOAD: Output length incorrect, should be 8"
	@test tmean[1] == DateTime(2016,11,23,14,4) || "SLTLOAD: First timestamp should be 2016-11-23T14:04:00"
	@test tmean[end] == DateTime(2016,11,23,17,34) || "SLTLOAD: Last timestamp should be 2016-11-23T17:34:00"

	src = splitdir(@__FILE__)[1]
	mindate = DateTime(2016,11,23,14,4)
	maxdate = DateTime(2016,11,23,18)
	dest = joinpath(splitdir(@__FILE__)[1],"temporary_files")
	isdir(dest) ? nothing : mkdir(dest)
	maxcols = 2 # Max analog columns
	ClimateDataIO.slt_trim(src,dest,mindate,maxdate,maxcols)
	Data = ClimateDataIO.slt_load(dest,mindate,maxdate)
	f = ["W20163281404.cfg",
		"W20163281404.csr",
		"W20163281404.csu",
		"W20163281404.csv",
		"W20163281404.flx",
		"W20163281404.log",
		"W20163281404.slt",
		"W20163281430.slt",
		"W20163281500.slt",
		"W20163281530.slt",
		"W20163281600.slt",
		"W20163281630.slt",
		"W20163281700.slt",
		"W20163281730.slt"]
	for i in f
		rm(joinpath(dest,i))
	end
	@test size(Data) == (141377,9) || "SLTLOAD: Output size incorrect, should be (141377,9)"
	@test Data.Time[1] == DateTime(2016,11,23,14,4) || "SLTLOAD: First timestamp should be 2016-11-23T14:04:00"
	@test Data.Time[end] == DateTime(2016,11,23,17,59,59,900) || "SLTLOAD: Last timestamp should be 2016-11-23T17:59:59.9"

	src = joinpath(splitdir(@__FILE__)[1],"W20163280930.slt")
	header = ClimateDataIO.slt_header(src,6,10)
	println("SLTLOAD Test Complete")

	# SLTREAD: Load one known file
	println("\n====  SLTREAD Tests  ====")
	src = splitdir(@__FILE__)[1]
	src = joinpath(src,"W20163281000.slt")
	analog_inputs = 6 # Number of analog inputs
	sample_frequency = 10 # Hz
	Time,Data = ClimateDataIO.slt_read(src,analog_inputs,sample_frequency)
	@test size(Data) == (18003,10) || "SLTREAD: Output size incorrect, should be (18003,10)"
	@test Time[1] == DateTime(2016,11,23,10) || "SLTREAD: First timestamp should be 2016-11-23T10:00:00"
	@test Time[end] == DateTime(2016,11,23,10,30,00,200) || "SLTREAD: Last timestamp should be 2016-11-23T10:30:00.2"
	println("SLTREAD Test Complete")

	# SLT_TIMESHIFT
	println("\n====  SLT_TIMESHIFT Tests  ====")
	src = splitdir(@__FILE__)[1]
	dest = joinpath(splitdir(@__FILE__)[1],"temporary_files")
	if !isdir(dest)
		mkdir(dest)
	end
	mindate = DateTime(2016,11,23,14,4)
	maxdate = DateTime(2016,11,23,17)
	dt = Dates.Hour(4)
	ClimateDataIO.slt_timeshift(src,dest,mindate,maxdate,dt)
	Data = ClimateDataIO.slt_load(dest,mindate + dt,maxdate + dt,verbose=true)
	f = ["W20163281804.cfg",
		"W20163281804.csr",
		"W20163281804.csu",
		"W20163281804.csv",
		"W20163281804.flx",
		"W20163281804.log",
		"W20163281804.slt",
		"W20163281830.slt",
		"W20163281900.slt",
		"W20163281930.slt",
		"W20163282000.slt",
		"W20163282030.slt"]
	for i in f
		println(joinpath(dest,i))
		rm(joinpath(dest,i))
	end
	@test Data.Time[1] == DateTime(2016,11,23,18,4) || "SLT_TIMESHIFT: Last timestamp should be 2016-11-23T18:04:00"
	@test Data.Time[end] == DateTime(2016,11,23,21,0,0,100) || "SLT_TIMESHIFT: Last timestamp should be 2016-11-23T21:00:00.1"
	println("SLT_TIMESHIFT Test Complete")

	println("\n===  SLT_WRITE  ===")
	dest = joinpath(splitdir(@__FILE__)[1],"temporary_files")
	if !isdir(dest)
		mkdir(dest)
	end
	src = joinpath(joinpath(splitdir(@__FILE__)[1],"temporary_files"),"X")
	t0 = DateTime(2017,1,31,14)
	ch = convert(Vector{Int8},[1,1])
	bm = convert(Vector{Int8},[1,1])
	l = 40000 # Number of lines
	windxyz = 2000 .*rand(l,3) # Wind Vectors
	sos = 50. *(10. *rand(l,1) .+ 320) # Speed of Sound
	analog_signals = 5000. *rand(l,2) # Voltages
	D0 = convert(Array{Int16},floor.([windxyz sos analog_signals])) # Data to save
	ClimateDataIO.slt_write(src,t0,ch,bm,D0)
	src = joinpath(joinpath(splitdir(@__FILE__)[1],"temporary_files"),"X20170311400.slt")
	t1,D1 = ClimateDataIO.slt_read(src,2,10)
	rm(src) # Delete Temporary File
	D1[:,1:3] = D1[:,1:3].*100
	D1[:,4] = D1[:,4].*50
	D1 = convert(Array{Int16},floor.(D1))
	D2 = abs.(D1 .- D0) .> 1
	@test isempty(LinearIndices(D2)[findall(D2)]) || "SLT_WRITE: The differences between the original array and the reloaded array should be no greater than 1 (rounding errors)."
	println("SLT_WRITE Test Complete")

	println("\n===  SLT_CONFIG  ===")
	src = splitdir(@__FILE__)[1]
	cfgs = ClimateDataIO.slt_config(src)
	@test get(cfgs[1],"Time",DateTime(0)) == DateTime(2016,11,23,14,4) || "SLT_CONFIG: cfgs[:Time] should be 2016-11-23T14:04:00"
	@test length(keys(cfgs[1])) == 23 || "SLT_CONFIG: length(keys(cfgs[1])) should be 23 (ie 23 Dict keys)"

	src = joinpath(src,"W20163281404.cfg")
	cfgs = ClimateDataIO.slt_config(src)
	@test get(cfgs[1],"Time",DateTime(0)) == DateTime(2016,11,23,14,4) || "SLT_CONFIG: cfgs[:Time] should be 2016-11-23T14:04:00"
	@test length(keys(cfgs[1])) == 23 || "SLT_CONFIG: length(keys(cfgs[1])) should be 23 (ie 23 Dict keys)"
	println("SLT_CONFIG Test Complete")
end

@testset "============================ Aerodyne TESTS  ============================" begin
	# STR_LOAD: Load one Aerodyne STR file
	println("\n====  STR_LOAD Tests  ====")
	src = splitdir(@__FILE__)[1]
	src = joinpath(src,"161210_000000.str")
	cols = ["NH3"]
	Time,Data = ClimateDataIO.str_load(src,verbose=true,cols=cols)
	@test Time[1] == DateTime(2016,12,10,0,0,0,525) || "STR_LOAD: First timestamp should be 2016-12-10T00:00:00.525"
	@test Time[end] == DateTime(2016,12,10,0,59,59,922) || "STR_LOAD: Last timestamp should be 2016-12-10T00:59:59.922"
	@test size(Data) == (35160,1) || "STR_LOAD: Output size incorrect, should be (35160,2)"

	println("\nMultiple files")
	src = splitdir(@__FILE__)[1]
	Time,Data = ClimateDataIO.str_load(src,verbose=true)
	@test Time[1] == DateTime(2016,12,10,0,0,0,525) || "STR_LOAD: First timestamp should be 2016-12-10T00:00:00.525"
	@test Time[end] == DateTime(2016,12,10,3,59,59,917) || "STR_LOAD: Last timestamp should be 2016-12-10T03:59:59.917"
	@test size(Data) == (140645,2) || "STR_LOAD: Output size incorrect, should be (140645,2)"
	println("STR_LOAD Test Complete")

	# STC_LOAD: Load one Aerodyne STC file
	println("\n====  STC_LOAD Tests  ====")
	src = splitdir(@__FILE__)[1]
	src = joinpath(src,"161210_000000.stc")
	Time,Data = ClimateDataIO.stc_load(src,verbose=true)
	@test Time[1] == DateTime(2016,12,10,0,0,0,622) || "STC_LOAD: First timestamp should be 2016-12-10T00:00:02.622"
	@test Time[end] == DateTime(2016,12,10,0,59,59,622) || "STC_LOAD: Last timestamp should be 2016-12-10T00:59:59.622"
	@test size(Data) == (3600,26) || "STC_LOAD: Output size incorrect, should be (3600,26)"

	println("\nMultiple files")
	src = splitdir(@__FILE__)[1]
	Time,Data = ClimateDataIO.stc_load(src,verbose=true)
	@test Time[1] == DateTime(2016,12,10,0,0,0,622) || "STC_LOAD: First timestamp should be 2016-12-10T00:00:02.622"
	@test Time[end] == DateTime(2016,12,10,3,59,59,117) || "STC_LOAD: Last timestamp should be 2016-12-10T03:59:59.117"
	@test size(Data) == (14399,26) || "STC_LOAD: Output size incorrect, should be (14399,26)"

	println("\nStatus Values")
	statuses = ClimateDataIO.AerodyneStatus.(Data.StatusW)
	temp = fill!(Array{Bool}(undef,14399),false)
	@test temp == getfield.(statuses,:Valve1) || "AERODYNESTATUS: All values in statuses.Valve1 should be FALSE"
	@test temp == getfield.(statuses,:Valve2) || "AERODYNESTATUS: All values in statuses.Valve2 should be FALSE"
	@test temp == getfield.(statuses,:Valve3) || "AERODYNESTATUS: All values in statuses.Valve3 should be FALSE"
	@test temp != getfield.(statuses,:Valve4) || "AERODYNESTATUS: All values in statuses.Valve4 should not all be FALSE"
	@test temp == getfield.(statuses,:Valve5) || "AERODYNESTATUS: All values in statuses.Valve5 should be FALSE"
	@test temp == getfield.(statuses,:Valve6) || "AERODYNESTATUS: All values in statuses.Valve6 should be FALSE"
	@test temp == getfield.(statuses,:Valve7) || "AERODYNESTATUS: All values in statuses.Valve7 should be FALSE"
	@test temp != getfield.(statuses,:Valve8) || "AERODYNESTATUS: All values in statuses.Valve8 should not all be FALSE"
	known_names = (:AutoBG,:AutoCal,:FrequencyLock,:BinomialFilter,:AltMode,:GuessLast,:PowerNorm,:ContRefLock,:AutoSpectSave,:PressureLock,:WriteData,:RS232,:ElectronicBGSub,:Valve1,:Valve2,:Valve3,:Valve4,:Valve5,:Valve6,:Valve7,:Valve8)
	@test known_names == propertynames(statuses[1]) || "AERODYNESTATUS: statuses field names should be as follows - \n21-element Array{Symbol,1}:\n  :AutoBG\n  :AutoCal\n  :FrequencyLock\n  :BinomialFilter\n  :AltMode\n  :GuessLast\n  :PowerNorm\n  :ContRefLock\n  :AutoSpectSave\n  :PressureLock\n  :WriteData\n  :RS232\n  :ElectronicBGSub\n  :Valve1\n  :Valve2\n  :Valve3\n  :Valve4\n  :Valve5\n  :Valve6\n  :Valve7\n  :Valve8"

	src = splitdir(@__FILE__)[1]
	src = joinpath(src,"161210_000000.stc")
	Time,Data = ClimateDataIO.stc_load(src,verbose=true,cols = ["Praw","time"])
	@test Symbol.(names(Data)) == [:Praw,:time] || "STC_LOAD: Fields should be [:Praw,:time]"
	println("STC_LOAD Test Complete")
end

@testset "====================== Campbell Scientific TESTS  ======================" begin
	# CSCI_TEXTREAD: Load one file
	println("\n====  CSCI_TEXTREAD Tests  ====")
	src = splitdir(@__FILE__)[1]
	src = joinpath(src,"CampbellScientific_HC2S3_20161201-000000.dat")
	Data = ClimateDataIO.csci_textread(src,verbose=true)
	@test Data.Timestamp[1] == DateTime(2016,12,1,0,0,0) || "CSCI_TEXTREAD: Data[:Timestamp] first timestamp should be 2016-12-01T00:00:00"
	@test Data.Timestamp[end] == DateTime(2016,12,4,23,59,30) || "CSCI_TEXTREAD: Data[:Timestamp] last timestamp should be 2016-12-04T23:59:30"
	@test size(Data) == (11520,10) || "CSCI_TEXTREAD: Data output size incorrect, should be (11520,10)"

	Data,loggerStr,colsStr,unitsStr,processingStr = ClimateDataIO.csci_textread(src,verbose=true,headeroutput=true)
	h1 = "\"TOA5\",\"CR1000 - IP\",\"CR1000\",\"E2948\",\"CR1000.Std.22\",\"CPU:LoggerCode.CR1\",\"35271\",\"Rotronics_HC2S3\""
	h2 = "\"TIMESTAMP\",\"RECORD\",\"AirTemp_HC2S3_S01\",\"RelHum_HC2S3_S01\",\"AirTemp_HC2S3_S02\",\"RelHum_HC2S3_S02\",\"AirTemp_HC2S3_S03\",\"RelHum_HC2S3_S03\",\"AirTemp_HC2S3_S04\",\"RelHum_HC2S3_S04\""
	h3 = "\"TS\",\"RN\",\"Deg C\",\"%\",\"Deg C\",\"%\",\"Deg C\",\"%\",\"Deg C\",\"%\""
	h4 = "\"\",\"\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\""
	@test loggerStr == h1 || "loggerStr should be \"TOA5\",\"CR1000 - IP\",\"CR1000\",\"E2948\",\"CR1000.Std.22\",\"CPU:LoggerCode.CR1\",\"35271\",\"Rotronics_HC2S3\""
	@test colsStr == h2 || "colsStr should be \"TIMESTAMP\",\"RECORD\",\"AirTemp_HC2S3_S01\",\"RelHum_HC2S3_S01\",\"AirTemp_HC2S3_S02\",\"RelHum_HC2S3_S02\",\"AirTemp_HC2S3_S03\",\"RelHum_HC2S3_S03\",\"AirTemp_HC2S3_S04\",\"RelHum_HC2S3_S04\""
	@test unitsStr == h3 || "unitsStr should be \"TS\",\"RN\",\"Deg C\",\"%\",\"Deg C\",\"%\",\"Deg C\",\"%\",\"Deg C\",\"%\""
	@test processingStr == h4 || "processingStr should be \"\",\"\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\""
	@test Data.Timestamp[1] == DateTime(2016,12,1,0,0,0) || "CSCI_TEXTREAD: Data[:Timestamp] first timestamp should be 2016-12-01T00:00:00"
	@test Data.Timestamp[end] == DateTime(2016,12,4,23,59,30) || "CSCI_TEXTREAD: Data[:Timestamp] last timestamp should be 2016-12-04T23:59:30"
	@test size(Data) == (11520,10) || "CSCI_TEXTREAD: Data output size incorrect, should be (11520,10)"
	println("CSCI_TEXTREAD Test Complete")

	Data = ClimateDataIO.csci_textread(src,verbose=true,headeroutput=false,select=[:Timestamp,:AirTemp_HC2S3_S01,:RelHum_HC2S3_S02])
	@test size(Data) == (11520,3) || "CSCI_TEXTREAD: Data output size incorrect, should be (11520,3)"
	@test names(Data) == ["Timestamp","AirTemp_HC2S3_S01","RelHum_HC2S3_S02"] || "CSCI_TEXTREAD: Columns should be \"Timestamp\",\"AirTemp_HC2S3_S01\",\"RelHum_HC2S3_S02\""

	Data = ClimateDataIO.csci_textread(src,verbose=true,headeroutput=false,select=[:AirTemp_HC2S3_S01,:RelHum_HC2S3_S02])
	@test size(Data) == (11520,3) || "CSCI_TEXTREAD: Data output size incorrect, should be (11520,3)"
	@test names(Data) == ["Timestamp","AirTemp_HC2S3_S01","RelHum_HC2S3_S02"] || "CSCI_TEXTREAD: Columns should be \"Timestamp\",\"AirTemp_HC2S3_S01\",\"RelHum_HC2S3_S02\""

	println("\n===  CSCI_TEXTLOAD  ===")
	src = splitdir(@__FILE__)[1]
	src = joinpath(src,"CampbellScientific_HC2S3_20161201-000000.dat")
	Data, HeaderInfo, HeaderColumns, HeaderUnits, HeaderSample = ClimateDataIO.csci_textload([src],verbose=true,headerlines=4,headeroutput=true)
	h1 = "\"TOA5\",\"CR1000 - IP\",\"CR1000\",\"E2948\",\"CR1000.Std.22\",\"CPU:LoggerCode.CR1\",\"35271\",\"Rotronics_HC2S3\""
	h2 = "\"TIMESTAMP\",\"RECORD\",\"AirTemp_HC2S3_S01\",\"RelHum_HC2S3_S01\",\"AirTemp_HC2S3_S02\",\"RelHum_HC2S3_S02\",\"AirTemp_HC2S3_S03\",\"RelHum_HC2S3_S03\",\"AirTemp_HC2S3_S04\",\"RelHum_HC2S3_S04\""
	h3 = "\"TS\",\"RN\",\"Deg C\",\"%\",\"Deg C\",\"%\",\"Deg C\",\"%\",\"Deg C\",\"%\""
	h4 = "\"\",\"\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\""
	@test Data.Timestamp[1] == DateTime(2016,12,1,0,0,0) || "CSCI_TEXTREAD: Data[:Timestamp] first timestamp should be 2016-12-01T00:00:00"
	@test Data.Timestamp[end] == DateTime(2016,12,4,23,59,30) || "CSCI_TEXTREAD: Data[:Timestamp] last timestamp should be 2016-12-04T23:59:30"
	@test size(Data) == (11520,10) || "CSCI_TEXTLOAD: Data output size incorrect, should be (11520,10)"
	@test h1 == HeaderInfo || "CSCI_TEXTLOAD: HeaderInfo should be \"\"TOA5\",\"CR1000 - IP\",\"CR1000\",\"E2948\",\"CR1000.Std.22\",\"CPU:LoggerCode.CR1\",\"35271\",\"Rotronics_HC2S3\"\""
	@test h2 == HeaderColumns || "CSCI_TEXTLOAD: HeaderColumns should be \"\"TIMESTAMP\",\"RECORD\",\"AirTemp_HC2S3_S01\",\"RelHum_HC2S3_S01\",\"AirTemp_HC2S3_S02\",\"RelHum_HC2S3_S02\",\"AirTemp_HC2S3_S03\",\"RelHum_HC2S3_S03\",\"AirTemp_HC2S3_S04\",\"RelHum_HC2S3_S04\"\""
	@test h3 == HeaderUnits || "CSCI_TEXTLOAD: HeaderUnits should be \"\"TS\",\"RN\",\"Deg C\",\"%\",\"Deg C\",\"%\",\"Deg C\",\"%\",\"Deg C\",\"%\"\""
	@test h4 == HeaderSample || "CSCI_TEXTLOAD: HeaderSample should be \"\"\",\"\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\",\"Smp\"\""

	src = splitdir(@__FILE__)[1]
	rootname = "CampbellScientific"
	Data = ClimateDataIO.csci_textload(src,rootname,DateTime(2016,12),DateTime(2017,1))
	@test Data.Timestamp[1] == DateTime(2016,12,1,0,0,0) || "CSCI_TEXTREAD: Data[:Timestamp] first timestamp should be 2016-12-01T00:00:00"
	@test Data.Timestamp[end] == DateTime(2016,12,4,23,59,30) || "CSCI_TEXTREAD: Data[:Timestamp] last timestamp should be 2016-12-04T23:59:30"
	@test size(Data) == (11520,10) || "CSCI_TEXTLOAD: Data output size incorrect, should be (11520,10)"
	println("CSCI_TEXTLOAD Test Complete")

	Data = ClimateDataIO.csci_textload(src,rootname,DateTime(2016,12),DateTime(2017,1),select=[:Timestamp,:AirTemp_HC2S3_S01,:RelHum_HC2S3_S02])
	@test size(Data) == (11520,3) || "CSCI_TEXTLOAD: Data output size incorrect, should be (11520,3)"
	@test names(Data) == ["Timestamp","AirTemp_HC2S3_S01","RelHum_HC2S3_S02"] || "CSCI_TEXTLOAD: Columns should be \"Timestamp\",\"AirTemp_HC2S3_S01\",\"RelHum_HC2S3_S02\""

	println("\n===  CSCI_TIMES  ===")
	src = splitdir(@__FILE__)[1]
	src = joinpath(src,"CampbellScientific_HC2S3_20161201-000000.dat")
	srcs = [src,src]
	mintime = DateTime(2016,12,1)
	maxtime = DateTime(2016,12,4,23,59,30)
	mintimes, maxtimes = ClimateDataIO.csci_times(srcs,headerlines=4)
	@test [mintime;mintime] == mintimes || "CSCI_TIME: Both minimum time values should be 2016-12-01T00:00:00"
	@test [maxtime;maxtime] == maxtimes || "CSCI_TIME: Both maximum time values should be 2016-12-04T23:59:30"
	println("CSCI_TIMES Test Complete")
end

@testset "============================= Licor TESTS  =============================" begin
	# GHGREAD: Load one file
	println("\n====  GHGREAD Tests  ====")
	src = splitdir(@__FILE__)[1]
	src = joinpath(src,"2016-12-11T203000_AIU-1359.ghg")
	Time,Data = ClimateDataIO.ghg_read(src,verbose=true,filetype="primary")
	@test Time[1] == DateTime(2016,12,11,20,30,0) || "GHGREAD: Time[1] first timestamp should be 2016-12-11T20:30:00"
	@test Time[end] == DateTime(2016,12,11,20,59,59,950) || "GHGREAD: Time[end] last timestamp should be 2016-12-04T20:59:59.95"
	@test size(Data) == (36000,48) || "GHGREAD: Data output size incorrect, should be (36000,48)"

	Time,Data = ClimateDataIO.ghg_read(src,verbose=true,filetype="primary",select=[:Dew_Point_C,:Flow_Drive_])
	@test size(Data) == (36000,4) || "GHGREAD: Data output size incorrect, should be (36000,4)"

	# GHGLOAD: Load one file
	println("\n====  GHGLOAD Tests  ====")
	src = splitdir(@__FILE__)[1]
	Time,Data = ClimateDataIO.ghg_load(src,verbose=true,average=false,filetype="primary")
	@test Time[1] == DateTime(2016,12,11,20,0,0) || "GHGLOAD: Time[1] first timestamp should be 2016-12-11T20:00:00"
	@test Time[end] == DateTime(2016,12,11,21,59,59,950) || "GHGLOAD: Time[end] last timestamp should be 2016-12-04T21:59:59.95"
	@test size(Data) == (144000,48) || "GHGLOAD: Data output size incorrect, should be (144000,48)"

	Time,Data = ClimateDataIO.ghg_load(src,verbose=true,average=true,filetype="primary")
	@test Time[1] == DateTime(2016,12,11,20,0,0) || "GHGLOAD: Time[1] first timestamp should be 2016-12-11T20:00:00"
	@test Time[end] == DateTime(2016,12,11,21,30) || "GHGLOAD: Time[end] last timestamp should be 2016-12-04T21:30:00"
	@test size(Data) == (4,48) || "GHGLOAD: Data output size incorrect, should be (4,48)"

	Time,Data = ClimateDataIO.ghg_load(src,verbose=true,average=true,filetype="primary",select=[:Dew_Point_C,:Flow_Drive_])
	@test size(Data) == (4,4) || "GHGLOAD: Data output size incorrect, should be (4,4)"

end

#global err = true
#global Data = DataFrame()
@testset "======================= Los Gatos Research TESTS  =======================" begin
	# LGR_LOAD: Load a file
	println("\n====  LGR_LOAD Tests  ====")
	src = splitdir(@__FILE__)[1]
	Data = ClimateDataIO.lgr_load(src,verbose=true)
	@test Data.Time[1] == DateTime(2015,1,30,15,13,47,994) || "LGR_LOAD: Time[1] first timestamp should be 2015-01-30T15:13:47.994"
	@test Data.Time[end] == DateTime(2015,1,31,15,13,55,727) || "LGR_LOAD: Time[end] last timestamp should be 2015-01-31T15:13:55.727"
	@test size(Data) == (78272,23) || "LGR_LOAD: Data output size incorrect, should be (78272,23)"

	# LGR_READ Error Check
	err = true
	try
		Data = ClimateDataIO.lgr_read("blahblahblah")
		err = false
	catch

	end
	@test err == true || "LGR_READ: Should throw an error, invalid file given"
end

@testset "============================== Licor TEST  ==============================" begin
	# LICOR_SPLIT
	println("\n====  LICOR_SPLIT Tests  ====")
	src = splitdir(@__FILE__)[1]
	dest = joinpath(src,"temporary_files")
	F1 = joinpath(src,"2016-12-11T213000_AIU-1359.ghg")
	F2 = joinpath(dest,"2016-12-11T213000_AIU-1359.data")
	F3 = joinpath(dest,"2016-12-11T213000.txt")
	if Sys.isunix()
		temp = read(`unzip -d $dest $F1`,String);
	elseif Sys.iswindows()
		temp = read(`$exe7z x $F1 -y -o$dest`,String);
	end

	# Prepare File
	println("\nGenerating test file...")
	epoch = DateTime(1970) + Dates.Hour(1)

	# Copy the old file
	cp(F2,F3)

	# Add Data
	fid = open(F2,"r")
	fid2 = open(F3,"a")
	#l = ""
	#t = []
	#ft = []
	#s = []
	for i=1:1:8
    	# Skip Header
    	readline(fid,keep=true)
	end
	p = position(fid) # File position
	for j=1:1:4 # 20 for another eight hours
		seek(fid,p) # Reset file position
		global l = ""
		t = []
		ft = []
    	while !eof(fid)
			l = readline(fid,keep=true)
        	ft = findall(x -> x == '\t',l)
        	t = epoch + Dates.Second(Meta.parse(l[ft[1]+1:ft[2]])) + Dates.Millisecond(round(parse(Float64,l[ft[2]+1:ft[3]-1])/10^6))
        	t = t + Dates.Second(1800*j)
			l = l[1:5] * string(parse(Int,l[ft[1]+1:ft[2]]) + 1800*j) * l[ft[2]:ft[6]] * Dates.format(t,"yyyy-mm-dd\tHH:MM:SS:sss") * l[ft[8]:end]
			write(fid2,l)
    	end
	end
	close(fid)
	close(fid2)
	rm(joinpath(dest,"2016-12-11T213000_AIU-1359.data"))
	rm(joinpath(dest,"2016-12-11T213000_AIU-1359.metadata"))

	println("\nSplitting test file...")
	ClimateDataIO.licor_split(dest,dest)
	rm(joinpath(dest,"2016-12-11T213000.txt"))

	println("\nLoading Split Files...")
	Time,Data = ClimateDataIO.ghg_load(dest,DateTime(2016,12,11,21,30),average=false,verbose=false,filetype="primary")
	rm(joinpath(dest,"2016-12-11T213000_AIU-1359.ghg"))
	#rm(joinpath(dest,"2016-12-12T000000_AIU-1359.ghg"))
	#rm(joinpath(dest,"2016-12-12T040000_AIU-1359.ghg"))
	@test Time[1] == DateTime(2016,12,11,21,30) || "LICOR_SPLIT: Time[1] first timestamp should be 2016-12-11T21:30:00"
	@test Time[end] == DateTime(2016,12,11,23,59,59,950) || "LICOR_SPLIT: Time[end] last timestamp should be 2016-12-11T23:59:59.95"
	@test size(Data) == (180000,48) || "LGR_SPLIT: Data output size incorrect, should be (180000,48)"

	src = splitdir(@__FILE__)[1]
	if isdir(joinpath(src,"temporary_files"))
		rm(joinpath(src,"temporary_files"))
	end
end

println("\n\nAll Tests Complete Successfully")

ext: jl

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module UDChains

import ..DebugMsg
DebugMsg.init()

using CompilerTools

import Base.show

"""
Contains the UDchains for one basic block.
"""
type UDInfo
    live_in  :: Dict{Symbol,Set}
    live_out :: Dict{Symbol,Set}

    function UDInfo()
        new(Dict{Symbol,Set}(), Dict{Symbol,Set}())
    end
end

"""
Get the set of definition blocks reaching this block for a given symbol "s".
"""
function getOrCreate(live :: Dict{Symbol, Set}, s :: Symbol)
    if !haskey(live, s)
        live[s] = Set()
    end
    return live[s]
end

"""
Get the UDInfo for a specified basic block "bb" or create one if it doesn't already exist.
"""
function getOrCreate(udchains :: Dict{CompilerTools.LivenessAnalysis.BasicBlock,UDInfo} , bb :: CompilerTools.LivenessAnalysis.BasicBlock)
    if !haskey(udchains, bb)
        udchains[bb] = UDInfo()
    end
    return udchains[bb]
end

"""
Print the set part of a live in or live out dictiononary in a nice way if the debug level is set high enough.
"""
function printSet(level, s)
    for j in s
        if typeof(j) == CompilerTools.LivenessAnalysis.BasicBlock
            dprint(level, " ", j.label)
        else
            dprint(level, " ", j)
        end
    end
end

"""
Print a live in or live out dictionary in a nice way if the debug level is set high enough.
"""
function printLabels(level, dict)
    for i in dict
       dprint(level, "\tSymbol: ", i[1], " From:")
       printSet(level, i[2])
       dprintln(level,"")
    end
end

"""
Print UDChains in a nice way if the debug level is set high enough.
"""
function printUDInfo(level, ud)
    for i in ud
        dprintln(level, i[1].label, " Live In:")
        printLabels(level, i[2].live_in)
        dprintln(level, i[1].label, " Live Out:")
        printLabels(level, i[2].live_out)
    end
end

"""
Get the Use-Definition chains at a basic block level given LivenessAnalysis.BlockLiveness as input in "bl".
"""
function getUDChains(bl :: CompilerTools.LivenessAnalysis.BlockLiveness)
    udchains = Dict{CompilerTools.LivenessAnalysis.BasicBlock,UDInfo}()

    @dprintln(3,"getUDChains: bl = ", bl)

    body_order = CompilerTools.LivenessAnalysis.getBbBodyOrder(bl)
    changed = true
    # Iterate until nothing changes.
    while changed
        @dprintln(3,"getUDChains: main loop")
        changed = false
        # For each basic block from beginning to end.
        for i = 1:length(body_order)
            bb = bl.basic_blocks[body_order[i]]
            # Get the UDChain info for this basic block.
            udinfo = getOrCreate(udchains, bb)
            @dprintln(3,"getUDChains: bb = ", bb, " udinfo = ", udinfo)

            for li in bb.live_in
                # Get the current set of blocks from which this symbol could be defined and reach here.
                li_set = getOrCreate(udinfo.live_in, li)
                @dprint(3,"getUDChains: li = ", li, " li_set = ")
                printSet(3,li_set)
                @dprintln(3,"")

                if isempty(bb.preds)
                    @dprintln(3,"getUDChains: no preds")
                    # Must be the starting block.
                    # Use "nothing" to represent the parameter set.
                    if !in(nothing, li_set)
                        @dprintln(3,"getUDChains: added nothing to li_set")
                        push!(li_set, nothing)
                        changed = true
                    end
                else
                    # Non-starting block.
                    for pred in bb.preds
                        pred_udinfo = getOrCreate(udchains, pred)
                        pred_lo_set = getOrCreate(pred_udinfo.live_out, li)
                        #@dprint(3,"getUDChains: pred = ", pred.label, " pred_udinfo = ", pred_udinfo, " pred_lo_set = ", pred_lo_set)
                        @dprint(3,"getUDChains: pred = ", pred.label, " pred_lo_set = ", pred_lo_set)
                        printSet(3,li_set)
                        @dprintln(3,"")
                        for pred_lo in pred_lo_set
                            if !in(pred_lo, li_set)
                                push!(li_set, pred_lo)
                                changed = true
                                @dprintln(3,"getUDChains: added ", pred_lo, " to li_set = ", li_set)
                            end
                        end
                    end
                end
            end

            # For each symbol live_out of this block.
            for lo in bb.live_out
                # Get the current set of blocks from which this symbol could be defined and reach here.
                lo_set = getOrCreate(udinfo.live_out, lo)
                @dprint(3,"getUDChains: lo = ", lo, " lo_set = ")
                printSet(3,lo_set)
                @dprintln(3,"")

                # If this live out was defined in this block...
                if in(lo, bb.def)
                    @dprintln(3,"getUDChains: lo def in block")
                    # ... then add this block to udchain for this symbol.
                    if !in(bb, lo_set)
                        @dprintln(3,"getUDChains: adding bb to lo_set")
                        push!(lo_set, bb)
                        changed = true
                    end
                else
                    # Not defined in this block so must be live_in.
                    li_set = getOrCreate(udinfo.live_in, lo)
                    @dprintln(3,"getUDChains: not def in block so using li_set = ", li_set)
                    for li_bb in li_set
                        if !in(li_bb, lo_set)
                            @dprintln(3,"getUDChains: adding ", typeof(li_bb)==CompilerTools.LivenessAnalysis.BasicBlock ? li_bb.label : li_bb, " to lo_set")
                            push!(lo_set, li_bb)
                            changed = true
                        end
                    end
                end
            end
        end
    end

    return udchains
end

end


text:
function Base.filter(f::Function, A::TimeAxisArray)

    n = size(A,1)
    predresults = Array{Bool}(n)
    tstamps = timestamp(A)

    for i in 1:n
        ts = tstamps[i]
        predvals[i] = f(ts, A[Axis{:Timestamp}(ts)])
    end #for

    return A[Axis{:Timestamp}(predvals)]
end #filter

filtertimestamps(f::Function, A::TimeAxisArray) = filter((ts,v)->f(ts), A)
filtervalues(f::Function, A::TimeAxisArray) = filter((ts,v)->f(v), A)

"""
    split(A::TimeAxisArray, f::Function)

Returns an array containing sequential fragments of `A`, split according to clusters of values in the mapping of `f` over the timestamps of `A`. `split(f, A)` is also defined so as to support do-notation.
"""
function Base.split(A::TimeAxisArray, f::Function)

    discriminants = map(f, timestamps(A))
    n = length(discriminants)

    clusterstart = 1
    splits = UnitRange{Int}[]
    for i in 2:n
        if discriminants[i] !== discriminants[i-1]
            push!(splits, clusterstart:i-1)
            clusterstart = i
        end #if
    end #for
    push!(splits, clusterstart:n)

    return map(s -> A[Axis{:Timestamp}(s)], splits)

end #split

# For use with do-notation
Base.split(f::Function, A::TimeAxisArray) = Base.split(A, f)

applytimewise(f::Function, A::AbstractArray) = if ndims(A) == 1
    return [f(A)]
else
    higherdimsizes = size(A)[2:end]
    result = Array{typeof(f(A[1]))}(1, higherdimsizes...)
    for inds in product(map(n->1:n, higherdimsizes)...)
        result[1, inds...] = f(A[:, inds...])
    end #for
    return result
end #applytimewise

"""
    collapse(A::TimeAxisArray, tsreducer::Function, reducer::Function=tsreducer)

Collapses the timestamps of `A` to a single observation in the time dimension using `tsreducer`. Data is collapsed in the time dimension using `reducer`, which defaults to `tsreducer`.
"""
collapse(A::TimeAxisArray, tsreducer::Function, reducer::Function=tsreducer; lift::Bool=true) = lift ?
    TimeAxisArray(applytimewise(reducer, A.data), [tsreducer(timestamps(A))], A.axes[2:end]...) :
    TimeAxisArray(reducer(A.data), [tsreducer(timestamps(A))], A.axes[2:end]...)

"""
    downsample(A::TimeAxisArray, splitter::Function, tsreducer::Function, reducer::Function=tsreducer)

Combines `split`, `collapse`, and `vcat` to partition `A` according to sequential values in the mapping of `splitter` over the timestamps of `A`, then collapses each of the split TimeAxisArrays according to `tsreducer` (for timestamps) and `reducer` (for data), before recombining the collapsed values.
"""
downsample(A::TimeAxisArray, splitter::Function, tsreducer::Function, reducer::Function=tsreducer; lift::Bool=true) =
    cat(1, map(a -> collapse(a, tsreducer, reducer, lift=lift), split(A, splitter))...)

# # TODO: Allow for time-based interval selection (rather than fixed integer range window)
"""
    moving(A::TimeAxisArray, reducer::Function, n::Int)

Applies a time-wise reduction specified by `reducer` to a moving window of `n` observations, storing the result at the last timestamp in the window.
"""
function moving{T}(A::TimeAxisArray{T}, reducer::Function, window::Int; lift::Bool=true)
    if lift
        data = applytimewise(reducer, A[Axis{:Timestamp}(1:window)].data)
        for i in 2:(size(A,1)+1-window)
            data = vcat(data, applytimewise(reducer, A[Axis{:Timestamp}(i:i+window-1)].data))
        end #for
    else
        data = reducer(A[Axis{:Timestamp}(1:window)].data)
        for i in 2:(size(A,1)+1-window)
            data = vcat(data, reducer(A[Axis{:Timestamp}(i:i+window-1)].data))
        end #for
    end #if
    return TimeAxisArray(data, timestamps(A)[window:end], A.axes[2:end]...)
end #moving

"""
    lag(A::TimeAxisArray, k::Int=1)

Shifts all observations in `A` later in time by `k` timestamps.
"""
function lag(A::TimeAxisArray, k::Int=1)
    n = size(A,1)
    return TimeAxisArray(A[Axis{:Timestamp}(1:n-k)].data, timestamps(A)[1+k:n], A.axes[2:end]...)
end #lag

"""
    lead(A::TimeAxisArray, k::Int=1)

Shifts all observations in `A` earlier in time by `k` timestamps.
"""
function lead(A::TimeAxisArray, k::Int=1)
    n = size(A,1)
    return TimeAxisArray(A[Axis{:Timestamp}(1+k:n)].data, timestamps(A)[1:n-k], A.axes[2:end]...)
end #lead

"""
    diff(A::TimeAxisArray, k::Int=1)

Perform `k`th order differencing across time observations in `A`.
"""
Base.diff(A::TimeAxisArray, k::Int=1) = if k == 1
    n = size(A,1)
    return TimeAxisArray(A[Axis{:Timestamp}(2:n)].data .- A[Axis{:Timestamp}(1:n-1)].data, timestamps(A)[2:n], A.axes[2:end]...)
else
    return diff(diff(A, k-1))
end #if

"""
    percentchange(A::TimeAxisArray; logdiff::Bool=false)

Computes the percent change between observations in time in `A`. If `logdiff` is true, returns the difference of log-transformed values.
"""
function percentchange(A::TimeAxisArray; logdiff::Bool=false)
    n = size(A,1)
    return logdiff ? diff(map(log,A)) : TimeAxisArray(A[Axis{:Timestamp}(2:n)] ./  A[Axis{:Timestamp}(1:n-1)] .- 1, timestamps(A)[2:end], A.axes[2:end]...)
end #percentchange

"""
    dropif(selector::Function, predicate::Function, A::TimeAxisArray)

Drops observations at timestamps where `selector` (e.g. `any`, `all`) data values statisfy `predicate`.
"""
function dropif(selector::Function, predicate::Function, A::TimeAxisArray)
    keepers = Int[]
    for i in 1:size(A,1)
        !selector(predicate, A[Axis{:Timestamp}(i)]) && push!(keepers, i)
    end #for
    return A[Axis{:Timestamp}(keepers)]
end #dropif

"""
    dropnan(selector::Function, A::TimeAxisArray)

Drops observations at timestamps where `selector` (e.g. `any`, `all`) data values are NaN. Equivalent to `dropif(selector, isnan, A)`
"""
dropnan(A::TimeAxisArray, selector::Function=all) = dropif(selector, isnan, A)

max_stars_repo_name: GordStephen/TimeAxisArrays

max_issues_repo_name: GordStephen/TimeAxisArrays

max_forks_repo_name: GordStephen/TimeAxisArrays


text:
module Xkcd221

export Xkcd221RNG

using Random: Random, AbstractRNG, Sampler, SamplerType, LessThan

import Random: rand, seed!

mutable struct Xkcd221RNG <: AbstractRNG
    seed::UInt64
    Xkcd221RNG(seed::Integer = 4) = new(seed)
end

function Base.copy!(dst::Xkcd221RNG, src::Xkcd221RNG)
    dst.seed = src.seed
    dst
end

Base.copy(src::Xkcd221RNG) = Xkcd221RNG(src.seed)

Base.:(==)(x::Xkcd221RNG, y::Xkcd221RNG) = x.seed == y.seed

Base.hash(rng::Xkcd221RNG, h::UInt) = hash(rng.seed, 0x2bad40c9493e0718 % UInt ⊻ h)

## Sampling

rand(rng::Xkcd221RNG, ::SamplerType{UInt64}) = rng.seed

for T = [Bool, Base.BitInteger_types...]
    T === UInt64 && continue
    @eval rand(rng::Xkcd221RNG, ::SamplerType{$T}) = rand(rng, UInt64) % $T
end

Random.rng_native_52(::Xkcd221RNG) = UInt64

function rand(rng::Xkcd221RNG, sp::LessThan)
    x = rand(rng, sp.s)
    sp.sup < rng.seed ? x % sp.sup  : x
end

## Range Sampling

struct RangeSampler{T} <: Sampler{T}
    s::T
    m::UInt64
end

function RangeSampler(r::AbstractUnitRange{T}) where {T}
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    m = (last(r) - first(r)) % unsigned(T) % UInt64
    RangeSampler(first(r), m)
end

function rand(rng::Xkcd221RNG, sp::RangeSampler{T}) where {T}
    x = !iszero(sp.m) ? ((sp.m>rng.seed ? rng.seed : (rng.seed%sp.m))) : 0
    (sp.s % UInt64 + x) % T
end

for T in Base.BitInteger_types
    @eval Sampler(::Type{Xkcd221RNG}, rg::AbstractUnitRange{$T}, ::Random.Repetition) = RangeSampler(rg)
end

end


text:
using Plots

using simulator.fresnelltools
using simulator.matrixcore
using simulator.materials.spesifics

save = false

n_1 = SiO2_core_Sellmeier
n_2 = Au_Johnson
n_3 = Pd052_unloaded
n_4 = Air
n_5 = Air

d_1 = 30e-9
d_2 = 30e-9
d_3 = 5e-9
d_4 = 10e-9
d_5 = 10e-9

λ = 1200e-9
θd = 45

step = 1e-10
expression_total(Up, Un) = abs(Up + Un)^2
expression_forward(Up, _) = abs(Up)^2
expression_backward(_, Un) = abs(Un)^2
expression_borders(u) = abs(u)^2

description = "silver silica 4"
θds = 0:1e-2:90
λs = 400e-9:1e-9:1500e-9

fig_path = "../rapport/figures/interrogator_example/"
fig_token = replace(description, " " => "_")

if save
    gr()
else
    plotly()
end
θ1 = θd * π/180
θs = θds .* (π/180)

function system_slabs(θ1, λ, select)
    bulk_1 = FresnellSlab(n_1(λ), 2*π/λ, d_1, θ1)
    i_1_2_te, i_1_2_tm, θ2 = FresnellBoundrary(n_1(λ), n_2(λ), θ1)
    bulk_2 = FresnellSlab(n_2(λ), 2*π/λ, d_2, θ2)
    i_2_3_te, i_2_3_tm, θ3 = FresnellBoundrary(n_2(λ), n_3(λ), θ2)
    bulk_3 = FresnellSlab(n_3(λ), 2*π/λ, d_3, θ3)
    i_3_4_te, i_3_4_tm, θ4 = FresnellBoundrary(n_3(λ), n_4(λ), θ3)
    bulk_4 = FresnellSlab(n_4(λ), 2*π/λ, d_4, θ4)
    i_4_5_te, i_4_5_tm, θ5 = FresnellBoundrary(n_4(λ), n_5(λ), θ4)
    bulk_5 = FresnellSlab(n_5(λ), 2*π/λ, d_5, θ5)

    system_te = CascadeScattering([bulk_1, i_1_2_te, bulk_2, i_2_3_te, bulk_3, i_3_4_te, bulk_4, i_4_5_te, bulk_5])
    system_tm = CascadeScattering([bulk_1, i_1_2_tm, bulk_2, i_2_3_tm, bulk_3, i_3_4_tm, bulk_4, i_4_5_tm, bulk_5])

    te_r = abs(system_te[2, 1])^2
    te_t = abs(system_te[1, 1])^2
    tm_r = abs(system_tm[2, 1])^2
    tm_t = abs(system_tm[1, 1])^2

    [te_r, tm_r, te_t, tm_t][select]
end

function inside_slabs(expression)
    bulk_1(d) = FresnellSlab(n_1(λ), 2*π/λ, d, θ1)

    i_1_2_te, i_1_2_tm, θ2 = FresnellBoundrary(n_1(λ), n_2(λ), θ1)
    i_1_2_te_f(d) = i_1_2_te
    i_1_2_tm_f(d) = i_1_2_tm

    bulk_2(d) = FresnellSlab(n_2(λ), 2*π/λ, d, θ2)

    i_2_3_te, i_2_3_tm, θ3 = FresnellBoundrary(n_2(λ), n_3(λ), θ2)
    i_2_3_te_f(d) = i_2_3_te
    i_2_3_tm_f(d) = i_2_3_tm

    bulk_3(d) = FresnellSlab(n_3(λ), 2*π/λ, d, θ3)

    i_3_4_te, i_3_4_tm, θ4 = FresnellBoundrary(n_3(λ), n_4(λ), θ3)
    i_3_4_te_f(d) = i_3_4_te
    i_3_4_tm_f(d) = i_3_4_tm

    bulk_4(d) = FresnellSlab(n_4(λ), 2*π/λ, d, θ4)

    i_4_5_te, i_4_5_tm, θ5 = FresnellBoundrary(n_4(λ), n_5(λ), θ4)
    i_4_5_te_f(d) = i_4_5_te
    i_4_5_tm_f(d) = i_4_5_tm

    bulk_5(d) = FresnellSlab(n_5(λ), 2*π/λ, d, θ5)

    u_te, d_te = Interrogator(
        [bulk_1, i_1_2_te_f, bulk_2, i_2_3_te_f, bulk_3, i_3_4_te_f, bulk_4, i_4_5_te_f, bulk_5],
        [d_1, 0, d_2, 0, d_3, 0, d_4, 0, d_5],
        step,
        expression
    )
    u_tm, d_tm = Interrogator(
        [bulk_1, i_1_2_tm_f, bulk_2, i_2_3_tm_f, bulk_3, i_3_4_tm_f, bulk_4, i_4_5_tm_f, bulk_5],
        [d_1, 0, d_2, 0, d_3, 0, d_4, 0, d_5],
        step,
        expression
    )

    u_bp_te, u_bn_te, d_b_te = BorderInterrogator(
        [bulk_1(d_1), i_1_2_te, bulk_2(d_2), i_2_3_te, bulk_3(d_3), i_3_4_te, bulk_4(d_4), i_4_5_te, bulk_5(d_5)],
        [d_1, 0, d_2, 0, d_3, 0, d_4, 0, d_5],
        expression_borders
    )
    u_bp_tm, u_bn_tm, d_b_tm = BorderInterrogator(
        [bulk_1(d_1), i_1_2_tm, bulk_2(d_2), i_2_3_tm, bulk_3(d_3), i_3_4_tm, bulk_4(d_4), i_4_5_tm, bulk_5(d_5)],
        [d_1, 0, d_2, 0, d_3, 0, d_4, 0, d_5],
        expression_borders
    )

    u_te, d_te, u_tm, d_tm, u_bp_te, u_bn_te, d_b_te, u_bp_tm, u_bn_tm, d_b_tm
end

# println("\nTiming for inside")
# total_u_te, total_d_te, total_u_tm, total_d_tm, u_bp_te, u_bn_te, d_b_te, u_bp_tm, u_bn_tm, d_b_tm = @time inside_slabs(expression_total)
# forward_u_te, forward_d_te, forward_u_tm, forward_d_tm, _, _, _, _, _, _ = @time inside_slabs(expression_forward)
# backward_u_te, backward_d_te, backward_u_tm, backward_d_tm, _, _, _, _, _, _ = @time inside_slabs(expression_backward)

# fig_inside = plot(title="Field inside " * description * "\nλ = $λ, θ = $θd", legend=:topleft, xlabel="distance in nm", ylabel="power", ticks=:native)
# plot!(fig_inside, total_d_te .* 1e9, total_u_te, label="total te")
# plot!(fig_inside, total_d_tm .* 1e9, total_u_tm, label="total tm")
# plot!(fig_inside, forward_d_te .* 1e9, forward_u_te, label="forward te")
# plot!(fig_inside, d_b_te .* 1e9, u_bp_te, seriestype =:scatter, label="forward te borders")
# plot!(fig_inside, forward_d_tm .* 1e9, forward_u_tm, label="forward tm")
# plot!(fig_inside, d_b_tm .* 1e9, u_bp_tm, seriestype =:scatter, label="forward tm borders")
# plot!(fig_inside, backward_d_te .* 1e9, backward_u_te, label="backward te")
# plot!(fig_inside, d_b_te .* 1e9, u_bn_te, seriestype =:scatter, label="backward te borders")
# plot!(fig_inside, backward_d_tm .* 1e9, backward_u_tm, label="backward tm")
# plot!(fig_inside, d_b_tm .* 1e9, u_bn_tm, seriestype =:scatter, label="backward tm borders")
# display(fig_inside)
# if save savefig(fig_inside, fig_path * "inside_" * fig_token * ".pdf") end

println("\nTiming for by angle")
ter_angles = @time system_slabs.(θs, λ, 1)
tmr_angles = @time system_slabs.(θs, λ, 2)
tet_angles = @time system_slabs.(θs, λ, 3)
tmt_angles = @time system_slabs.(θs, λ, 4)
fig_angles = plot(title="By angle " * description * "\nλ = $λ", legend=:left, xlabel="angle in degrees", ylabel="power", ticks=:native)
plot!(fig_angles, θs .* (180/π), ter_angles, label="te reflection")
plot!(fig_angles, θs .* (180/π), tmr_angles, label="tm reflection")
# plot!(fig_angles, θs .* (180/π), tet_angles, label="te transmission")
# plot!(fig_angles, θs .* (180/π), tmt_angles, label="tm transmission")
display(fig_angles)
if save savefig(fig_angles, fig_path * "angles_" * fig_token * ".pdf") end

println("\nTiming for by wavelength")
ter_wavelength = @time system_slabs.(θ1, λs, 1)
tmr_wavelength = @time system_slabs.(θ1, λs, 2)
tet_wavelength = @time system_slabs.(θ1, λs, 3)
tmt_wavelength = @time system_slabs.(θ1, λs, 4)
fig_wavelengths = plot(title="By wavelength " * description * "\nθ = $θd", legend=:right, xlabel="wavelength in nanometer", ylabel="power", ticks=:native)
plot!(fig_wavelengths, λs .* 1e9, ter_wavelength, label="te reflection")
plot!(fig_wavelengths, λs .* 1e9, tmr_wavelength, label="tm reflection")
# plot!(fig_wavelengths, λs .* 1e9, tet_wavelength, label="te transmission")
# plot!(fig_wavelengths, λs .* 1e9, tmt_wavelength, label="tm transmission")
display(fig_wavelengths)
if save savefig(fig_wavelengths, fig_path * "wavelengths_" * fig_token * ".pdf") end


text:
@defstruct RandomMaskLayer Layer (
  name :: AbstractString = "random-mask",
  (ratio :: AbstractFloat = 0.5, 0 < ratio < 1),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) > 0)
)
@characterize_layer(RandomMaskLayer,
  can_do_bp => true,
  is_inplace => true,
)

type RandomMaskLayerState <: LayerState
  layer      :: RandomMaskLayer

  dropouts   :: Vector{DropoutLayerState}
end

function setup(backend::Backend, layer::RandomMaskLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  dropouts = Array(DropoutLayerState, length(inputs))
  for i = 1:length(inputs)
    dropout_layer = DropoutLayer(name="$(layer.name)-dropout-$i", auto_scale=false, ratio=layer.ratio,
        bottoms=Symbol[symbol("$(layer.bottoms[i])-$i")])
    dropouts[i] = setup(backend, dropout_layer, Blob[inputs[i]], Blob[diffs[i]])
  end
  return RandomMaskLayerState(layer, dropouts)
end

function shutdown(backend::Backend, state::RandomMaskLayerState)
  for dropout_state in state.dropouts
    shutdown(backend, dropout_state)
  end
end

function forward(backend::Backend, state::RandomMaskLayerState, inputs::Vector{Blob})
  for i = 1:length(inputs)
    forward(backend, state.dropouts[i], Blob[inputs[i]])
  end
end

function backward(backend::Backend, state::RandomMaskLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  for i = 1:length(inputs)
    backward(backend, state.dropouts[i], Blob[inputs[i]], Blob[diffs[i]])
  end
end


text:
# Julia wrapper for header: /usr/include/scip/pub_bandit_epsgreedy.h
# Automatically generated using Clang.jl wrap_c

function SCIPcreateBanditEpsgreedy(scip, epsgreedy, priorities, eps, preferrecent, decayfactor, avglim, nactions, initseed)
    ccall((:SCIPcreateBanditEpsgreedy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_BANDIT}}, Ptr{Cdouble}, Cdouble, UInt32, Cdouble, Cint, Cint, UInt32), scip, epsgreedy, priorities, eps, preferrecent, decayfactor, avglim, nactions, initseed)
end

function SCIPgetWeightsEpsgreedy(epsgreedy)
    ccall((:SCIPgetWeightsEpsgreedy, libscip), Ptr{Cdouble}, (Ptr{SCIP_BANDIT},), epsgreedy)
end

function SCIPsetEpsilonEpsgreedy(epsgreedy, eps)
    ccall((:SCIPsetEpsilonEpsgreedy, libscip), Cvoid, (Ptr{SCIP_BANDIT}, Cdouble), epsgreedy, eps)
end


text:
cd(@__DIR__)

using Pkg; Pkg.activate(".."); Pkg.instantiate()

using CSV, DataFrames, DelimitedFiles, ProgressMeter, Pipe

#---

dsts = readdlm("../data/levenshteinLanguageDistances.csv", ',')

languages = string.(sort(unique(vcat(dsts[1,2:end], dsts[2:end,1]))))

nonProto = filter(x -> !occursin("PROTO", x), languages)

nPairs = length(nonProto)*(length(nonProto)-1) ÷ 2

dstArray = zeros(length(languages), length(languages))
for j in 1:length(languages)
    for i in 2:length(languages)
        if j>=i
            dstArray[i-1,j] = dsts[i,j]
        end
    end
end

prInd = Tuple.(findall(dstArray .<= 0.7))
filter!(x -> x[1] < x[2], prInd)

probablyRelated = @pipe [
          [languages[x[1]], languages[x[2]]]
          for
          x in prInd if languages[x[1]] ∈ nonProto && languages[x[2]] ∈ nonProto
      ] |>
      hcat(_...) |>
      permutedims

writedlm("../data/probablyRelated.csv", probablyRelated, ',')

#---
d = CSV.File("../data/asjp19wide.csv") |> DataFrame!

prI = indexin(probablyRelated, d.longname)

wordPairs = Tuple{String, String}[]

@showprogress for i in 1:size(prI,1)
    j,k = prI[i,:]
    wl1 = Vector(d[j,2:end])
    wl2 = Vector(d[k,2:end])
    for (w1,w2) in zip(wl1, wl2)
        if !ismissing(w1) && !ismissing(w2)
            for x in split(w1,"-")
                for y in split(w2,"-")
                    push!(wordPairs, (x,y))
                end
            end
        end
    end
end

@time wordPairs = permutedims(hcat(collect.(wordPairs)...))

writedlm("../data/potentialCognates.csv", wordPairs, ',')


text:
function DfpField(arg0::jint)
    return DfpField((jint,), arg0)
end

function clear_ieee_flags(obj::DfpField)
    return jcall(obj, "clearIEEEFlags", void, ())
end

function compute_exp(::Type{DfpField}, arg0::Dfp, arg1::Dfp)
    return jcall(DfpField, "computeExp", Dfp, (Dfp, Dfp), arg0, arg1)
end

function compute_ln(::Type{DfpField}, arg0::Dfp, arg1::Dfp, arg2::Dfp)
    return jcall(DfpField, "computeLn", Dfp, (Dfp, Dfp, Dfp), arg0, arg1, arg2)
end

function equals(obj::DfpField, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_deg_to_rad(obj::DfpField)
    return jcall(obj, "getDegToRad", Dfp, ())
end

function get_e(obj::DfpField)
    return jcall(obj, "getE", Dfp, ())
end

function get_e_split(obj::DfpField)
    return jcall(obj, "getESplit", Vector{Dfp}, ())
end

function get_extended_field(obj::DfpField, arg0::jint, arg1::jboolean)
    return jcall(obj, "getExtendedField", DfpField, (jint, jboolean), arg0, arg1)
end

function get_ieee_flags(obj::DfpField)
    return jcall(obj, "getIEEEFlags", jint, ())
end

function get_ln10(obj::DfpField)
    return jcall(obj, "getLn10", Dfp, ())
end

function get_ln2(obj::DfpField)
    return jcall(obj, "getLn2", Dfp, ())
end

function get_ln2_split(obj::DfpField)
    return jcall(obj, "getLn2Split", Vector{Dfp}, ())
end

function get_ln5(obj::DfpField)
    return jcall(obj, "getLn5", Dfp, ())
end

function get_ln5_split(obj::DfpField)
    return jcall(obj, "getLn5Split", Vector{Dfp}, ())
end

function get_one(obj::DfpField)
    return jcall(obj, "getOne", Dfp, ())
end

function get_pi(obj::DfpField)
    return jcall(obj, "getPi", Dfp, ())
end

function get_pi_split(obj::DfpField)
    return jcall(obj, "getPiSplit", Vector{Dfp}, ())
end

function get_rad_to_deg(obj::DfpField)
    return jcall(obj, "getRadToDeg", Dfp, ())
end

function get_radix_digits(obj::DfpField)
    return jcall(obj, "getRadixDigits", jint, ())
end

function get_rounding_mode(obj::DfpField)
    return jcall(obj, "getRoundingMode", DfpField_RoundingMode, ())
end

function get_runtime_class(obj::DfpField)
    return jcall(obj, "getRuntimeClass", Class, ())
end

function get_sqr2(obj::DfpField)
    return jcall(obj, "getSqr2", Dfp, ())
end

function get_sqr2_reciprocal(obj::DfpField)
    return jcall(obj, "getSqr2Reciprocal", Dfp, ())
end

function get_sqr2_split(obj::DfpField)
    return jcall(obj, "getSqr2Split", Vector{Dfp}, ())
end

function get_sqr3(obj::DfpField)
    return jcall(obj, "getSqr3", Dfp, ())
end

function get_sqr3_reciprocal(obj::DfpField)
    return jcall(obj, "getSqr3Reciprocal", Dfp, ())
end

function get_two(obj::DfpField)
    return jcall(obj, "getTwo", Dfp, ())
end

function get_zero(obj::DfpField)
    return jcall(obj, "getZero", Dfp, ())
end

function hash_code(obj::DfpField)
    return jcall(obj, "hashCode", jint, ())
end

function new_dfp(obj::DfpField)
    return jcall(obj, "newDfp", Dfp, ())
end

function new_dfp(obj::DfpField, arg0::Dfp)
    return jcall(obj, "newDfp", Dfp, (Dfp,), arg0)
end

function new_dfp(obj::DfpField, arg0::JString)
    return jcall(obj, "newDfp", Dfp, (JString,), arg0)
end

function new_dfp(obj::DfpField, arg0::jbyte)
    return jcall(obj, "newDfp", Dfp, (jbyte,), arg0)
end

function new_dfp(obj::DfpField, arg0::jbyte, arg1::jbyte)
    return jcall(obj, "newDfp", Dfp, (jbyte, jbyte), arg0, arg1)
end

function new_dfp(obj::DfpField, arg0::jdouble)
    return jcall(obj, "newDfp", Dfp, (jdouble,), arg0)
end

function new_dfp(obj::DfpField, arg0::jint)
    return jcall(obj, "newDfp", Dfp, (jint,), arg0)
end

function new_dfp(obj::DfpField, arg0::jlong)
    return jcall(obj, "newDfp", Dfp, (jlong,), arg0)
end

function set_ieee_flags(obj::DfpField, arg0::jint)
    return jcall(obj, "setIEEEFlags", void, (jint,), arg0)
end

function set_ieee_flags_bits(obj::DfpField, arg0::jint)
    return jcall(obj, "setIEEEFlagsBits", void, (jint,), arg0)
end

function set_rounding_mode(obj::DfpField, arg0::DfpField_RoundingMode)
    return jcall(obj, "setRoundingMode", void, (DfpField_RoundingMode,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end


text:
""" Visualization for Delta Sets and their duals.

This module provides wrapper functions for plotting embedded delta sets
using the Makie functionality. This includes the wireframe, mesh, and scatter
plotting functions.

It is important to note that the end-user will need to import one of the Makie
backends (CairoMakie, GLMakie, or WGLMakie) in order for plotting to work
correctly.
"""

module MeshGraphics

using Catlab.CategoricalAlgebra
using ..SimplicialSets
using ..MeshInterop
#import ...AbstractPlotting: wireframe, wireframe!, mesh, mesh!, scatter, scatter!
#
#export wireframe, wireframe!, mesh, mesh!, scatter, scatter!

import ...AbstractPlotting
import ...AbstractPlotting: convert_arguments

""" This extends the "Mesh" plotting recipe for embedded deltasets by converting
an embedded deltaset into arguments that can be passed into AbstractPlotting.mesh
"""

function convert_arguments(P::Union{Type{<:AbstractPlotting.Wireframe},
                                    Type{<:AbstractPlotting.Mesh},
                                    Type{<:AbstractPlotting.Scatter}},
                           dset::AbstractACSet)
  convert_arguments(P, Mesh(dset))
end

""" This extends the "LineSegments" plotting recipe for embedded deltasets by converting
an embedded deltaset into arguments that can be passed into AbstractPlotting.linesegments
"""
function convert_arguments(P::Type{<:AbstractPlotting.LineSegments}, dset::EmbeddedDeltaSet2D)
  edge_positions = zeros(ne(dset)*2,3)
  for e in edges(dset)
    edge_positions[2*e-1,:] = point(dset, dset[e,:src])
    edge_positions[2*e,:] = point(dset, dset[e,:tgt])
  end
  convert_arguments(P, edge_positions)
end

plottype(::EmbeddedDeltaSet2D) = Mesh

end


text:
struct LoggerPlay
    io::IO
end

struct LoggerPlayMinder <: Actors.AbsMinder end

hear(s::Scene{LoggerPlayMinder}, ::Left!) = say(s, stage(s), Leave!())

function hear(s::Scene{LoggerPlay}, ::Genesis!)
    log = invite!(s, Actors.Logger(my(s).io))
    m1 = invite!(s, LoggerPlayMinder())
    m2 = enter!(s, Actors.PassiveMinder(log), m1)

    delegate(s, m2) do s
        @say_info s "Noise"
        error("Drama")
    end
end

@testset "Logger Test" begin
    io = Base.BufferStream()

    play!(LoggerPlay(io))

    close(io)
    s = String(take!(io.buffer))

    @test occursin("Noise", s)
    @test occursin("Drama", s)
end


text:
# Module-scope global variables

# TODO: improve this interface:
# - should be more dict-like: get and setindex(::name), haskey(::name)
# - globals(::Type)?

export
    CuGlobal, get, set

"""
    CuGlobal{T}(mod::CuModule, name::String)

Acquires a typed global variable handle from a named global in a module.
"""
struct CuGlobal{T}
    buf::Mem.DeviceBuffer

    function CuGlobal{T}(mod::CuModule, name::String) where T
        ptr_ref = Ref{CuPtr{Cvoid}}()
        nbytes_ref = Ref{Csize_t}()
        cuModuleGetGlobal_v2(ptr_ref, nbytes_ref, mod, name)
        if nbytes_ref[] != sizeof(T)
            throw(ArgumentError("size of global '$name' does not match type parameter type $T"))
        end
        buf = Mem.DeviceBuffer(ptr_ref[], nbytes_ref[])

        return new{T}(buf)
    end
end

Base.cconvert(::Type{CuPtr{Cvoid}}, var::CuGlobal) = var.buf

Base.:(==)(a::CuGlobal, b::CuGlobal) = a.handle == b.handle
Base.hash(var::CuGlobal, h::UInt) = hash(var.ptr, h)

"""
    eltype(var::CuGlobal)

Return the element type of a global variable object.
"""
Base.eltype(::Type{CuGlobal{T}}) where {T} = T

"""
    Base.getindex(var::CuGlobal)

Return the current value of a global variable.
"""
function Base.getindex(var::CuGlobal{T}) where T
    val_ref = Ref{T}()
    cuMemcpyDtoH_v2(val_ref, var, var.buf.bytesize)
    return val_ref[]
end
# TODO: import Base: get?

"""
    Base.setindex(var::CuGlobal{T}, val::T)

Set the value of a global variable to `val`
"""
function Base.setindex!(var::CuGlobal{T}, val::T) where T
    val_ref = Ref{T}(val)
    cuMemcpyHtoD_v2(var, val_ref, var.buf.bytesize)
end


text:
function update_sig2!(i::Int, s::State, c::Constants, d::Data)
  ss = 0.0

  cardinality = 0

  for j in 1:d.J
    for n in 1:d.N[i]
      k = s.lam[i][n]
      if k > 0
        z = s.Z[j, k]
        l = s.gam[i][n, j]
        ss += (s.y_imputed[i][n, j] - mus(z, l, s, c, d)) ^ 2
        cardinality += 1
      end
    end
  end

  newShape = shape(c.sig2_prior) + cardinality / (2 * c.temper)
  newScale = scale(c.sig2_prior) + ss / (2 * c.temper)

  if c.sig2_range == [0, Inf]
    s.sig2[i] = rand(InverseGamma(newShape, newScale))
  else
    lower, upper = c.sig2_range
    s.sig2[i] = rand_uptrunc_ig(InverseGamma(newShape, newScale), upper)
  end
end

function update_sig2!(s::State, c::Constants, d::Data)
  for i in 1:d.I
    update_sig2!(i, s, c, d)
  end
end


text:
module R_sep_Model
    using JuMP
    # using LinearAlgebra

    srcDir = dirname(dirname(@__FILE__))
    include(srcDir*"\\Constraints\\R_Constraints.jl")

    using .R_constraints
    const rcon = R_constraints

    export R_Modelξₜˢᵉᵖ
    """The model"""
    function Modelξₜˢᵉᵖ(ρ,d,t;con_list="S1G")
        model = Model()
        @variable(model, Lx[rcon.make_mon_expo_keys(sum(d),t[1])] ) # Create variables
        function set_con(c)
            for k in keys(c)
                 size(c[k]) == (1, 1) ?
                 @constraint(model, c[k] .>= 0) :
                 @constraint(model, c[k] in PSDCone())
            end
        end
## PSD Moment matrix
        PSD_MM_con = rcon.make_PSD_con(d,t,Lx)
        set_con(PSD_MM_con) # L(ηᵦ) ⪰ 0 ∀ ηᵦ ∈  ⨁ᵦ ηᵦ = L([x,y]ₜ[x,y]ₜᵀ)
## L(xxᵀ ⊗ yyᵀ) = ρ
        @constraint(model, rcon.make_ord4_con(d,Lx) .==  ρ)
## Zero propagation: ρ = aaᵀ⊗bbᵀ , ρᵢⱼᵢⱼ = 0 ⟹ aᵢbⱼ = 0 ⟹ L(xᵢyⱼ) = 0
        # zero_moms = rcon.zeroprop(ρ,d,t,Lx)
        # @constraint(model,zero_moms .== zeros(size(zero_moms)))
## Localizing g constraint: L ≥ 0 on M₂ₜ(S)
        if     occursin("S1",con_list)
            println("----------S1 constraints active----------")
            set_con(rcon.make_loc_cons_S1(ρ,d,t,Lx))
        elseif occursin("S2",con_list)
            println("----------S2 constraints active----------")
            loc_con, g₂  =  rcon.make_loc_cons_S2(ρ,d,t,Lx)
            for k in keys(g₂)
                @constraint(model, g₂[k] .==  zeros(size(g₂[k])))
            end
            set_con(loc_con)
        elseif occursin("S3",con_list)
            println("----------S3 constraints active----------")
            loc_con, g₂  =  rcon.make_loc_cons_S3(ρ,d,t,Lx)
            for k in keys(g₂)
                @constraint(model, g₂[k] .==  zeros(size(g₂[k])))
            end
            set_con(loc_con)
        end

## G Constraints
        if  occursin("G",con_list)
            println("----------G-constraints are active----------")
            set_con(rcon.make_G_con(ρ,d,t,Lx)) # L((ρ-xxᵀ⊗yyᵀ)⊗([x,y]ₜ₋₂[x,y]ᵀₜ₋₂))
        end

        @objective(model, Min, Lx[zeros(Int64,sum(d))]) #  Set objective
        return model
    end
end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
mutable struct BDSVI{T}
	dobs::AbstractArray{T,3}
	dcal::AbstractArray{T,3}
	ϕ::Array{T,2}
	x::Array{T,2}
	is0::Int
	bdr::BD{T}
	bds::BD{T}
	opt::Solver
end

function BDSVI(dobs, is0, opt::Solver)
	nt,nr,ns=size(dobs)
	ntslag=nt
	T=eltype(dobs)

	bdr=FBD.BD(nt, nt, nr, 2*ntslag-1, dobs=randn(T, nt,nr), slags=[ntslag-1,ntslag-1],
	    	fft_threads=true, fftwflag=FFTW.MEASURE);
	bds=FBD.BD(2*ntslag-1, nt, ns, nt, dobs=randn(T, nt,ns), glags=[ntslag-1,ntslag-1],
	   slags=[nt-1, 0],
	    fft_threads=true, fftwflag=FFTW.MEASURE);

	ϕ=zeros(T,2*ntslag-1,ns)
	x=zeros(T,nt,nr)

	if(typeof(opt) == UseGradNMF)
		dobs=abs.(dobs)
	end

	pa=BDSVI(dobs,zero(dobs),ϕ,x,is0,bdr,bds,opt)

	initialize!(pa)
	return pa
end

# initialize
function initialize!(pa::BDSVI)
	nt,nr,ns=size(pa.dobs)
	ntslag=nt

	if((ns-pa.is0)*(1-pa.is0) ≤ 0)
		for ir in 1:nr, it in 1:nt
			pa.x[it,ir]=pa.dobs[it,ir,pa.is0]
		end
	else
		randn!(pa.x)
	end

	randn!(pa.ϕ)
	if(typeof(pa.opt) == UseGradNMF)
		pa.ϕ=1.0 .+ abs.(randn(size(pa.ϕ)))
	end

	if((ns-pa.is0)*(1-pa.is0) ≤ 0)
		for it in 1:2*ntslag-1
			pa.ϕ[it,pa.is0]=0
		end
		pa.ϕ[ntslag,pa.is0]=1.0
	end

	# compute other phi
	if(typeof(pa.opt) == UseGradNMF)
		for i in 1:10  # just perform 10 iterations, change later?!
			fitϕ!(pa, pa.opt)
		end
	else
		fitϕ!(pa, UseIterativeSolvers())
	end
	return pa
end

# for each receiver, update x
function fitx!(pa::BDSVI, opt::Solver)

	nt,nr,ns=size(pa.dobs)
	J0=0.0
	for ir in 1:nr
		d=pa.bds.optm.obs.d
		dobs=view(pa.dobs,:,ir,:)
		for i in eachindex(d)
			d[i]=dobs[i]
		end
		for i in eachindex(pa.ϕ)
			pa.bds.optm.cal.g[i]=pa.ϕ[i]
		end

		# initial value
		xs=view(pa.x,:,ir)
		for i in eachindex(xs)
			pa.bds.optm.cal.s[i]=xs[i]
		end

		J1=FBD.update!(pa.bds, pa.bds.sx.x, FBD.S(), opt)

		J0 += J1

		for is in 1:ns, it in 1:nt
			pa.dcal[it,ir,is]=pa.bds.optm.cal.d[it,is]
		end

		copyto!(xs,pa.bds.sx.x)

	end
	return J0
end

function fitϕ!(pa::BDSVI, opt::Solver)

	nt,nr,ns=size(pa.dobs)
	J0=0.0
	for iss in 1:ns
		d=pa.bdr.optm.obs.d
		dobs=view(pa.dobs,:,:,iss)
		for i in eachindex(d)
			d[i]=dobs[i]
		end
		for i in eachindex(pa.x)
			pa.bdr.optm.cal.g[i]=pa.x[i]
		end

		# initial value of ϕ
		ϕs=view(pa.ϕ,:,iss)
		for i in eachindex(pa.bdr.optm.cal.s)
			pa.bdr.optm.cal.s[i]=ϕs[i]
		end

		if(((pa.is0-ns)*(1-pa.is0) ≤ 0) & (iss == pa.is0))
			ntm=div(length(ϕs),2)
			for it in vcat(1:ntm, ntm+2:length(ϕs))
				ϕs[it]=0.0
			end
			J1=0.
		else
			J1 = FBD.update!(pa.bdr, pa.bdr.sx.x, FBD.S(), opt)
			copyto!(ϕs,pa.bdr.sx.x)

		end
		J0 += J1

		for ir in 1:nr, it in 1:nt
			pa.dcal[it,ir,iss]=pa.bdr.optm.cal.d[it,ir]
		end

	end
	return J0
end

function fit!(pa::BDSVI, io=stdout;
		     max_roundtrips=100,
		     max_reroundtrips=1,
		     ParamAM_func=nothing,
		     roundtrip_tol=1e-6,
		     optim_tols=[1e-6, 1e-6], verbose=true,
		     log_file=false,
		     )

	global to, optG, optS
	reset_timer!(to)

	if(ParamAM_func===nothing)
		ParamAM_func=xx->ParamAM(xx, optim_tols=optim_tols,name="Blind Decon",
				    roundtrip_tol=roundtrip_tol, max_roundtrips=max_roundtrips,
				    max_reroundtrips=max_reroundtrips,
				    min_roundtrips=10,
				    verbose=verbose,
				    reinit_func=xxx->initialize!(pa),
				    #after_roundtrip_func=x->(write(io,string(TimerOutputs.flatten(to),"\n"))),
				    #after_roundtrip_func=x->(err!(pa, nothing)),
				    )
	end

	# create alternating minimization parameters
	f1=x-> fitϕ!(pa, pa.opt)
	f2=x-> fitx!(pa, pa.opt)
	paam=ParamAM_func([f1, f2])

	if(io===nothing)
		logfilename=joinpath(pwd(),string("XBDSVI",Dates.now(),".log"))
		io=open(logfilename, "a+")
	end

	go(paam, io)  # run alternative minimization

	# save log file
	#if(log_file)
	#	AMlogfile=joinpath(pwd(),string("AM",Dates.now(),".log"))
	#	CSV.write(AMlogfile, paam.log)
	#end

	write(io,string(to))
	write(io, "\n")

	if(io === nothing)
		close(io)
	end
end


text:
export Codecov

"""
Coverage.Codecov Module

This module provides functionality to push coverage information to the CodeCov.io
web service. It exports the `submit` and `submit_local` methods.
"""
module Codecov
    using HTTP
    using Coverage
    using CoverageTools
    using JSON
    using LibGit2

    export submit, submit_token, submit_local, submit_generic

    #=
    JSON structure for Codecov.io
    https://codecov.io/api#post-report

    {
      "coverage": {
        "path/to/file.py": [null, 1, 0, null, true, 0, 0, 1, 1],
        "path/to/other.py": [null, 0, 1, 1, "1/3", null]
      },
      "messages": {
        "path/to/other.py": {
          "1": "custom message for line 1"
        }
      }
    }
    =#

    # Turn vector of filename : coverage pairs into a dictionary
    function to_json(fcs::Vector{FileCoverage})
        cov = Dict()
        for fc in fcs
            cov[fc.filename] = vcat(nothing, fc.coverage)
        end
        return Dict("coverage" => cov)
    end

    """
    kwargs provides default values to insert into args_array, only if they are
    not already specified in args_array.
    """
    function set_defaults(args_array; kwargs...)
        defined_names = keys(args_array)
        is_args_array = Pair{Symbol, Any}[]
        is_args_array = append!(is_args_array, args_array)
        for kwarg in kwargs
            if !(kwarg[1] in defined_names)
                push!(is_args_array, kwarg)
            end
        end
        return is_args_array
    end

    """
        submit(fcs::Vector{FileCoverage})

    Takes a vector of file coverage results (produced by `process_folder`),
    and submits them to Codecov.io. Assumes that this code is being run
    on TravisCI or AppVeyor. If running locally, use `submit_local`.
    """
    function submit(fcs::Vector{FileCoverage}; kwargs...)
        submit_generic(fcs; add_ci_to_kwargs(;kwargs...)...)
    end

    function add_ci_to_kwargs(;kwargs...)
        if lowercase(get(ENV, "APPVEYOR", "false")) == "true"
            appveyor_pr = get(ENV, "APPVEYOR_PULL_REQUEST_NUMBER", "")
            appveyor_job = join(
                [
                    ENV["APPVEYOR_ACCOUNT_NAME"],
                    ENV["APPVEYOR_PROJECT_SLUG"],
                    ENV["APPVEYOR_BUILD_VERSION"],
                ],
                "%2F",
            )
            kwargs = set_defaults(kwargs,
                service      = "appveyor",
                branch       = ENV["APPVEYOR_REPO_BRANCH"],
                commit       = ENV["APPVEYOR_REPO_COMMIT"],
                pull_request = appveyor_pr,
                job          = appveyor_job,
                slug         = ENV["APPVEYOR_REPO_NAME"],
                build        = ENV["APPVEYOR_JOB_ID"],
            )
        elseif lowercase(get(ENV, "TRAVIS", "false")) == "true"
            kwargs = set_defaults(kwargs,
                service      = "travis-org",
                branch       = ENV["TRAVIS_BRANCH"],
                commit       = ENV["TRAVIS_COMMIT"],
                pull_request = ENV["TRAVIS_PULL_REQUEST"],
                job          = ENV["TRAVIS_JOB_ID"],
                slug         = ENV["TRAVIS_REPO_SLUG"],
                build        = ENV["TRAVIS_JOB_NUMBER"],
            )
        elseif lowercase(get(ENV, "CIRCLECI", "false")) == "true"
            circle_slug = join(
                [
                    ENV["CIRCLE_PROJECT_USERNAME"],
                    ENV["CIRCLE_PROJECT_REPONAME"],
                ],
                "%2F",
            )
            kwargs = set_defaults(kwargs,
                service      = "circleci",
                branch       = ENV["CIRCLE_BRANCH"],
                commit       = ENV["CIRCLE_SHA1"],
                pull_request = get(ENV, "CIRCLE_PR_NUMBER", "false"),  # like Travis
                build_url    = ENV["CIRCLE_BUILD_URL"],
                slug         = circle_slug,
                build        = ENV["CIRCLE_BUILD_NUM"],
            )
        elseif lowercase(get(ENV, "JENKINS", "false")) == "true"
            kwargs = set_defaults(kwargs,
                service      = "jenkins",
                branch       = ENV["GIT_BRANCH"],
                commit       = ENV["GIT_COMMIT"],
                job          = ENV["JOB_NAME"],
                build        = ENV["BUILD_ID"],
                build_url    = ENV["BUILD_URL"],
                jenkins_url  = ENV["JENKINS_URL"],
            )
        elseif haskey(ENV, "BUILD_BUILDURI") # Azure Pipelines
            ref = get(ENV, "SYSTEM_PULLREQUEST_TARGETBRANCH", ENV["BUILD_SOURCEBRANCHNAME"])
            branch = startswith(ref, "refs/heads/") ? ref[12:end] : ref
            kwargs = set_defaults(kwargs,
                service      = "azure_pipelines",
                branch       = branch,
                commit       = ENV["BUILD_SOURCEVERSION"],
                pull_request = get(ENV, "SYSTEM_PULLREQUEST_PULLREQUESTNUMBER", ""),
                job          = ENV["BUILD_DEFINITIONNAME"],
                slug         = ENV["BUILD_REPOSITORY_NAME"],
                build        = ENV["BUILD_BUILDID"],
            )
        elseif haskey(ENV, "GITHUB_ACTION") # GitHub Actions
            kwargs = set_defaults(kwargs,
                service      = "custom",
                commit       = ENV["GITHUB_SHA"],
                slug         = ENV["GITHUB_REPOSITORY"],
            )
        else
            error("No compatible CI platform detected")
        end

        return kwargs
    end

    """
        submit_local(fcs::Vector{FileCoverage}, dir::AbstractString=pwd())

    Take a `Vector` of file coverage results (produced by `process_folder`),
    and submit them to Codecov.io. Assumes the submission is being made from
    a local git installation, rooted at `dir`. A repository token should be specified by a
    `token` keyword argument or the `CODECOV_TOKEN` environment variable.
    """
    function submit_local(fcs::Vector{FileCoverage}, dir::AbstractString=pwd(); kwargs...)
        submit_generic(fcs; add_local_to_kwargs(dir; kwargs...)...)
    end

    function add_local_to_kwargs(dir; kwargs...)
        LibGit2.with(LibGit2.GitRepoExt(dir)) do repo
            LibGit2.with(LibGit2.head(repo)) do headref
                branch_name = LibGit2.shortname(headref) # this function returns a String
                commit_oid  = LibGit2.GitHash(LibGit2.peel(headref))
                kwargs = set_defaults(kwargs,
                    commit = string(commit_oid),
                    branch = branch_name
                    )
            end
        end

        if haskey(ENV, "REPO_TOKEN")
            @warn "the environment variable REPO_TOKEN is deprecated, use CODECOV_TOKEN instead"
            kwargs = set_defaults(kwargs, token = ENV["REPO_TOKEN"])
        end

        return kwargs
    end

    @deprecate submit_token submit_local

    """
        submit_generic(fcs::Vector{FileCoverage})

    Takes a vector of file coverage results (produced by `process_folder`),
    and submits them to a Codecov.io instance. Keyword arguments are converted
    into a generic Codecov.io API uri.  It is essential that the keywords and
    values match the Codecov upload/v2 API specification.
    The `codecov_url` keyword argument or the CODECOV_URL environment variable
    can be used to specify the base path of the uri.
    The `dry_run` keyword can be used to prevent the http request from
    being generated.
    """
    function submit_generic(fcs::Vector{FileCoverage}; kwargs...)
        @assert length(kwargs) > 0
        dry_run = get(kwargs, :dry_run, false)
        if haskey(kwargs, :verbose)
            Base.depwarn("The verbose keyword argument is deprecated, set the environment variable " *
                         "JULIA_DEBUG=Coverage for verbose output", :submit_generic)
            verbose = kwargs[:verbose]
        else
            verbose = false
        end
        uri_str = construct_uri_string(;kwargs...)

        verbose && @info "Submitting data to Codecov..."
        verbose && @debug "Codecov.io API URL:\n" * mask_token(uri_str)

        if !dry_run
            heads   = Dict("Content-Type" => "application/json")
            data    = to_json(fcs)
            req     = HTTP.post(uri_str; body = JSON.json(data), headers = heads)
            @debug "Result of submission:" * String(req)
        end
    end

    function construct_uri_string(;kwargs...)
        if haskey(ENV, "CODECOV_URL")
            kwargs = set_defaults(kwargs, codecov_url = ENV["CODECOV_URL"])
        end

        if haskey(ENV, "CODECOV_TOKEN")
            kwargs = set_defaults(kwargs, token = ENV["CODECOV_TOKEN"])
        end

        codecov_url = "https://codecov.io"
        for (k,v) in kwargs
            if k == :codecov_url
                codecov_url = v
            end
        end
        @assert codecov_url[end] != "/" "the codecov_url should not end with a /, given url $(codecov_url)"

        uri_str = "$(codecov_url)/upload/v2?"
        for (k,v) in kwargs
            # add all except a few special key/value pairs to the URL
            # (:verbose is there for backwards compatibility with versions
            # of this code that treated it in a special way)
            if k != :codecov_url && k != :dry_run && k != :verbose
                uri_str = "$(uri_str)&$(k)=$(v)"
            end
        end

        return uri_str
    end

    function mask_token(uri_string)
        return replace(uri_string, r"token=[^&]*" => "token=<HIDDEN>")
    end

end  # module Codecov


text:
#2020 Klaus Stricker - Tests for ReadWriteDlm2
#License is MIT: http://julialang.org/license

# runtests.jl

using DelimitedFiles
using ReadWriteDlm2
using Test
using Random
using Dates
import Tables

import ReadWriteDlm2.dfregex

include("rwd2tests_1.jl")
include("rwd2tests_2.jl")
include("rwd2tests_3.jl")
include("rwd2tests_4.jl")
include("rwd2tests_5.jl")
include("rwd2tests_6.jl")


text:
include("helper.jl")

const ExamplesDir = joinpath(dirname(@__FILE__()), "..", "examples")
examples = sort(filter(fn -> isfile(joinpath(ExamplesDir, fn)) && occursin(r"\.jl$", fn),
    readdir(ExamplesDir)))

function run_examples_as_test(examplefile::String, dirpath::String)
    testcode = quote
        printstyled("Running example: ", $examplefile, "\n", color=:blue)
        @testset $examplefile begin
            err = nothing
            try
                include(joinpath($dirpath, $examplefile))
            catch _err
                err = _err
            end
            @test err == nothing
        end
    end
    eval(testcode)
end

@testset "Examples as tests" begin
    map(efn -> run_examples_as_test(efn, ExamplesDir), examples)
end


text:
using Benchmark

#------------------------------------------------------------
# On my laptop
#
# julia:
# | Row | Function             | Average   | Relative | Replications |
# |-----|----------------------|-----------|----------|--------------|
# | 1   | "pool_single_proc"   | 0.0320213 | 1.0      | 50           |
# | 2   | "pool_multiple_proc" | 0.067476  | 2.10722  | 50           |
#
# julia -p 2:
# | Row | Function             | Average   | Relative | Replications |
# |-----|----------------------|-----------|----------|--------------|
# | 1   | "pool_single_proc"   | 0.0333833 | 1.0      | 50           |
# | 2   | "pool_multiple_proc" | 4.48791   | 134.436  | 50           |
#
#------------------------------------------------------------
# On a server with
# julia -p 8:
# | Row | Function             | Average   | Relative | Replications |
# |-----|----------------------|-----------|----------|--------------|
# | 1   | "pool_single_proc"   | 0.0387865 | 1.0      | 50           |
# | 2   | "pool_multiple_proc" | 5.88758   | 151.795  | 50           |

require("parallel-pool-module")

input = rand(28, 28, 20, 128)
kernel_dim = (5, 5)
stride_dim = (2, 2)

output = pool_single_proc(input, :max, kernel_dim, stride_dim)
output2 = pool_multiple_proc(input, :max, kernel_dim, stride_dim)
@assert all(abs.(output-output2) .< 1e-10)

pool_single_proc() = pool_single_proc(input, :max, kernel_dim, stride_dim)
pool_multiple_proc() = pool_multiple_proc(input, :max, kernel_dim, stride_dim)

df = compare([pool_single_proc, pool_multiple_proc], 50)
println("$df")


text:
@testset "Partition" begin
    p = Parquet.orderedPartition(5, 2)
    expect = [[4, 1], [1, 4], [2, 3], [3, 2]]
    @test Set(p) == Set(expect)

    p = Parquet.orderedPartition(3, 2, 0)
    expect = [[3, 0], [0, 3], [1, 2], [2, 1]]
    @test Set(p) == Set(expect)
end

@testset "FindFirstIdx" begin
    function testLoopIdx(partition, firstidx, expected)
        firstLoopIdx, total = Parquet.findFirstLoopIdx(partition, firstidx)
        @test firstLoopIdx == expected
        totalExp = sum(partition) + firstidx - 1
        @test total == totalExp
    end

    testLoopIdx([1, 1, 2, 1], 1, [1, 2, 3, 5])
    testLoopIdx([1, 1, 2, 1], 0, [0, 1, 2, 4])
    testLoopIdx([1, 0, 2, 0], 1, [1, 2, 2, 4])
    testLoopIdx([1,], 1, [1,])

    function testTauIdx(partition, isG, firstidx, tauNum, expected)
        firstIdx, total = Parquet.findFirstTauIdx(partition, isG, firstidx, tauNum)
        @test firstIdx == expected
    end
    tauNum = 1
    # isG = [false, true, false, true]
    isG = [Ver4Diag, GreenDiag, Ver4Diag, GreenDiag]
    testTauIdx([1, 1, 2, 1], isG, 1, tauNum, [1, 3, 4, 7])
    testTauIdx([1, 1, 2, 1], isG, 0, tauNum, [0, 2, 3, 6])
    testTauIdx([1, 0, 2, 0], isG, 1, tauNum, [1, 3, 3, 6])

end

@testset "Filter" begin

    # for G irreducible diagrams, only 0-loop G is allowed
    @test Parquet.isValidG([Girreducible,], 0) == true
    @test Parquet.isValidG([Girreducible,], 1) == false
    @test Parquet.isValidG([Girreducible,], 2) == false

    # for Fock irreducible diagrams, only 0-loop or 2, 3, 4...-loop G is allowed
    @test Parquet.isValidG([NoFock,], 0) == true
    @test Parquet.isValidG([NoFock,], 1) == false
    @test Parquet.isValidG([NoFock,], 2) == true

    # for G irreducible diagrams, no sigma subdiagram is allowed
    @test Parquet.isValidSigma([Girreducible,], 0, true) == false
    @test Parquet.isValidSigma([Girreducible,], 1, true) == false
    @test Parquet.isValidSigma([Girreducible,], 2, true) == false

    @test Parquet.isValidSigma([Girreducible,], 0, false) == false
    @test Parquet.isValidSigma([Girreducible,], 1, false) == true
    @test Parquet.isValidSigma([Girreducible,], 2, false) == true

    # for Fock irreducible diagrams, no Fock sigma subdiagram is allowed
    @test Parquet.isValidSigma([NoFock,], 0, true) == false
    @test Parquet.isValidSigma([NoFock,], 1, true) == false
    @test Parquet.isValidSigma([NoFock,], 2, true) == true

    @test Parquet.isValidSigma([NoFock,], 0, false) == false
    @test Parquet.isValidSigma([NoFock,], 1, false) == true
    @test Parquet.isValidSigma([NoFock,], 2, false) == true
end

function evalG(K, τin, τout)
    # println(τBasis, ", ", varT)
    kF, β = 1.0, 1.0
    ϵ = dot(K, K) / 2 - kF^2
    τ = τout - τin
    if τ ≈ 0.0
        return Spectral.kernelFermiT(-1e-8, ϵ, β)
    else
        return Spectral.kernelFermiT(τ, ϵ, β)
    end
end
evalV(K) = 8π / (dot(K, K) + 1)

evalGfixK(K, τin, τout) = evalG(zero(K), τin, τout)
evalVfixK(K) = 1.0

evalFakeG(K, τin, τout) = 1.0
evalFakeV(K) = 1.0

################## api for expression tree ##############################
evalPropagator(id::BareGreenId, K, extT, varT) = evalG(K, varT[extT[1]], varT[extT[2]])
evalPropagator(id::BareInteractionId, K, extT, varT) = evalV(K)
evalPropagatorfixK(id::BareGreenId, K, extT, varT) = evalGfixK(K, varT[extT[1]], varT[extT[2]])
evalPropagatorfixK(id::BareInteractionId, K, extT, varT) = evalVfixK(K)
evalFakePropagator(id::PropagatorId, K, extT, varT) = 1.0

@testset "ParquetNew Ver4" begin
    Benchmark = Parquet.Benchmark

    function getfunction(type)
        if type == :physical
            return evalG, evalV, evalPropagator
        elseif type == :fixK
            return evalGfixK, evalVfixK, evalPropagatorfixK
        elseif type == :fake
            return evalFakeG, evalFakeV, evalFakePropagator
        else
            error("not implemented")
        end
    end

    function testVertex4(loopNum, chan, type::Symbol; filter = [], timing = false, toeval = true)
        println("$(Int.(chan)) Channel Test")
        Kdim, spin = 3, 2
        interactionTauNum = 1
        isFermi = true

        K0 = zeros(loopNum + 2)
        KinL, KoutL, KinR, KoutR = deepcopy(K0), deepcopy(K0), deepcopy(K0), deepcopy(K0)
        KinL[1] = KoutL[1] = 1
        KinR[2] = KoutR[2] = 1
        legK = [KinL, KoutL, KinR, KoutR]

        blocks = ParquetBlocks(phi = [PHEr, PPr], ppi = [PHr, PHEr])

        para = GenericPara(
            diagType = Ver4Diag,
            loopDim = Kdim,
            isFermi = isFermi,
            hasTau = true,
            innerLoopNum = loopNum,
            totalLoopNum = length(KinL),
            totalTauNum = (loopNum + 1) * interactionTauNum,
            spin = spin,
            weightType = Float64,
            firstLoopIdx = 3,
            firstTauIdx = 1,
            filter = union(filter, [Girreducible,]), #ver4 evaluation only support one-particle-irreducible diagram
            transferLoop = KinL - KoutL,
            interaction = [Interaction(ChargeCharge, Instant),],
            extra = blocks
        )

        varK = rand(Kdim, para.totalLoopNum)
        varT = [rand() for i in 1:para.totalTauNum]

        #################### DiagTree ####################################
        diags = Parquet.vertex4(para, legK, chan)
        diags = mergeby(diags, :response)
        # DiagTreeNew.plot_tree(diags[1])
        # DiagTreeNew.plot_tree(diags[2])

        ################### ExprTree ###################################
        tree = ExprTree.build(diags.diagram)
        # println("root", root)

        ################### original Parquet builder ###################################
        ver4 = Benchmark.Ver4{Benchmark.Weight}(para, Int.(chan), Int.(blocks.phi), Int.(blocks.ppi))

        if toeval

            evalG, evalV, evalPropagator = getfunction(type)

            # w1 = DiagTree.evalNaive(diag, varK, varT, evalPropagator)
            evalDiagTree!(diags, varK, varT, evalPropagator)
            w1 = [diags.diagram[1].weight, diags.diagram[2].weight]
            if timing
                printstyled("naive DiagTree evaluator cost:", color = :green)
                @time evalDiagTree!(diags, varK, varT, evalPropagator)
            end

            ExprTree.evalNaive!(tree, varK, varT, evalPropagator)
            w1e = [tree[1], tree[2]]
            if timing
                printstyled("naive ExprTree cost:", color = :green)
                @time ExprTree.evalNaive!(tree, varK, varT, evalPropagator)
            end

            optdiags = DiagTree.optimize(diags.diagram)
            opttree = ExprTree.build(optdiags)
            ExprTree.evalNaive!(opttree, varK, varT, evalPropagator)
            w1eopt = [opttree[1], opttree[2]]

            if timing
                printstyled("naive optimized ExprTree cost:", color = :green)
                @time ExprTree.evalNaive!(opttree, varK, varT, evalPropagator)
            end

            ##################### lower level subroutines  #######################################

            KinL, KoutL, KinR, KoutR = varK[:, 1], varK[:, 1], varK[:, 2], varK[:, 2]
            legK = [KinL, KoutL, KinR, KoutR]
            # Benchmark.eval(para, ver4, varK, varT, [KinL, KoutL, KinR, KoutR], evalG, evalV, true)
            Benchmark.eval(para, ver4, varK, varT, legK, evalG, evalV, true)

            if timing
                printstyled("parquet evaluator cost:", color = :green)
                # @btime sin(p, ver4, var) setup = (x = rand())
                # @time Benchmark.eval(para, ver4, varK, varT, [KinL, KoutL, KinR, KoutR], evalG, evalV, true)
                @time Benchmark.eval(para, ver4, varK, varT, legK, evalG, evalV, true)
                # @btime Benchmark.eval(p, v4, vK, vT, lK, eG, eV, flag) setup = (p = para, v4 = ver4, vK = varK, vT = varT, l = legK, eG = evalG, eV = evalV, flag = true)
            end

            w2 = ver4.weight[1]

            # println(w1, " vs ", w1e, " vs ", w2)

            @assert w1 ≈ w1e
            @assert w1 ≈ w1eopt

            # The upup channel of charge-charge vertex4 == Direct + exchange
            @test w1[1] ≈ w2[1] + w2[2]
            # The updown channel of charge-charge vertex4 == Direct
            @test w1[2] ≈ w2[1]

        end

        return para, diags, ver4
    end

    function testEval(type)
        for l = 1:3
            testVertex4(l, [PHr,], type)
            testVertex4(l, [PHEr,], type)
            testVertex4(l, [PPr,], type)
            testVertex4(l, [PHr, PHEr, PPr], type; timing = true)
        end
    end

    # testEval(:fake)
    # testEval(:fixK)
    testEval(:physical)

    #test only proper diagrams are generated if the switch is turned on
    # para, diag, ver4 = testVertex4(3, [Parquet.T, Parquet.U, Parquet.S], :physical; filter = [Builder.Proper], eval = false)
    # for i in 1:length(diag.basisPool)
    #     @test (diag.basisPool.basis[:, i] ≈ para.transferLoop) == false
    # end
end

@testset "Parquet Sigma" begin
    function getSigma(loopNum; Kdim = 3, spin = 2, interactionTauNum = 1, filter = [], isFermi = true, subdiagram = false)
        println("LoopNum =$loopNum Sigma Test")

        para = GenericPara(
            diagType = SigmaDiag,
            loopDim = Kdim,
            hasTau = true,
            innerLoopNum = loopNum,
            totalLoopNum = loopNum + 1,
            totalTauNum = loopNum * interactionTauNum,
            isFermi = isFermi,
            spin = spin,
            weightType = Float64,
            firstLoopIdx = 2,
            firstTauIdx = 1,
            filter = filter,
            interaction = [Interaction(ChargeCharge, Instant),],
            extra = ParquetBlocks(phi = [PHEr, PPr], ppi = [PHr, PHEr])
        )

        extK = zeros(para.totalLoopNum)
        extK[1] = 1.0

        varK = rand(Kdim, para.totalLoopNum)
        varT = [rand() for i in 1:para.totalTauNum]

        #################### DiagTree ####################################
        diag = Parquet.sigma(para, extK, subdiagram)
        diag = mergeby(diag)
        # print_tree(diag.diagram[1])

        return para, diag.diagram[1], varK, varT
    end

    function testDiagramNumber(para, diag, varK, varT)
        # w = DiagTree.evalNaive(diag, varK, varT, evalFakePropagator)
        w = evalDiagTree!(diag, varK, varT, evalFakePropagator)
        # plot_tree(diag, maxdepth = 7)
        factor = (1 / (2π)^para.loopDim)^para.innerLoopNum
        num = w / factor
        @test num * (-1)^(para.innerLoopNum) ≈ Parquet.Benchmark.count_sigma_G2v(para.innerLoopNum, para.spin)
    end

    ##################  G^2*v expansion #########################################
    for l = 1:4
        # ret = getSigma(l, spin = 1, isFermi = false, filter = [Builder.Girreducible,])
        # testDiagramNumber(ret...)
        ret = getSigma(l, spin = 2, isFermi = false, filter = [Girreducible,])
        testDiagramNumber(ret...)
    end

    # para, diag, varK, varT = getSigma(1, spin = 2, isFermi = false, filter = [Builder.NoFock,], subdiagram = true)
    # @test isempty(diag.root)

end

@testset "Green" begin
    function buildG(loopNum, extT; Kdim = 3, spin = 2, interactionTauNum = 1, filter = [], isFermi = true)
        para = GenericPara(
            diagType = GreenDiag,
            loopDim = Kdim,
            hasTau = true,
            innerLoopNum = loopNum,
            isFermi = isFermi,
            spin = spin,
            filter = filter,
            interaction = [Interaction(ChargeCharge, Instant),]
        )
        extK = zeros(para.totalLoopNum)
        extK[1] = 1.0
        G = Parquet.green(para, extK, extT)
        return G
    end
    # diag, Gidx = buildG(2, [1, 2], 3; filter = [])
    # DiagTree.showTree(diag, Gidx)

    # If G is irreducible, then only loop-0 G exist for main diagram, and no G exist for subdiagram
    G = buildG(0, [1, 2]; filter = [Girreducible,])
    @test G isa Diagram
    G = buildG(1, [1, 2]; filter = [Girreducible,])
    @test isnothing(G)
    G = buildG(2, [1, 2]; filter = [Girreducible,])
    @test isnothing(G)

    # If Fock diagram is not allowed, then one-loop G diagram should not be exist for subdiagram
    G = buildG(0, [1, 2]; filter = [NoFock,])
    @test G isa Diagram
    G = buildG(1, [1, 2]; filter = [NoFock,])
    @test isnothing(G)
    G = buildG(2, [1, 2]; filter = [NoFock,]) #high order subdiagram is allowed
    @test G isa Diagram

end

@testset "Parquet Vertex3" begin
    function getGamma3(loopNum; Kdim = 3, spin = 2, interactionTauNum = 1, filter = [Girreducible, Proper,], isFermi = true, subdiagram = false)
        println("LoopNum =$loopNum Vertex3 Test")

        para = GenericPara(
            diagType = Ver3Diag,
            loopDim = Kdim,
            innerLoopNum = loopNum,
            isFermi = isFermi,
            hasTau = true,
            filter = filter,
            interaction = [Interaction(ChargeCharge, Instant),]
        )

        K0 = zeros(para.totalLoopNum)
        KinL, Q = deepcopy(K0), deepcopy(K0)
        Q[1] = 1
        KinL[2] = 1
        legK = [Q, KinL]

        varK = rand(Kdim, para.totalLoopNum)
        varT = [rand() for i in 1:para.totalTauNum]

        #################### DiagTree ####################################
        vertex3 = Parquet.vertex3(para, legK)
        diag = mergeby(vertex3)
        # print_tree(diag.diagram[1])

        return para, diag.diagram[1], varK, varT
    end

    function testDiagramNumber(para, diag, varK, varT)
        # w = DiagTree.evalNaive(diag, varK, varT, evalFakePropagator)
        w = evalDiagTree!(diag, varK, varT, evalFakePropagator)
        # plot_tree(diag, maxdepth = 9)
        factor = (1 / (2π)^para.loopDim)^para.innerLoopNum
        num = w / factor
        @test num * (-1)^(para.innerLoopNum) ≈ Parquet.Benchmark.count_ver3_G2v(para.innerLoopNum, para.spin)
    end

    ##################  G^2*v expansion #########################################
    for l = 1:3
        # ret = getSigma(l, spin = 1, isFermi = false, filter = [Builder.Girreducible,])
        # testDiagramNumber(ret...)
        ret = getGamma3(l, isFermi = false, filter = [Girreducible, Proper])
        testDiagramNumber(ret...)
    end

    # para, diag, varK, varT = getSigma(1, spin = 2, isFermi = false, filter = [Builder.NoFock,], subdiagram = true)
    # @test isempty(diag.root)

end

@testset "Parquet Polarization" begin
    function getPolar(loopNum; Kdim = 3, spin = 2, interactionTauNum = 1, filter = [Girreducible,], isFermi = true, subdiagram = false)
        println("LoopNum =$loopNum Polarization Test")

        para = GenericPara(
            diagType = PolarDiag,
            loopDim = Kdim,
            innerLoopNum = loopNum,
            isFermi = isFermi,
            hasTau = true,
            filter = filter,
            interaction = [Interaction(ChargeCharge, Instant),]
        )

        Q = zeros(para.totalLoopNum)
        Q[1] = 1

        varK = rand(Kdim, para.totalLoopNum)
        varT = [rand() for i in 1:para.totalTauNum]

        #################### DiagTree ####################################
        diag = Parquet.polarization(para, Q)
        # print_tree(diag.diagram[1])
        return para, diag, varK, varT
    end

    ##################  G^2*v expansion #########################################
    for l = 1:4
        para, diag, varK, varT = getPolar(l, isFermi = false, filter = [Girreducible,])
        diag = mergeby(diag).diagram[1]
        w = evalDiagTree!(diag, varK, varT, evalFakePropagator)
        factor = (1 / (2π)^para.loopDim)^para.innerLoopNum
        num = w / factor
        # println(num * para.spin)
        @test num * para.spin * (-1)^(para.innerLoopNum - 1) ≈ Parquet.Benchmark.count_polar_G2v(para.innerLoopNum, para.spin)
    end

    ##################  g^2*v expansion #########################################
    for l = 1:4
        para, diag, varK, varT = getPolar(l, isFermi = false, filter = [NoFock,])
        diag = mergeby(diag).diagram[1]
        w = evalDiagTree!(diag, varK, varT, evalFakePropagator)
        factor = (1 / (2π)^para.loopDim)^para.innerLoopNum
        num = w / factor
        # println(num * para.spin)
        @test num * para.spin * (-1)^(para.innerLoopNum - 1) ≈ Parquet.Benchmark.count_polar_g2v_noFock(para.innerLoopNum, para.spin)
    end

    ##################  g^2*v expansion for the upup polarization #########################################
    for l = 1:4
        para, diag, varK, varT = getPolar(l, isFermi = false, filter = [NoFock,])
        w = evalDiagTree!(diag.diagram[1], varK, varT, evalFakePropagator)
        factor = (1 / (2π)^para.loopDim)^para.innerLoopNum
        num = w / factor
        # println(num * para.spin)
        @test num * para.spin * (-1)^(para.innerLoopNum - 1) ≈ Parquet.Benchmark.count_polar_g2v_noFock_upup(para.innerLoopNum, para.spin)
    end
end