text:
push!(LOAD_PATH, pwd())
include("startup.jl")
using NetCDF

# function definitions start ##########################################################################################
function createBCData(filepath :: String, bcheader)
filefilter = r"\.CSV$"
files = filter(s->occursin(filefilter, s), readdir(filepath))
nFiles = size(files,1)
for j=1:nFiles
    print(j)
    totalPath = joinpath(filepath, files[j])
    if j < nFiles
      totalPrecachePath = joinpath(filepath, files[j+1])
    else
      totalPrecachePath = ""
    end
    bc = readdlm(totalPath,',')
    bcdatetime = Dates.Date.(bc[:,2],"dd-u-yy").+Dates.Year(2000) .+ Dates.Time.(Dates.DateTime.(bc[:,3],"H:M"))
    bc[findall(x -> typeof(x) != Float64 && typeof(x) != Int64 && typeof(x) != DateTime, bc)] .= NaN
    bcunixtime = Dates.datetime2unix.(bcdatetime)
    bcnew = hcat(bcunixtime, bcdatetime)
    bcnew = hcat(bcnew, bc[:,4:10])
    bcheader = vcat(bcheader,bcnew)
  end
  return bcheader
end

function createACSMData(ACSMf :: String)
	(acsm, acsmHeader) = readdlm(ACSMf,'\t', header = true)
	ACSMdatetime = Dates.DateTime.(acsm[:,2],"d.m.Y H:M")+Hour(1) # time-> unix2datetime - 1h diff
	ACSMunixtime = Dates.datetime2unix.(ACSMdatetime)
	return (hcat(["unixTime"], acsmHeader), hcat(ACSMunixtime, acsm))
end

function createMeteoDataCumbre(filepath :: String)
    (meteo, meteoHeader) = readdlm(filepath, header = true) # time in local time
    MeteoDatetime = DateTime.(meteo[:,1], meteo[:,2], meteo[:,3], meteo[:,4], meteo[:,5])
    MeteoUnixtime = Dates.datetime2unix.(MeteoDatetime)
  return (hcat(["unixTime"], meteoHeader), hcat(MeteoUnixtime, meteo))
end

# function definitions end ********************************************************************************************

# loading of data and definition of variables #########################################################################

if plotMeteo == true && !(@isdefined Temp)
	println("  -> load Meteo Data")
	meteoFile = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/CHC_cumbre_15M_2017-2018_saltena_isabel/CHC_cumbre_15M_2017-2018_saltena_noHeader.txt"
	(meteohead, meteo) = createMeteoDataCumbre(meteoFile)
	meteo[findall(x -> x == -999.9, meteo)] .= NaN
	meteo_UT = meteo[:,1]
	meteo_DT = Dates.unix2datetime.(meteo_UT)
	SWd = meteo[:,8]
	WS = meteo[:,13]
	WD = meteo[:,14]
	WDstdv = meteo[:,15]
	Temp = meteo[:,16]
	RH = meteo[:,19]
	#AH = (exp.(77.34 .- (7235 ./ Temp)- 8.2 .* log.(Temp) .+ 0.005711.*Temp) .* (RH./100)) ./ (Temp.*461.5)

	meteoS = readdlm("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/meteo_chc/results/meteoData.txt")
	stationUT = convert(Array{Float64,1}, meteoS[2:end,1]) #.-3600*4)
	stationDT = Dates.unix2datetime.(stationUT)
	TEMPstation = Float64.(meteoS[2:end,4])
	RHstation = Float64.(meteoS[2:end,6])
	PRESSUREstation = Float64.(meteoS[2:end,5])
	AHstation = Array{Float64}(meteoS[2:end,9])*1000	#

	if selectedTime == true
		SWd = SWd[DatasetFunctions.timeRange(toi[1],toi[2],meteo_DT)]
		WS = WS[DatasetFunctions.timeRange(toi[1],toi[2],meteo_DT)]
		WD = WD[DatasetFunctions.timeRange(toi[1],toi[2],meteo_DT)]
		WDstdv = WDstdv[DatasetFunctions.timeRange(toi[1],toi[2],meteo_DT)]
		Temp = Temp[DatasetFunctions.timeRange(toi[1],toi[2],meteo_DT)]
		RH = RH[DatasetFunctions.timeRange(toi[1],toi[2],meteo_DT)]
		meteo_UT = meteo_UT[DatasetFunctions.timeRange(toi[1],toi[2],meteo_DT)]
		meteo_DT = meteo_DT[DatasetFunctions.timeRange(toi[1],toi[2],meteo_DT)]

		TEMPstation = TEMPstation[(toi[1] .< stationDT .< toi[2])]
		RHstation = RHstation[(toi[1] .< stationDT .< toi[2])]
		AHstation = AHstation[(toi[1] .< stationDT .< toi[2])]
		PRESSUREstation = PRESSUREstation[(toi[1] .< stationDT .< toi[2])]
		stationUT = stationUT[(toi[1] .< stationDT .< toi[2])]
		stationDT = stationDT[(toi[1] .< stationDT .< toi[2])]
	end
	if interpolateTraces == true
		SWd = InterpolationFunctions.interpolate(intpolT, meteo_UT, SWd)
		WS = InterpolationFunctions.interpolate(intpolT, meteo_UT, WS)
		WD = InterpolationFunctions.interpolate(intpolT, meteo_UT, WD)
		WDstdv = InterpolationFunctions.interpolate(intpolT, meteo_UT, WDstdv)
		Temp = InterpolationFunctions.interpolate(intpolT, meteo_UT, Temp)
		RH = InterpolationFunctions.interpolate(intpolT, meteo_UT, RH)
		TEMPstation = InterpolationFunctions.interpolate(intpolT, stationUT, TEMPstation)
		RHstation = InterpolationFunctions.interpolate(intpolT, stationUT, RHstation)
		AHstation = InterpolationFunctions.interpolate(intpolT, stationUT, AHstation)
		PRESSUREstation = InterpolationFunctions.interpolate(intpolT, stationUT, PRESSUREstation)
	end
end

if plotBC == true && !(@isdefined BC_sum)
	println("  -> load Black Carbon")
	BCfp = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/BC_data_SALTENA/Chacaltaya/"
	BChead = ["unixTime" "DateTime" "conc370nm" "conc470nm" "conc520nm" "conc590nm" "conc660nm" "conc880nm" "conc950nm";]
	BC = createBCData(BCfp, BChead)[2:end,:]
	BC_UT = BC[:,1]
	BC_DT = DateTime.(BC[:,2])
	BC_sum = DatasetFunctions.nansum(BC[:,3:9],2)
	BC_UT = BC_UT[DatasetFunctions.timeRange(toi[1],toi[2],BC_DT)]
	BC_sum = BC_sum[DatasetFunctions.timeRange(toi[1],toi[2],BC_DT)]
	BC_DT = BC_DT[DatasetFunctions.timeRange(toi[1],toi[2],BC_DT)]
	BC_sum = InterpolationFunctions.interpolate(intpolT, BC_UT, BC_sum)
end

if plotACSM == true && !(@isdefined ACSM_org)
	println("  -> load ACSM Data")
	ACSMfile = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/ACSM/ACSM_MarMay2018_CHC_Wiebke-1.csv"
	(ACSMhead, ACSM) = createACSMData(ACSMfile)
	ACSM_UT = ACSM[:,1] .- 5*3600 	#UTC to local
	ACSM_DT = Dates.unix2datetime.(ACSM_UT)
	ACSM_sulf = ACSM[:,4]
	ACSM_org = ACSM[:,5]
	ACSM_nitr = ACSM[:,6]
	ACSM_ammon = ACSM[:,7]
	ACSM_chlor = ACSM[:,8]
	ACSM_all = ACSM[:,4:8]
	ACSM_org = ACSM_org[DatasetFunctions.timeRange(toi[1],toi[2],ACSM_DT)]
	ACSM_sulf = ACSM_sulf[DatasetFunctions.timeRange(toi[1],toi[2],ACSM_DT)]
	ACSM_nitr = ACSM_nitr[DatasetFunctions.timeRange(toi[1],toi[2],ACSM_DT)]
	ACSM_ammon = ACSM_ammon[DatasetFunctions.timeRange(toi[1],toi[2],ACSM_DT)]
	ACSM_chlor = ACSM_chlor[DatasetFunctions.timeRange(toi[1],toi[2],ACSM_DT)]
	ACSM_all = ACSM_all[DatasetFunctions.timeRange(toi[1],toi[2],ACSM_DT),:]
	ACSM_UT = ACSM_UT[DatasetFunctions.timeRange(toi[1],toi[2],ACSM_DT)]
	ACSM_DT = ACSM_DT[DatasetFunctions.timeRange(toi[1],toi[2],ACSM_DT)]
	if length(ACSM_org) > 0
		ACSM_org = InterpolationFunctions.interpolate(intpolT, ACSM_UT, ACSM_org)
		ACSM_sulf = InterpolationFunctions.interpolate(intpolT, ACSM_UT, ACSM_sulf)
		ACSM_nitr = InterpolationFunctions.interpolate(intpolT, ACSM_UT, ACSM_nitr)
		ACSM_ammon = InterpolationFunctions.interpolate(intpolT, ACSM_UT, ACSM_ammon)
		ACSM_chlor = InterpolationFunctions.interpolate(intpolT, ACSM_UT, ACSM_chlor)
		ACSM_all = InterpolationFunctions.interpolateMatrix_1D(intpolT, ACSM_UT, ACSM_all)
	else plotACSM = false
	end
end

if plotSMPS == true && !(@isdefined SMPS_binconc)	# can be done faster by loading old interpolated dataset, if available!
	println("  -> load SMPS Data")
	SMPSfile = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/Level_1_data_particle_instruments/SMPS_chacaltaya_data/chacaltaya_new.sum"
	SMPS = readdlm(SMPSfile)
	SMPS_DT = DatasetFunctions.matlab2datetime(SMPS[2:end,1])
	SMPS_UT = DatasetFunctions.matlab2unixtime(SMPS[2:end,1])
	SMPS_binsize = SMPS[1,3:end]
	SMPS_totconc = SMPS[2:end,2]
	SMPS_binconc = SMPS[2:end,3:end]
	SMPS_binconc = SMPS_binconc[DatasetFunctions.timeRange(toi[1],toi[2],SMPS_DT),:]
	SMPS_totconc = SMPS_totconc[DatasetFunctions.timeRange(toi[1],toi[2],SMPS_DT)]
	SMPS_UT = SMPS_UT[DatasetFunctions.timeRange(toi[1],toi[2],SMPS_DT)]
	SMPS_DT = SMPS_DT[DatasetFunctions.timeRange(toi[1],toi[2],SMPS_DT)]
	SMPS_binconc = InterpolationFunctions.interpolateMatrix_1D(intpolT, SMPS_UT, SMPS_binconc)
	SMPS_totconc = InterpolationFunctions.interpolate(intpolT, SMPS_UT, SMPS_totconc)
end

if plotNO3cims == true && !(@isdefined Nitrophenol)
	println("  -> load NO3cims Data")
	(no3cims, no3cimshead) = readdlm("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/George_Data_CI.csv", ',',header = true)
	no3cims[1,[2,3,4,6]] .= NaN
	no3cims_DT = Dates.DateTime.(no3cims[:,1],"d/m/YTH:M:S") .+ Dates.Hour(1)	# weird time to local time
	no3cims_UT = Dates.datetime2unix.(no3cims_DT)
	no3cims_SA_DT = Dates.DateTime.(no3cims[:,5][findall(x -> x != "", no3cims[:,5])],"d/m/YTH:M:S")
	no3cims_SA_UT = Dates.datetime2unix.(no3cims_SA_DT)
	(no3cimsHOMs, no3cimsheadHOMs) = readdlm("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/SA_HOM/May_SA_Total HOM_betterforReaddlm.csv", '\t',header = true)
	HOMs = no3cimsHOMs[:,3]
	HOMs_UT = no3cimsHOMs[:,1] .- 4*3600
	Nitrophenol = no3cims[2:end,2]
	SO5trace = no3cims[2:end,3]
	MSAtrace = no3cims[2:end,4]
	SAtrace = no3cims[findall(x -> x != "", no3cims[:,5]),6]
	SAtrace[findall(x -> typeof(x) == SubString{String}, SAtrace)] .= NaN
	SAtrace = SAtrace./2.47e4
	if selectedTime == true
		SO5trace = SO5trace[DatasetFunctions.timeRange(toi[1],toi[2],no3cims_DT)]
		MSAtrace = MSAtrace[DatasetFunctions.timeRange(toi[1],toi[2],no3cims_DT)]
		Nitrophenol = Nitrophenol[DatasetFunctions.timeRange(toi[1],toi[2],no3cims_DT)]
		SAtrace = SAtrace[DatasetFunctions.timeRange(toi[1],toi[2],no3cims_SA_DT)]
		no3cims_UT = no3cims_UT[DatasetFunctions.timeRange(toi[1],toi[2],no3cims_DT)]
		no3cims_DT = no3cims_DT[DatasetFunctions.timeRange(toi[1],toi[2],no3cims_DT)]
		no3cims_SA_UT = no3cims_SA_UT[DatasetFunctions.timeRange(toi[1],toi[2],no3cims_SA_DT)]
		no3cims_SA_DT = no3cims_SA_DT[DatasetFunctions.timeRange(toi[1],toi[2],no3cims_SA_DT)]
	end
	if interpolateTraces == true && selectedTime == true
		MSAtrace = InterpolationFunctions.interpolate(intpolT, no3cims_UT, MSAtrace)
		Nitrophenol = InterpolationFunctions.interpolate(intpolT, no3cims_UT, Nitrophenol)
		SO5trace = InterpolationFunctions.interpolate(intpolT, no3cims_UT, SO5trace)
		SAtrace = InterpolationFunctions.interpolate(intpolT, no3cims_SA_UT, SAtrace)
		HOMs = InterpolationFunctions.interpolate(intpolT, HOMs_UT, HOMs)
	end
end

if plotSRR == true  && !(@isdefined SRRs)
	println("  -> load Source Receptor Relationships")
	SRRfilename = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/Flexpart/v03/cluster_series_v3.nc"
	# SRRdata = ncread(SRRfilename, "conc_lab_nc18", z_column="ALL", normalized=1)
	SRRdata = ncread(SRRfilename, "conc_lab_nc18")[:,1,2,:]
	SRRheader = ncread(SRRfilename, "lab_nc18")
	SRRheader = ncread(SRRfilename, "conc_all")
	SRRutctimes = Dates.Hour.(ncread(SRRfilename, "releases")) .+ Dates.DateTime(2017,12,06)
	SRR_DT = SRRutctimes .- Dates.Hour(4)	# UTC -> localTime
	SRR_UT = Dates.datetime2unix.(SRR_DT)
	SRRage = ncread(SRRfilename, "age_all")[:,1]
	SRRheight = ncread(SRRfilename, "ZSL_all")[:,1]
	SRRheight_AG = ncread(SRRfilename, "ZGL_all")[:,1]
	SRRclockdir = ncread(SRRfilename, "clock_dir_all")[:,1]
	SRRconcAll = ncread(SRRfilename, "conc_all")[:,1,1]
	SRRs_06 = ncread(SRRfilename, "conc_lab_nc06")[:,1,2,:]
	SRRs = InterpolationFunctions.interpolateMatrix_1D(intpolT, SRR_UT, SRRdata)
	SRRs_06 = InterpolationFunctions.interpolateMatrix_1D(intpolT, SRR_UT, SRRs_06)
	SRRage = InterpolationFunctions.interpolate(intpolT, SRR_UT, SRRage)
	SRRheight = InterpolationFunctions.interpolate(intpolT, SRR_UT, SRRheight)
	SRRheight_AG = InterpolationFunctions.interpolate(intpolT, SRR_UT, SRRheight_AG)
	SRRclockdir = InterpolationFunctions.interpolate(intpolT, SRR_UT, SRRclockdir)
	SRRconcAll = InterpolationFunctions.interpolate(intpolT, SRR_UT, SRRconcAll)
	SRR_z0 = ncread(SRRfilename, "conc_all")[:,:,2][:,3]
	SRR_z0 = InterpolationFunctions.interpolate(intpolT, SRR_UT, SRR_z0)
	SRR_bl = ncread(SRRfilename, "conc_all")[:,:,2][:,2]
	SRR_bl = InterpolationFunctions.interpolate(intpolT, SRR_UT, SRR_bl)
	SRRnames = SRRheader[1:end]
	SRR06names = ncread(SRRfilename, "lab_nc06")
end

function createTotalPTR3data()
	fp1 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3OnoRF/results/calibratedWithHexanone/"
	#(traces1, header1) = readdlm("$(fp1)ptr3tracesCORR.csv", header = true)
	(traces1, header1) = readdlm("$(fp1)ptr3traces.csv", header = true)
	(peaks1, peakheader1) = readdlm("$(fp1)ptr3compositions.txt", header = true)
	times1 = traces1[:,1]
	traces1 = traces1[:,2:end]
	#fp2 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/WetAtmo/results/calibrated_drySulfMasslist/"
	#(traces2, header2) = readdlm("$(fp2)_ptr3traces-hexanone_drySulfML.csv", header = true)
	#(peaks2, peakheader2) = readdlm("$(fp2)_ptr3compositions-hexanone_drySulfML.txt", header = true)
	#times2 = traces2[:,1]
	#traces2 = traces2[:,2:end]
	#fp3 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/DryAtmo/results/calibratedWithHexanone/"
	#(traces3, header3) = readdlm("$(fp3)ptr3traces.csv", header = true)
	#(peaks3, peakheader3) = readdlm("$(fp3)ptr3compositions.txt", header = true)
	#(traces3sulf, header3sulf) = readdlm("$(fp3)Sulfur_ptr3traces.csv", header = true)
	#(peaks3sulf, peakheader3sulf) = readdlm("$(fp3)Sulfur_ptr3compositions.txt", header = true)
	#times3 = traces3[:,1]
	#traces3 = traces3[:,2:end]
	#traces3sulf = traces3sulf[:,2:end]
	fp2 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/400Vpp_all/results/"
	(traces2, header2) = readdlm("$(fp2)_ptr3traces-hexanone_all400V.csv", header = true)
	(peaks2, peakheader2) = readdlm("$(fp2)_ptr3compositions-hexanone_all400V.txt", header = true)
	times2 = traces2[2:end,1]
	traces2 = traces2[2:end,2:end]
	fp4 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF600Vpp_final/results/"
	(traces4, header4) = readdlm("$(fp4)_ptr3traces-hexanone_drySulfML.csv", header = true)
	(peaks4, peakheader4) = readdlm("$(fp4)_ptr3compositions-hexanone_drySulfML.txt", header = true)
	times4 = traces4[:,1]
	times4 = times4[findall(!(in(times2)), times4)]
	traces4 = traces4[findall(!(in(times2)), times4),2:end]

	allPeaks = vcat(peaks1, peaks2, peaks4) # peaks3, peaks3sulf, peaks4)
	allPeaks = unique(allPeaks; dims = 1)
	allPeaks = allPeaks[sortperm(allPeaks[:,1]),:]
	allTraces = Array{Float64}(undef, (size(traces1)[1] + size(traces2)[1] + size(traces4)[1]), size(allPeaks)[1]) # + size(traces3)[1]
	allTraces .= NaN
	allTimes = vcat(times1, times2, times4) # times3,
	timemask1 = findall(in(times1), allTimes)
	timemask2 = findall(in(times2), allTimes)
	#timemask3 = findall(in(times3), allTimes)
	timemask4 = findall(in(times4), allTimes)
	peakmask1 = findall(in(peaks1[:,1]), allPeaks[:,1])
	peakmask2 = findall(in(peaks2[:,1]), allPeaks[:,1])
	#peakmask3 = findall(in(peaks3[:,1]), allPeaks[:,1])
	#peakmask3sulf = findall(in(peaks3sulf[:,1]), allPeaks[:,1])
	peakmask4 = findall(in(peaks4[:,1]), allPeaks[:,1])
	allTraces[timemask1,peakmask1] = traces1
	allTraces[timemask2,peakmask2] = traces2
	allTraces[timemask4[:],peakmask4[:]] = traces4
	#allTraces[timemask3[:],peakmask3sulf[:]] = traces3sulf
	#allTraces[timemask3[:],peakmask3[:]] = traces3
	maskingTimes = (allTimes .< Dates.datetime2unix(DateTime(2018,5,2))) .+
			(Dates.datetime2unix(DateTime(2018,5,3,14,40)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,5,14,20))) .+
			(Dates.datetime2unix(DateTime(2018,5,3,14,40)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,5,14,20))) .+
			(Dates.datetime2unix(DateTime(2018,5,8,9,30)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,8,14,40))) .+
			(Dates.datetime2unix(DateTime(2018,5,9,6,30)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,9,9,45))) .+
			(Dates.datetime2unix(DateTime(2018,5,9,23,30)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,10,0,45))) .+
			(Dates.datetime2unix(DateTime(2018,5,17,3,00)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,17,11))) .+
			(Dates.datetime2unix(DateTime(2018,5,22,8,00)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,23,8))) .+
			(Dates.datetime2unix(DateTime(2018,5,24,16,00)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,28,11))) .+
			(Dates.datetime2unix(DateTime(2018,5,30,0,00)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,30,23,59)))
	allTraces[findall(x -> x > 0, maskingTimes), :]	.= NaN
	f = open("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/TotalTraces_new2020.csv", "w")
	writedlm(f, hcat(vcat(["unixTimes"], allTimes), vcat(permutedims(allPeaks[:,2]), allTraces)))
	close(f)
	g = open("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/TotalPeaks_new2020.csv", "w")
	writedlm(g, vcat(peakheader1, allPeaks))
	close(g)
	return allPeaks, allTimes, allTraces
end

function createTotalPTR3dataNonOxidized()
	fp1 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3OnoRF/results/calibratedWithAcetonitril/"
	#(traces1, header1) = readdlm("$(fp1)ptr3tracesCORR.csv", header = true)
	(traces1, header1) = readdlm("$(fp1)ptr3traces.csv", header = true)
	(peaks1, peakheader1) = readdlm("$(fp1)ptr3compositions.txt", header = true)
	times1 = traces1[:,1]
	traces1 = traces1[:,2:end]
	fp2 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/WetAtmo/results/calibrated_drySulfMasslist/"
	(traces2, header2) = readdlm("$(fp2)_ptr3traces-acetonitrile_drySulfML.csv", header = true)
	(peaks2, peakheader2) = readdlm("$(fp2)_ptr3compositions-acetonitrile_drySulfML.txt", header = true)
	times2 = traces2[:,1]
	traces2 = traces2[:,2:end]
	fp3 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/DryAtmo/results/calibratedWithAcetonitril/"
	(traces3, header3) = readdlm("$(fp3)ptr3traces.csv", header = true)
	(peaks3, peakheader3) = readdlm("$(fp3)ptr3compositions.txt", header = true)
	(traces3sulf, header3sulf) = readdlm("$(fp3)Sulfur_ptr3traces.csv", header = true)
	(peaks3sulf, peakheader3sulf) = readdlm("$(fp3)Sulfur_ptr3compositions.txt", header = true)
	times3 = traces3[:,1]
	traces3 = traces3[:,2:end]
	traces3sulf = traces3sulf[:,2:end]
	fp4 = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF600Vpp_final/results/"
	(traces4, header4) = readdlm("$(fp4)_ptr3traces-acetonitrile_drySulfML.csv", header = true)
	(peaks4, peakheader4) = readdlm("$(fp4)_ptr3compositions-acetonitrile_drySulfML.txt", header = true)
	times4 = traces4[:,1]
	times4 = times4[findall(!(in(times3)), times4)]
	traces4 = traces4[findall(!(in(times3)), times4),2:end]

	allPeaks = vcat(peaks1, peaks2, peaks3, peaks3sulf, peaks4)
	allPeaks = unique(allPeaks; dims = 1)
	allPeaks = allPeaks[sortperm(allPeaks[:,1]),:]
	allTraces = Array{Float64}(undef, (size(traces1)[1] + size(traces2)[1] + size(traces3)[1]) + size(traces4)[1], size(allPeaks)[1])
	allTraces .= NaN
	allTimes = vcat(times1, times2, times3, times4)
	timemask1 = findall(in(times1), allTimes)
	timemask2 = findall(in(times2), allTimes)
	timemask3 = findall(in(times3), allTimes)
	timemask4 = findall(in(times4), allTimes)
	peakmask1 = findall(in(peaks1[:,1]), allPeaks[:,1])
	peakmask2 = findall(in(peaks2[:,1]), allPeaks[:,1])
	peakmask3 = findall(in(peaks3[:,1]), allPeaks[:,1])
	peakmask3sulf = findall(in(peaks3sulf[:,1]), allPeaks[:,1])
	peakmask4 = findall(in(peaks4[:,1]), allPeaks[:,1])
	allTraces[timemask1,peakmask1] = traces1
	allTraces[timemask2,peakmask2] = traces2
	allTraces[timemask4[:],peakmask4[:]] = traces4
	allTraces[timemask3[:],peakmask3sulf[:]] = traces3sulf
	allTraces[timemask3[:],peakmask3[:]] = traces3
	maskingTimes = (allTimes .< Dates.datetime2unix(DateTime(2018,5,2))) .+
			(Dates.datetime2unix(DateTime(2018,5,3,14,40)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,5,14,20))) .+
			(Dates.datetime2unix(DateTime(2018,5,3,14,40)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,5,14,20))) .+
			(Dates.datetime2unix(DateTime(2018,5,8,9,30)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,8,14,40))) .+
			(Dates.datetime2unix(DateTime(2018,5,9,6,30)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,9,9,45))) .+
			(Dates.datetime2unix(DateTime(2018,5,9,23,30)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,10,0,45))) .+
			(Dates.datetime2unix(DateTime(2018,5,17,3,00)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,17,11))) .+
			(Dates.datetime2unix(DateTime(2018,5,22,8,00)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,23,8))) .+
			(Dates.datetime2unix(DateTime(2018,5,24,16,00)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,28,11))) .+
			(Dates.datetime2unix(DateTime(2018,5,30,0,00)) .< allTimes .< Dates.datetime2unix(DateTime(2018,5,30,23,59)))
	allTraces[findall(x -> x > 0, maskingTimes), :]	.= NaN
	f = open("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/TotalTraces_AcNitcalib.csv", "w")
	writedlm(f, hcat(vcat(["unixTimes"], allTimes), vcat(permutedims(allPeaks[:,2]), allTraces)))
	close(f)
	g = open("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/TotalPeaks_AcNitcalib.csv", "w")
	writedlm(g, vcat(peakheader1, allPeaks))
	close(g)
	return allPeaks, allTimes, allTraces
end

# variable definitions end ********************************************************************************************


text:
"""
    NISTConst{T} <: AbstractPhysicsConst{T}

NIST CODATA Physics Constants (2014).
"""
struct NISTConst{T} <: AbstractPhysicsConst{T}
    quantity::String
    value::T
    uncertainty::T
    unit::String
end

value(x::NISTConst) = x.value
eltype(::Type{NISTConst{T}}) where T = T

function NISTConst(d::Dict{String})
    value = Compat.replace(d["Value"], " " => "")
    uncertainty = 0.0

    m = match(r"([0-9]+\.[0-9]+)(\.{3})(.*)", value)
    if m != nothing
        value = Meta.parse(m[1] * m[3])
    elseif d["Uncertainty"] == "(exact)"
        value = Meta.parse(value)
    else
        value = Meta.parse(value)
        uncertainty = Meta.parse(Compat.replace(d["Uncertainty"], " " => ""))
    end

    NISTConst(d["Quantity "], promote(value, uncertainty)..., d["Unit"])
end

json(c::NISTConst) = throw(ErrorException("NISTConst should not be serialized"))

function parse(::Type{NISTConst}, file::DataFile{JSONFormat})
    raw = parse(file)["constant"]
    data = Dict{String, Any}()
    for each in raw
        data[each["Quantity "]] = NISTConst(each)
    end
    return data
end

ext: jl


text:
import Statistics, DataFrames, StatsBase

"""
    standardize(X; robust::Bool=false)

Standardize (scale and reduce, Z-score) X so that the variables are expressed in terms of standard deviation (*i.e.*, mean = 0, SD = 1).

# Arguments
- `X`: Array or DataFrame.
- `robust::Bool`: If true, the standardization will be based on `median` and `mad` instead of `mean` and `sd` (default).

!!! note

    **Ideas / help required:**
    - Deal with missing values (See [#4](https://github.com/neuropsychology/Psycho.jl/issues/4))

# Examples
```jldoctest
standardize([1, 2, 3])

# output

3-element Array{Float64,1}:
 -1.0
  0.0
  1.0
```
"""
function standardize end

"""
    standardize!(X; robust::Bool=false)

In-place version of [`standardize`](@ref).
"""
function standardize! end

function standardize(X::Array{<:Number}; robust::Bool=false)

    if robust == false
        z = (X .- Statistics.mean(X, dims=1)) ./ Statistics.std(X, dims=1)
    else
        z = (X .- Statistics.median(X, dims=1)) ./ mapslices(x -> StatsBase.mad(x, normalize=true), X, dims=1)
    end
    return z
end

function standardize(X::Vector{<:Union{Missing, Number}}; robust::Bool=false)

    complete = collect(skipmissing(X))

    if robust == false
        z = (X .- Statistics.mean(complete)) ./ Statistics.std(complete)
    else
        z = (X .- Statistics.median(complete)) ./ StatsBase.mad(complete, normalize=true)
    end
    return z
end

function standardize!(X::Array{<:Union{Missing, Number}}; robust=false)
    X[:] = standardize(X, robust=robust)
end

function standardize!(X::DataFrames.DataFrame; robust=false)
    for (colname, col) in DataFrames.eachcol(X)
        if eltype(col) <: Union{Missing, Number}
            if length(unique(col)) > 2
                X[colname] = standardize(col)
            end
        end
    end
end

function standardize(X::DataFrames.DataFrame; robust=false)
    df = copy(X)
    standardize!(df; robust=robust)
    return df
end


text:
# validation
export is_povm, is_povm_element

# types
export Measurement
export POVMel, POVM, PVM

"""
    Measurement{T} <: AbstractVector{T}

The abstract type representing a quantum measurement. A measurement consists of a
complete set of possible outcomes.
"""
abstract type Measurement{T} <: AbstractVector{T} end
Base.size( Π :: Measurement) = size(Π.V)
Base.getindex( Π :: Measurement, id::Int) = getindex(Π.V, id...)
Base.setindex!( Π :: Measurement, val, id::Int) = (Π.V[id...] = val)

"""
    is_povm_element(M :: AbstractMatrix; atol=ATOL :: Float64) :: Bool

Returns `true` if matrix `M` satisfies the following constraints:
* `M` is Hermitian
* `M` is positive semi-definite
"""
function is_povm_element(M :: AbstractMatrix; atol=ATOL :: Float64) :: Bool
    if !is_hermitian(M, atol=atol)
        return false
    elseif !is_positive_semidefinite(M, atol=atol)
        return false
    end

    return true
end

"""
    is_povm( Π :: Vector; atol=ATOL :: Float64 ) :: Bool

Returns `true` if `Π` satisfies the following constraints
* The POVM is complete: `sum(Π) == I`
* Each POVM element is hermitian
* Each POVM element positive semi-definite
"""
function is_povm(Π::Vector{<:AbstractMatrix}; atol=ATOL :: Float64) :: Bool
    if !is_complete(Π, atol=atol)
        return false
    elseif !all(M -> is_povm_element(M, atol=atol), Π)
        return false
    end

    return true
end

"""
    POVMel( M :: AbstractMatrix{T<:Number}; atol=ATOL :: Float64) <: Operator{T}

A [`POVM`](@ref) element. A `DomainError` is thrown if matrix `M` is not hermitian
and positive semi-definite within absolute tolerance `atol`.
"""
struct POVMel{T<:Number} <: Operator{T}
    M :: Matrix{T}
    atol :: Float64
    POVMel(
        M :: Matrix{<:Number}; atol=ATOL :: Float64
    ) = is_povm_element(M,atol=atol) ? new{eltype(M)}(M,atol) : throw(DomainError(M, "POVM element M is invalid"))
    POVMel(
        ρ :: State; atol=ATOL :: Float64
    ) = is_povm_element(ρ, atol=atol) ? new{eltype(ρ.M)}(ρ.M,atol) : throw(DomainError(ρ.M, "POVM element M is invalid"))
end
is_povm_element(::POVMel) :: Bool = true
kron(Π_els :: Vararg{POVMel}) = POVMel(kron(map(Π_el -> Π_el.M, Π_els)...))

"""
    POVM( Π :: Vector{Matrix{T}} ) <: Measurement{T}

Positve-operator valued measures (POVMs) represent a general quantum measurement.
Each POVM-element is a hermitian, positive-semidefinite matrix. The sum of all
POVM-elements yields the identity matrix. The constructor, `POVM(Π)` throws a
`DomainError` if the provided array of matrices, `Π` is not a valid POVM.
"""
struct POVM{T<:Number} <: Measurement{POVMel{T}}
    V :: Vector{POVMel{T}}
    atol :: Float64
    POVM(
        Π :: Vector{Matrix{T}};
        atol=ATOL :: Float64
    ) where T <: Number = is_povm(Π, atol=atol) ? new{T}(POVMel.(Π, atol=atol), atol) : throw(DomainError(Π, "povm Π is invalid"))
    POVM(
        Π :: Vector{<:State};
        atol=ATOL :: Float64
    ) = is_povm(Π, atol=atol) ? new{eltype(Π[1].M)}(POVMel.(Π, atol=atol), atol) : throw(DomainError(Π, "povm Π is invalid"))
end
is_povm(::POVM) :: Bool = true

# print out matrix forms when POVM types are displayed
function show(io::IO, mime::MIME{Symbol("text/plain")}, Π :: POVM)
    summary(io, Π)
    for i in 1:length(Π)
        print(io,"\nΠ[",i,"] : ")
        show(io, mime, Π[i])
    end
end

"""
    PVM( Π :: Vector{Vector{T}}; atol=ATOL :: Float64 ) <: Measurement{T}

The concret type for a projector-valued measure. The projectors are represented
as a set of orthonormal basis vectors

A `DomainError` is thrown if `Π` does not contain an orthonormal basis.
"""
struct PVM{T<:Number} <: Measurement{T}
    V :: Vector{Ket{T}}
    atol :: Float64
    PVM(
        Π :: Vector{Vector{T}};
        atol=ATOL :: Float64
    ) where T <: Number = is_orthonormal_basis(Π,  atol=atol) ? new{T}(Ket.(Π, atol=atol), atol) : throw(DomainError(Π, "povm Π is invalid"))
    PVM(
        Π :: Vector{Ket{T}};
        atol=ATOL :: Float64
    ) where  T <: Number = is_orthonormal_basis(Π,  atol=atol) ? new{T}(Π, atol) : throw(DomainError(Π, "povm Π is invalid"))
end
is_orthonormal_basis(::PVM) = true

# TODO: Observable
    # hermitian matrix M
    # M = ∑_j m_j * P_j
    # P_j are rank one projectors
    # m_j are real eigenvalues
    # P_j are the orthonormal eigenvectors of M
# TODO: POVM and PVM conversions
# TODO: kron extensions for measurements


text:
using DBFTables
using Test
using Tables
using DataFrames

test_dbf_path = joinpath(@__DIR__, "test.dbf")
dbf = DBFTables.Table(test_dbf_path)
df = DataFrame(dbf)
row, st = iterate(dbf)

@testset "DBFTables" begin
    @testset "DataFrame indexing" begin
        @test size(df, 1) == 7 # records
        @test size(df, 2) == 6 # fields
        @test df[2, :CHAR] == "John"
        @test df[1, :DATE] == "19900102"
        @test df[3, :BOOL] == false
        @test df[1, :FLOAT] == 10.21
        @test df[2, :NUMERIC] == 12.21
        @test df[3, :INTEGER] == 102
    end

    @testset "missing entries" begin
        @test ismissing(df[4, :BOOL])
        @test ismissing(df[5, :FLOAT])
        @test ismissing(df[6, :NUMERIC])
        @test ismissing(df[7, :INTEGER])
    end

    @testset "header" begin
        h = DBFTables.Header(open(test_dbf_path))
        @test h.version == 3
        @test h.last_update == "20140806"
        @test h.records == 7
        @test length(h.fields) == 6
    end

    @testset "show" begin
        @test sprint(show, row) === sprint(show, NamedTuple(row))
        # use replace to update to julia 1.4 union printing
        @test replace(sprint(show, dbf), r"\} +" => "}") ===
              "DBFTables.Table with 7 rows and 6 columns\nTables.Schema:\n :CHAR     Union{Missing, String}\n :DATE     Union{Missing, String}\n :BOOL     Union{Missing, Bool}\n :FLOAT    Union{Missing, Float64}\n :NUMERIC  Union{Missing, Float64}\n :INTEGER  Union{Missing, $Int}\n"
    end

    @testset "iterate" begin
        @test st === 2
        @test haskey(row, :CHAR)
        @test row.CHAR === "Bob"
        @test row[2] === "19900102"
        @test_throws ArgumentError row.nonexistent_field
        firstrow = (
            CHAR = "Bob",
            DATE = "19900102",
            BOOL = false,
            FLOAT = 10.21,
            NUMERIC = 11.21,
            INTEGER = 100,
        )
        @test NamedTuple(row) === firstrow
        @test row isa DBFTables.Row
        @test row isa Tables.AbstractRow
        @test length(row) === 6
        @test size(row) === (6,)
        @test size(row, 1) === 6
        @test_throws BoundsError size(row, 2)
        @test DBFTables.getrow(row) === 1
        @test DBFTables.gettable(row) === dbf
        @test sum(1 for row in dbf) === 7
        @test sum(1 for cell in row) === 6
        @test propertynames(dbf) == [:CHAR, :DATE, :BOOL, :FLOAT, :NUMERIC, :INTEGER]
        @test propertynames(row) == [:CHAR, :DATE, :BOOL, :FLOAT, :NUMERIC, :INTEGER]
    end

    @testset "column" begin
        @test size(dbf) === (7, 6)
        @test size(dbf, 2) === 6

        @test length(dbf.CHAR) === 7
        @test dbf.CHAR isa Vector{Union{String,Missing}}
        @test dbf.INTEGER isa Vector{Union{Int,Missing}}
        @test_throws ArgumentError row.nonexistent_field
        @test dbf.INTEGER[2] === 101
        @test ismissing(dbf.INTEGER[7])
        @test dbf.CHAR[2] === "John"
        @test ismissing(dbf.CHAR[7])

        @test DBFTables.isdeleted(dbf) isa BitVector
        @test all(.!DBFTables.isdeleted(dbf))
        @test !DBFTables.isdeleted(dbf, 3)
    end

end  # testset "DBFTables"


text:
# ============================================================================ #
# RESULT TYPES
# ============================================================================ #
struct IntrinsicGenerator{T}
    triangle ::Triangle{T}
    path     ::Vector{Int}
    lengths  ::Vector{T}
end

#IntrinsicGenerator(t::Triangle{T}, path, lens) where T =
#    IntrinsicGenerator{T}(t, path, lens)

function Base.show(io::IO, gen::IntrinsicGenerator{T}) where T
    compact = get(io, :compact, false)
    if !compact
        print(io, "IntrinsicGenerator{$T}(")
        showcompact(io, gen.triangle)
        print(io, ", path: $(gen.path), lengths: $(gen.lengths))")
    else
        print(io, "IntrinsicGenerator(")
        showcompact(io, gen.triangle)
        print(io, ", n = $(length(gen.path)), l = $(round(sum(gen.lengths), 5)))")
    end
end

Base.length(gen::IntrinsicGenerator) = sum(gen.lengths)

struct IntrinsicPersistenceResult{T, D, G<:AbstractGraph}
    generators ::Vector{IntrinsicGenerator{T}}
    graph      ::G
    data       ::D

function IntrinsicPersistenceResult(gens  ::Vector{IntrinsicGenerator{T}},
                                    graph ::G,
                                    data  ::D = nothing
                                    ) where {T, G<:AbstractGraph, D}

    new{T, D, G}(gens, graph, data)
end
end

function Base.show(io::IO, res::IntrinsicPersistenceResult{T, D}) where {T, D}
    println(io, "IntrinsicPersistenceResult{$T,$D}:")
    println(io, "  Graph: $(res.graph)")
    print(io,   "  Generators ($(length(res.generators))):")
    for gen in res.generators
        print(io, "\n    ")
        showcompact(io, gen)
    end
end

function Base.convert(::Type{PersistenceBarcode},
                      res::IntrinsicPersistenceResult{T}) where T
    pairs = PersistencePair{T, Void}[]
    for gen in res.generators
        p1 = maximum(gen.lengths)
        p2 = sum(gen.lengths)
        push!(pairs, PersistencePair(p1, p2))
    end
    PersistenceBarcode(PersistencePair{T, Void}[], pairs)
end

# Associate some data with the result. This is used to add points to a result
# that was calculated from a spanner graph.
add_data(res::IntrinsicPersistenceResult{<:Any, Void, <:Any}, data) =
    IntrinsicPersistenceResult(res.generators, res.graph, data)

# ============================================================================ #
# Plots
# ============================================================================ #
# Get x, y, and z coordinates from matrix or arrav of SVectors.
function get_xyz(pts::AbstractMatrix{<:Any}, i=Colon())
    dim = size(pts, 1)
    dim > 3 && warn("Data is in $(dim)d. Only using the first three coordinates.",
                    once = true)
    if dim == 1
        x,y,z = pts[1, i], zeros(size(pts, 2)), zeros(size(pts, 2))
    elseif dim == 2
        x,y,z = pts[1, i], pts[2, i], zeros(size(pts, 2))
    else
        x,y,z = pts[1, i], pts[2, i], pts[3, i]
    end
    Float64.(x), Float64.(y), Float64.(z)
end

function get_xyz(pts::AbstractVector{<:StaticArray}, i=Colon())
    dim = length(eltype(pts))
    dim > 3 && warn("Data is in $(dim)d. Only using the first three coordinates.",
                    once = true)
    pts = pts[i]
    if dim == 1
        x,y,z = getindex.(pts, 1), zeros(length(pts)), zeros(length(pts))
    elseif dim == 2
        x,y,z = getindex.(pts, 1), getindex.(pts, 2), zeros(length(pts))
    else
        x,y,z = getindex.(pts, 1), getindex.(pts, 2), getindex.(pts, 3)
    end
    Float64.(x), Float64.(y), Float64.(z)
end

function get_xyz(::T) where T
    throw(ArgumentError("Can't plot $T!"))
end

@recipe function f(res::IntrinsicPersistenceResult{<:Any}, graph = 0)
    aspect_ratio := 1

    graph ≥ 2 && @series begin
        xs = Float64[]; ys = Float64[]; zs = Float64[]
        for Δ in equilaterals(res.graph.weights, tol = Inf)
            i,j,k = Δ.vertices
            xyz = get_xyz(res.data, [i,j,k,i])
            append!(xs, xyz[1]); push!(xs, NaN)
            append!(ys, xyz[2]); push!(ys, NaN)
            append!(zs, xyz[3]); push!(zs, NaN)
        end

        seriestype := :shape
        alpha := 0.1
        label := "2-skeleton"

        if all(z -> z == 0 || z == NaN, zs)
            xs, ys
        else
            xs, ys, zs
        end
    end

    graph ≥ 1 && @series begin
        xs = Float64[]; ys = Float64[]; zs = Float64[]
        for e in edges(res.graph)
            p1 = src(e)
            p2 = dst(e)
            xyz = get_xyz(res.data, [p1, p2])
            append!(xs, xyz[1]); push!(xs, NaN)
            append!(ys, xyz[2]); push!(ys, NaN)
            append!(zs, xyz[3]); push!(zs, NaN)
        end

        markersize := :none
        seriestype := :path
        linewidth --> 0.5
        if graph ≥ 2
            color := :black
        end
        label := "1-skeleton"

        if all(z -> z == 0 || z == NaN, zs)
            xs, ys
        else
            xs, ys, zs
        end
    end

    @series begin
        seriestype := :scatter
        label --> :points
        markersize --> 0.5

        xs, ys, zs = get_xyz(res.data)
        if all(z -> z == 0 || z == NaN, zs)
            xs, ys
        else
            xs, ys, zs
        end
    end

    for (i, gen) in enumerate(res.generators)
        @series begin
            l = sum(gen.lengths)

            seriestype := :path
            linewidth --> 5
            label --> "gen$(i); l = $(round(l, 3))"

            path = copy(gen.path)
            push!(path, path[1])

            xs, ys, zs = get_xyz(res.data, path)
            if all(z -> z == 0 || z == NaN, zs)
                xs, ys
            else
                xs, ys, zs
            end
        end
    end
end


text:
using JuMP, Cbc, CSV

file = "compatibility_score.csv"

scores_import = CSV.read(file, header=false)
scores = Matrix{Int}(scores_import)
dims = size(scores,1)

dateModel = Model(Cbc.Optimizer)

@variable(dateModel, matches[1:dims, 1:dims], binary=true)
@objective(dateModel, Max, sum(matches[:,:].*scores[:,:]))

for i = 1:dims
    @constraint(dateModel, sum(matches[i,:]) <=1 )
    @constraint(dateModel, sum(matches[:,i]) <=1 )
end

stats = JuMP.optimize!(dateModel)

print("Solution status: ")
println(termination_status(dateModel))
print("Objective value = ")
println(objective_value(dateModel))
matchList = findall(m->m==1,value.(matches))
println("Enter matchList to console to display list of matches")

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
"""
    get_distance(A::Int, B::Int,
                 nodes::Dict{Int,T} ,
                 vertices_to_nodes::Vector{Int}) where T<:Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF}

Auxiliary function - takes two vertices of graph and return the distance between them.
Used to compute straight line distance heuristic for A* algorithm.

**Arguments**

* `A` : start vertex
* `B` : end vertex
* `nodes` : dictionary of .osm nodes ID's and correspoding points coordinates
* `vertices_to_nodes` : dictionary mapping graph vertices to .osm file nodes
"""
function get_distance(A::Int, B::Int,
                     nodes::Dict{Int,T},
                     vertices_to_nodes::Vector{Int}) where T<:Union{OpenStreetMapX.ENU,OpenStreetMapX.ECEF}
    A,B = vertices_to_nodes[A], vertices_to_nodes[B]
    OpenStreetMapX.distance(nodes[A],nodes[B])
end

function extract_a_star_route(parents::Vector{Int},s::Int, u::Int)
    route = Int[]
    index = u
    push!(route,index)
    while index != s
        index = parents[index]
        push!(route, index)
    end
    reverse!(route)
end

"""
    a_star_algorithm(g::AbstractGraph{U},
                    s::Integer,
                    t::Integer,
                    distmx::AbstractMatrix{T}=Graphs.weights(g),
                    heuristic::Function = (u,v) -> zero(T)) where {T, U}

High level function - implementation of A star search algorithm:
(https://en.wikipedia.org/wiki/A*_search_algorithm).
Based on the implementation in Graphs library,
however significantly improved in terms of performance.

**Arguments**

* `g` : graph object
* `S` : start vertex
* `t` : end vertex
* `distmx` : distance matrix
* `heuristic` : search heuristic function; by default returns zero
"""
function a_star_algorithm(g::Graphs.AbstractGraph{U},  # the g
                          s::Integer,           # the start vertex
                          t::Integer,           # the end vertex
                          distmx::AbstractMatrix{T}=Graphs.weights(g),
                          heuristic::Function = (u,v) -> zero(T)) where {T, U}
    nvg = nv(g)
    checkbounds(distmx, Base.OneTo(nvg), Base.OneTo(nvg))
    frontier = DataStructures.PriorityQueue{U,T}()
    # The value should be `heuristic(s, t)` but it does not matter since it will
    # be `dequeue!`d in the first iteration independently of the value.
    frontier[U(s)] = zero(T)
    dists = fill(typemax(T), nvg)
    dists[s] = zero(T)
    parents = zeros(U, nvg)
    colormap = zeros(UInt8, nvg)
    colormap[s] = 1
    @inbounds while !isempty(frontier)
        u = dequeue!(frontier)
        cost_so_far = dists[u]
        u == t && (return OpenStreetMapX.extract_a_star_route(parents,s,u), cost_so_far)
        for v in Graphs.outneighbors(g, u)
            col = colormap[v]
            if col < UInt8(2)
                dist = distmx[u, v]
                colormap[v] = 1
                path_cost = cost_so_far + dist
                if iszero(col)
                    parents[v] = u
                    dists[v] = path_cost
                    enqueue!(frontier, v, path_cost + heuristic(v,t))
                elseif path_cost < dists[v]
                    parents[v] = u
                    dists[v] = path_cost
                    frontier[v] = path_cost + heuristic(v,t)
                end
            end
        end
        colormap[u] = 2
    end
    Vector{U}(), Inf
end

"""
    a_star_algorithm(m::MapData,
                    s::Integer,
                    t::Integer)
A star search algorithm with straight line distance heuristic

**Arguments**

* `m` : MapData object
* `S` : start vertex
* `t` : end vertex
* `distmx` : distance matrix
"""
function a_star_algorithm(m::MapData,
                    s::Integer,
                    t::Integer)
    heuristic(u,v) = OpenStreetMapX.get_distance(u, v, m.nodes, m.n)
	OpenStreetMapX.a_star_algorithm(m.g,s,t,m.w,heuristic)
end


text:
module RobotModels

using Parameters
using Vec, LinearAlgebra
using GraphUtils

using ..Utils
using ..Trajectories
using ..GridPaths

export
    forward_euler_integration,
    sat,

    UnicycleModel,
    unicycle_model,
    UnicycleState,
    UnicycleAction,

    get_state,
    dynamics,
    convert_to_wheel_velocities,

    Controller,
    TrackingController,
    StabilizeController,
    PivotController,
    BiDirectionalController,
    SwitchingController,
    BiDirectionalSwitchingController,

    UnicycleController,
    construct_unicycle_controller,

    get_action,
    simulate

"""
    `forward_euler_integration`
"""
function forward_euler_integration(model,state,u,dt;nsteps=10)
    s = deepcopy(state)
    for i in 1:nsteps
        sdot = dynamics(model,s,u)
        s += sdot*dt/nsteps
    end
    return s
end

""" saturation function """
sat(x,δ) = abs(x) <= δ ? x : sign(x)*δ

abstract type AbstractRobotModel end

"""
    `UnicycleModel`

    Default parameters are based on the Khepera IV Robot in the Webots simulator
"""
@with_kw struct UnicycleModel <: AbstractRobotModel
    track_width ::Float64 = 0.1054
    wheel_radius::Float64 = 0.021
end
const UnicycleState = Vector{R} where {R<:Real}
const UnicycleAction = Vector{R} where {R<:Real}
function get_state(model::UnicycleModel,pt::TrajectoryPoint)
    state = Vector{Float64}([pt.pos.x, pt.pos.y, atan(pt.heading)])
end
function get_state(model::UnicycleModel,traj::T,t::Float64) where {T<:AbstractTrajectory}
    ref = get_trajectory_point_by_time(traj,t)
    get_state(model,ref)
end
function dynamics(model::UnicycleModel,state,cmd)
    x = state[1]
    y = state[2]
    θ = state[3]
    w = cmd[1]
    v = cmd[2]
    sdot = [v*cos(θ),v*sin(θ),w]
end
function convert_to_wheel_velocities(model::UnicycleModel,cmd::UnicycleAction)
    w = cmd[1]
    v = cmd[2]
    L = model.track_width
    r = model.wheel_radius

    v_wheels = [v - w*(L/2), v + w*(L/2)] # velocity at wheels
    w_wheels = v_wheels / r # rotational velocity of wheels
end
################################################################################
################################## Controllers #################################
################################################################################

abstract type Controller end

################################################################################
############################## TrackingController ##############################
################################################################################
"""
    `TrackingController`

    To be employed when the robot needs to track a time-varying reference signal

    From Lee et al., "Tracking Control of
    Unicycle-Modeled Mobile Robots Using a Saturation Feedback Controller"
    [link]{https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=911382}
"""
@with_kw struct TrackingController <: Controller
    dt  ::Float64 = 0.0
    a   ::Float64 = 4.0
    k0  ::Float64 = 4.0
    μ   ::Float64 = 1.0
    γ   ::Float64 = 0.1
    ϵ   ::Float64 = 0.1
    k1  ::Float64 = 4.0
    b   ::Float64 = 4.0
    # helper
    eps ::Float64 = 10e-9
end

################################################################################
############################### PivotController ################################
################################################################################
"""
    `PivotController`

    To be employed when the robot needs to turn in place.
"""
@with_kw struct PivotController <: Controller
    kw::Float64 = 0.5
    kp::Float64 = 2.0
end
################################################################################
############################# StabilizeController ##############################
################################################################################
"""
    `StabilizeController`

    To be employed when the Robot needs to stabilize about a particular point
"""
@with_kw struct StabilizeController <: Controller
    kw::Float64 = 2.0
    kp::Float64 = 4.0
    kϕ::Float64 = 15.0
end

################################################################################
############################# SwitchingController ##############################
################################################################################
"""
    `SwitchingController`

    A composite controller that switches between operating modes based on the
    active trajectory type.
"""
@with_kw struct SwitchingController{T,P,S} <: Controller
    tracker     ::T    = TrackingController()
    pivoter     ::P       = PivotController()
    stabilizer  ::S   = StabilizeController()
end

"""
    `UnicycleController`
"""
struct UnicycleController{T<:AbstractTrajectory,C<:Controller} <: Controller
    model::UnicycleModel
    traj::T
    controller::C
end
function construct_unicycle_controller(model::UnicycleModel,traj::T,controller::C) where {T<:AbstractTrajectory,C<:Controller}
    UnicycleController(model,traj,controller)
end
function construct_unicycle_controller(model::UnicycleModel,path::GridWorldPath,controller::C) where {C <: Controller}
    base_traj = construct_trajectory(path)
    traj = optimize_velocity_profile_traj_only(base_traj)
    UnicycleController(model,traj,controller)
end
function get_action(controller::C,state::UnicycleState,t::Float64) where {C <: UnicycleController}
    try
        cmd = get_action(controller.controller,controller.traj,state,t)
        w_cmd = convert_to_wheel_velocities(controller.model,cmd)
    catch e
        println(e)
        rethrow()
    end
end

"""
    This function follows Lee et al., "Tracking Control of
    Unicycle-Modeled Mobile Robots Using a Saturation Feedback Controller"

    [link]{https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=911382}

    Inputs:
    * `target` - the current reference point [xr, yr, θr] of the tracked trajectory
    * `ff`     - the current feedforward command [vr, wr] of the tracked trajectory
    * `state`  - the current state [x, y, θ] of the robot
"""
function get_action(controller::TrackingController,target::UnicycleState,ff::UnicycleAction,state::UnicycleState,t::Float64)
    # println("get_action(controller::TrackingController,target::UnicycleState,ff::UnicycleAction,state::UnicycleState,t::Float64)")
    a = controller.a   # 0 < a < vmax - sup_{t >= 0} abs(vr)
    k0 = controller.k0 # k0 > 0
    μ = controller.μ   # μ = 0 or 1
    γ = controller.γ   # 0 < γ < 1
    ϵ = controller.ϵ   # 0 < ϵ < 1/(1+γ)
    k1 = controller.k1 # k1 > 0
    b = controller.b   # b > 0
    eps = controller.eps
    # state (global frame)
    x,y,θ = state[1],state[2],wrap_to_pi(state[3])
    # reference (global frame)
    xr,yr,θr = target[1],target[2],wrap_to_pi(target[3])
    wr,vr = ff[1],ff[2]
    # errors in robot frame
    e = [ cos(θ) sin(θ) 0;
         -sin(θ) cos(θ) 0;
              0      0  1]*[xr-x,yr-y,get_angular_offset(θ,θr)]
    xe,ye,θe = e[1],e[2],e[3]
    # dirty trick to avoid problems with waiting
    # if (abs(vr) <= 0.000001) && (abs(wr) <= 0.000001)
    #     return [0.0,0.0]
    # end
    # change of coordinates
    x0 = θe
    x1 = ye
    x2 = -xe
    # u0 = wr - w
    # u1 = v - vr*cos(x0)
    V1 = x1^2 + x2^2
    u1 = -sat(k0*x2,a)
    h = 1 + γ*cos(μ*t)
    h_dot = -μ*γ*sin(μ*t)
    x0_ = x0 + ϵ*h*x1 / (1 + V1^(1/2))
    α = 1 - ϵ*h*x2 / (1 + V1^(1/2))
    β = ϵ*(
        (h_dot*x1+h*wr*x2+h*vr*sin(x0))/(1 + V1^(1/2))
        - (h*x1/( ((1 + V1^(1/2))^2)*V1^(1/2) +  eps))
        * (x1*vr*sin(x0) - sat(k0*x2,a)*x2)
    )
    u0 = -β/α - sat(k1*x0_,b)
    # transform back
    w = -u0 + wr
    v = u1 + vr*cos(x0)
    if isnan(w) || isnan(v)
        println("NaN command encountered",controller, target, ff, state, t, w, v)
    end
    return [w, v]
end
"""
    `get_action(controller::PivotController,target,ff,state,t)`

    A simple control law to track about a time-varying heading while reducing
    the position error when possible. The idea is to track the heading signal
    and drive forward or backward to reduce the position error.

    Inputs:
    * `target` - the current reference point [xr, yr, θr] of the tracked trajectory
    * `ff`     - the current feedforward command [vr, wr] of the tracked trajectory
    * `state`  - the current state [x, y, θ] of the robot
"""
function get_action(controller::PivotController,target::UnicycleState,ff::UnicycleAction,state::UnicycleState,t::Float64)
    # println("get_action(controller::PivotController,target::Vector{Float64},ff::Vector{Float64},state::Vector{Float64},t::Float64)")
    kw = controller.kw
    kp = controller.kp
    # state (global frame)
    x,y,θ = state[1],state[2],state[3]
    # reference (global frame)
    xr,yr,θr = target[1],target[2],target[3]
    wr,vr = ff[1],ff[2]
    # control law
    dp = dot([xr-x,yr-y],[cos(θ),sin(θ)]) # distance to target position projected onto heading vector
    dθ = get_angular_offset(θ,θr) # heading error
    w = wr + kw*dθ
    v = vr + kp*dp
    if isnan(w) || isnan(v)
        println("NaN command encountered",controller, target, ff, state, t, w, v)
    end
    return [w, v]
end
"""
    `get_action(controller::stabilizeController,target,ff,state,t)`

    A simple control law to stabilize about a static pose.

    Inputs:
    * `target` - the current reference point [xr, yr, θr] of the tracked trajectory
    * `ff`     - the current feedforward command [vr, wr] of the tracked trajectory
    * `state`  - the current state [x, y, θ] of the robot
"""
function get_action(controller::StabilizeController,target::UnicycleState,ff::UnicycleAction,state::UnicycleState,t::Float64)
    # println("get_action(controller::StabilizeController,target::Vector{Float64},ff::Vector{Float64},state::Vector{Float64},t::Float64)")
    kw = controller.kw
    kp = controller.kp
    kϕ = controller.kϕ
    # state (global frame)
    x,y,θ = state[1],state[2],state[3]
    # reference (global frame)
    xr,yr,θr = target[1],target[2],target[3]
    wr,vr = ff[1],ff[2]
    # control law
    dx = xr-x
    dy = yr-y
    dp = dot([dx,dy],[cos(θ),sin(θ)]) # distance to target position projected onto heading vector
    dθ = ([cos(θ),sin(θ),0] × [cos(θr),sin(θr),0])[end] # heading error
    ϕ = atan(dy,dx)
    dϕ = ([cos(θ),sin(θ),0] × [cos(ϕ),sin(ϕ),0])[end] # angular displacement from target position
    w = wr + kϕ*dϕ*dp + kw*dθ
    v = vr + kp*dp
    if isnan(w) || isnan(v)
        println("NaN command encountered",controller, target, ff, state, t, w, v)
    end
    return [w, v]
end

function get_action(controller::C,ref::TrajectoryPoint,state::UnicycleState,t::Float64) where {C<:Controller}
    # println("get_action(controller::C,ref::TrajectoryPoint,state::UnicycleState,t::Float64) where {C<: Controller}")
    target = get_state(UnicycleModel(),ref) # [ref.pos.x, ref.pos.y, atan(ref.heading)]
    ff = [ref.yaw_rate, norm(ref.vel)]
    get_action(controller,target,ff,state,t)
end
function get_action(controller::C,traj::T,state::UnicycleState,t::Float64) where {C<:Controller,T<:AbstractTrajectory}
    # println("get_action(controller::C,traj::T,state::UnicycleState,t::Float64) where {C<: Controller, T <: AbstractTrajectory}")
    ref = get_trajectory_point_by_time(traj,t)
    get_action(controller,ref,state,t)
end
const RevTraj = Union{TT,DenseTrajectory{TT}} where {TT<:ReverseTrajectory}
function get_action(controller::C,traj::T,state::UnicycleState,t::Float64) where {C<:Controller,T<:RevTraj}
    ref = get_trajectory_point_by_time(traj,t)
    reversed_state = [state[1], state[2], pi+state[3]]
    action = get_action(controller,ref,reversed_state,t)
    return [action[1], -1.0*action[2]]
end

function get_action(controller::SwitchingController,traj::T,state::UnicycleState,t::Float64) where {T<:Trajectory}
    # println("get_action(controller::SwitchingController,traj::Trajectory,state::UnicycleState,t::Float64)")
    get_action(controller,get_active_segment(traj,t),state,t)
end
const TrackingTraj = Union{TT,DenseTrajectory{TT},ReverseTrajectory{TT},DenseTrajectory{ReverseTrajectory{TT}}} where {TT<:AbstractTrackingTrajectory}
function get_action(controller::SwitchingController,traj::T,state::UnicycleState,t::Float64) where {T<:TrackingTraj}
    # println("get_action(controller::SwitchingController,traj::DenseTrajectory,state::UnicycleState,t::Float64)")
    get_action(controller.tracker,traj,state,t)
end
const WaitTraj = Union{TT,DenseTrajectory{TT},ReverseTrajectory{TT},DenseTrajectory{ReverseTrajectory{TT}}} where {TT<:AbstractStabilizingTrajectory}
function get_action(controller::SwitchingController,traj::T,state::UnicycleState,t::Float64) where {T<:WaitTraj}
    get_action(controller.stabilizer,traj,state,t)
end
const PivotTraj = Union{TT,DenseTrajectory{TT},ReverseTrajectory{TT},DenseTrajectory{ReverseTrajectory{TT}}} where {TT<:AbstractPivotingTrajectory}
function get_action(controller::SwitchingController,traj::T,state::UnicycleState,t::Float64) where {T<:PivotTraj}
    get_action(controller.pivoter,traj,state,t)
end

################################################################################
################################## Simulation ##################################
################################################################################

function simulate(model::UnicycleModel,controller,traj,state,t0,tf,dt)
    states = [state]
    cmds = Vector{Vector{Float64}}()
    t = t0
    while t <= tf
        # get target state in global frame
        u           = get_action(controller,traj,state,t)
        # target_pt   = get_trajectory_point_by_time(traj, t)
        # u           = get_action(controller,target_pt,state,t)
        state     = forward_euler_integration(model,state,u,dt)
        push!(states, state)
        push!(cmds,u)
        t += dt
    end
    return states, cmds
end

function optimize_controller_params(sim_model,traj,initial_state;
        dt=0.1,
        k0_range=1.0:0.25:4.0,
        k1_range=1.0:0.25:4.0
    )
    t0 = get_start_time(traj)
    tf = get_end_time(traj)

    opt_vec = [NaN,NaN]
    opt_cost = Inf
    for k0 in k0_range
        for k1 in k1_range
            controller = TrackingController(k0=k0,k1=k1)
            states, cmds = simulate(sim_model,controller,traj,initial_state,t0,tf,dt)
            err = 0.0
            for (t,state) in zip(collect(t0:dt:tf),states)
                target_pt = get_trajectory_point_by_time(traj, t)
                target = [target_pt.pos.x, target_pt.pos.y, atan(target_pt.heading)]
                err += norm(target[1:2] - state[1:2])^2
            end
            if err < opt_cost
                opt_vec = [k0,k1]
            end
        end
    end
    opt_vec
end

end # end RobotModels


text:
############################################################
# judiJacobian ###############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiJacobian, judiJacobianException, subsample

############################################################

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections
struct judiJacobian{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    srcGeometry::Geometry
    recGeometry::Geometry
    source
    options::Options
    fop::Function              # forward
    fop_T::Union{Function, Nothing}  # transpose
end

mutable struct judiJacobianException <: Exception
    msg :: String
end

############################################################
## Constructor
"""
    judiJacobian(F,q)

Create a linearized modeling operator from the non-linear modeling operator `F` and \\
the source `q`. `F` is a full modeling operator including source/receiver projections.

Examples
========

1) `F` is a modeling operator without source/receiver projections:

    J = judiJacobian(Pr*F*Ps',q)

2) `F` is the combined operator `Pr*F*Ps'`:

    J = judiJacobian(F,q)

"""
function judiJacobian(F::judiPDEfull, source::judiVector; DDT::DataType=Float32, RDT::DataType=DDT)
# JOLI wrapper for nonlinear forward modeling
    compareGeometry(F.srcGeometry, source.geometry) == true || judiJacobianException("Source geometry mismatch")
    (DDT == Float32 && RDT == Float32) || throw(judiJacobianException("Domain and range types not supported"))
    if typeof(F.recGeometry) == GeometryOOC
        m = sum(F.recGeometry.nsamples)
    else
        m = 0
        for j=1:F.info.nsrc m += length(F.recGeometry.xloc[j])*F.recGeometry.nt[j] end
    end
    n = F.info.n
    if F.info.nsrc > 1
        srcnum = 1:F.info.nsrc
    else
        srcnum = 1
    end
    return J = judiJacobian{Float32,Float32}("linearized wave equation", m, n, F.info, F.model, F.srcGeometry, F.recGeometry, source.data, F.options,
                                           v -> time_modeling(F.model, F.srcGeometry, source.data, F.recGeometry, nothing, v, srcnum, 'J', 1, F.options),
                                           w -> time_modeling(F.model, F.srcGeometry, source.data, F.recGeometry, w.data, nothing, srcnum, 'J', -1, F.options)
                                           )
end

############################################################
## overloaded Base functions

# conj(judiJacobian)
conj(A::judiJacobian{DDT,RDT}) where {DDT,RDT} =
    judiJacobian{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.source,A.options,
        A.fop,
        A.fop_T
        )

# transpose(judiJacobian)
transpose(A::judiJacobian{DDT,RDT}) where {DDT,RDT} =
    judiJacobian{DDT,RDT}("adjoint linearized wave equation",A.n,A.m,A.info,A.model,A.srcGeometry,A.recGeometry,A.source,A.options,
        A.fop_T,
        A.fop
        )

# adjoint(judiJacobian)
adjoint(A::judiJacobian{DDT,RDT}) where {DDT,RDT} =
    judiJacobian{DDT,RDT}("adjoint linearized wave equation",A.n,A.m,A.info,A.model,A.srcGeometry,A.recGeometry,A.source,A.options,
        A.fop_T,
        A.fop
        )

############################################################
## overloaded Base *(...judiJacobian...)

# *(judiJacobian,vec)
function *(A::judiJacobian{ADDT,ARDT},v::AbstractVector{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v,1) || throw(judiJacobianException("shape mismatch"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiJacobian,vec):",A.name,typeof(A),vDT]," / "))
    V = A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiJacobian,vec):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(judiJacobian,judiVector)
function *(A::judiJacobian{ADDT,ARDT},v::judiVector{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v,1) || throw(judiJacobianException("shape mismatch"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiJacobian,judiVector):",A.name,typeof(A),vDT]," / "))
    compareGeometry(A.recGeometry,v.geometry) == true || throw(judiJacobianException("Geometry mismatch"))
    V = A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiJacobian,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(num,judiJacobian)
function *(a::Number,A::judiJacobian{ADDT,ARDT}) where {ADDT,ARDT}
    return judiJacobian{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.source,A.options,
                                v1 -> jo_convert(ARDT,a*A.fop(v1),false),
                                v2 -> jo_convert(ADDT,a*A.fop_T(v2),false)
                                )
end

function A_mul_B!(x::judiVector,J::judiJacobian,y::Array)
    z = J*y
    x.data = z.data
end

function Ac_mul_B!(x::Array,J::judiJacobian,y::judiVector)
    x[:] = adjoint(J)*y
end

############################################################
## overloaded Bases +(...judiJacobian...), -(...judiJacobian...)

# +(judiJacobian,num)
function +(A::judiJacobian{ADDT,ARDT},b::Number) where {ADDT,ARDT}
    return judiJacobian{ADDT,ARDT}("("*A.name*"+N)",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.source,A.options,
                            v1 -> A.fop(v1)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
                            v2 -> A.fop_T(v2)+joConstants(A.n,A.m,b;DDT=ADDT,RDT=ARDT)*v2
                            )
end

# -(judiJacobian,num)
function -(A::judiJacobian{ADDT,ARDT},b::Number) where {ADDT,ARDT}
    return judiJacobian{ADDT,ARDT}("("*A.name*"-N)",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.source,A.options,
                            v1 -> A.fop(v1)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
                            v2 -> A.fop_T(v2)-joConstants(A.n,A.m,b;DDT=ADDT,RDT=ARDT)*v2
                            )
end

# -(judiJacobian)
-(A::judiJacobian{DDT,RDT}) where {DDT,RDT} =
    judiJacobian{DDT,RDT}("(-"*A.name*")",A.m,A.n,A.info,A.model,A.srcGeometry,A.recGeometry,A.source,A.options,
                    v1 -> -A.fop(v1),
                    v2 -> -A.fop_T(v2)
                    )

############################################################
## Additional overloaded functions

# Subsample Jacobian
function subsample(J::judiJacobian{ADDT,ARDT}, srcnum) where {ADDT,ARDT}

    srcGeometry = subsample(J.srcGeometry,srcnum)       # Geometry of subsampled data container
    recGeometry = subsample(J.recGeometry,srcnum)

    info = Info(J.info.n, length(srcnum), J.info.nt[srcnum])
    Fsub = judiModeling(info, J.model, srcGeometry, recGeometry; options=J.options)
    qsub = judiVector(srcGeometry, J.source[srcnum])
    return judiJacobian(Fsub, qsub; DDT=ADDT, RDT=ARDT)
end

getindex(J::judiJacobian,a) = subsample(J,a)


text:
@testset "fn" begin
    @testset "fn cons" begin
        f = (Int ⇒ Int)(x -> x + 1)
        @test f(1) === 2
    end
    @testset "fn des" begin

        f = (Int ⇒ Int)(x -> x + 1)
        g = (Int ⇒ String)(x -> "$x")
        test_des(inp) = @match inp begin
            ::(Int ⇒ Int) => true
            _              => false
        end
        @test test_des(f)
        @test !test_des(g)

    end
end


text:
# Julia wrapper for header: /usr/include/openssl/md5.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Int32 MD5_Init (Ptr{MD5_CTX},) libcrypto
@c Int32 MD5_Update (Ptr{MD5_CTX}, Ptr{None}, size_t) libcrypto
@c Int32 MD5_Final (Ptr{Uint8}, Ptr{MD5_CTX}) libcrypto
@c Ptr{Uint8} MD5 (Ptr{Uint8}, size_t, Ptr{Uint8}) libcrypto
@c None MD5_Transform (Ptr{MD5_CTX}, Ptr{Uint8}) libcrypto


text:
import Base.show

export loo_compare

"""
    ModelComparison

A struct containing the results of model comparison.

# Fields

  - `pointwise::KeyedArray`:
  - `estimates::KeyedArray`: A table containing the results of model comparison, with the
    following columns --
      + `cv_est`: The difference in *total* leave-one-out cross validation scores
        between models.
      + `se_loo_diff`: The standard error for each row of `cv_est`.
      + `weight`: A set of Akaike-like weights assigned to each model, which can be used in
        pseudo-Bayesian model averaging.

# Example

```
┌───────┬──────────┬───────┬────────┐
│       │   cv_est │ se_cv │ weight │
├───────┼──────────┼───────┼────────┤
│ m5_1t │     0.00 │  0.00 │   0.67 │
│ m5_3t │    -0.69 │  0.42 │   0.33 │
│ m5_2t │    -6.68 │  4.74 │   0.00 │
└───────┴──────────┴───────┴────────┘
```

See also: [`PsisLoo`](@ref)
"""
struct ModelComparison
    pointwise::KeyedArray
    estimates::KeyedArray
end

"""
    function loo_compare(
        cv_results::PsisLoo...;
        sort_models::Bool=true,
        [, model_names::Tuple{Symbol}]
    ) -> ModelComparison

Construct a model comparison table from several [`PsisLoo`](@ref) objects.

# Arguments

    - `cv_results`: One or more [`PsisLoo`](@ref) objects to be compared. Alternatively,
    a tuple or named tuple of `PsisLoo` objects can be passed. If a named tuple is passed,
    these names will be used to label each model.
    - `model_names`: A vector or tuple of strings or symbols used to identify models. If
    none, models are numbered using the order of the arguments.
    - `sort_models=true`: Sort models by total score (best first).

See also: [`ModelComparison`](@ref), [`PsisLoo`](@ref), [`psis_loo`](@ref)
"""
function loo_compare(
    cv_results::PsisLoo...;
    model_names::Stringlike=[Symbol("model_$i") for i in 1:n_models],
    sort_models::Bool=true,
) where {Stringlike <: Base.AbstractVecOrTuple{<:Union{AbstractString, Symbol}}}

    model_names = Symbol.(model_names)
    model_names = [model for model in model_names]
    nmodels = length(cv_results)
    if nmodels != length(model_names)
        throw(ArgumentError("Number of model names does not match size of `cv_results`."))
    end

    # Extract relevant values from PsisLoo objects

    estimates = [cv_results[i].estimates(:cv_est, :total) for i in 1:nmodels]
    pointwise = [cv_results[i].pointwise for i in 1:nmodels]
    cv_pointwise = [cv_results[i].pointwise(:cv_est) for i in 1:nmodels]

    if sort_models
        ind = sortperm([estimates[i][1] for i in 1:nmodels]; rev=true)
        cv_results = cv_results[ind]
        estimates = estimates[ind]
        pointwise = pointwise[ind]
        cv_pointwise = cv_pointwise[ind]
        model_names = model_names[ind]
    end

    # Compute differences between models
    cv_est = [estimates[i] - estimates[1] for i in 1:nmodels]
    se_cv_est = [
        sqrt(length(cv_pointwise[1]) * var(cv_pointwise[1] - cv_pointwise[i])) for
        i in 1:nmodels
    ]

    table = cv_est
    table = hcat(table, se_cv_est)

    sumval = sum([exp(estimates[i]) for i in 1:nmodels])
    weight = [exp(estimates[i]) / sumval for i in 1:nmodels]
    table = hcat(table, weight)

    table = wrapdims(table; model=model_names, statistic=[:cv_est, :se_cv_est, :weight])

    # Construct pointwise differences

    pointwise = cat(pointwise...; dims=3)
    pointwise = KeyedArray(
        pointwise;
        data=1:size(pointwise, :data),
        statistic=[:cv_est, :naive_est, :overfit, :mcse, :pareto_k],
        model=model_names,
    )

    # Subtract the effective number of params and elpd ests; leave mcse+pareto_k the same
    base_case = pointwise[data=:, statistic=1:3, model=1]
    @inbounds @simd for model_number in axes(pointwise, :model)
        @. pointwise[:, 1:3, model_number] =
            pointwise[:, 1:3, model_number] - base_case
    end

    return ModelComparison(pointwise, table)

end

function loo_compare(cv_results::NamedTuple; kwargs...)
    return loo_compare(values(cv_results); model_names=keys(cv_results), kwargs...)
end

function loo_compare(cv_results::Tuple; kwargs...)
    return loo_compare(cv_results...; kwargs...)
end

function Base.show(io::IO, ::MIME"text/plain", model_comparison::ModelComparison)
    estimates = model_comparison.estimates
    return pretty_table(
        estimates;
        compact_printing=false,
        header=estimates.statistic,
        row_names=estimates.model,
        formatters=ft_printf("%5.2f"),
        alignment=:r,
    )
end


text:
module C_constraints
    include(dirname(@__FILE__)*"\\Utils_cons.jl")
    include(dirname(dirname(@__FILE__))*"\\Moments.jl")
    using .Utils_cons ; const uc = Utils_cons
    using .Moments ; const mom = Moments
    using LinearAlgebra ; const la = LinearAlgebra

    export make_mon_expo_keys,
           make_PSD_con,
           make_ord4_con,
           make_loc_cons_S1,
           make_loc_cons_S2,
           make_loc_cons_S3,
           make_Gᴿ_con

    function make_mon_expo_keys(d,t::Int;noBlock=false)
        if noBlock
            return mom.make_mon_expo(d,t*2)
        else
            # Man this is a bootleg solution
            _,MMᴿexᴿ = mom.get_ℂ_block_diag(d,t,noBlock=false)
            BB = [keys(MMᴿexᴿ)...]
            global boko = vcat([[b...] for b in [MMᴿexᴿ[BB[1]]...]]...)
            for j in 2:size(BB)[1]
                push!(boko,vcat([[b...] for b in [MMᴿexᴿ[BB[j]]...]]...)...)
            end
            return unique(boko)
        end
    end

    """L([xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜᵀ) ⪰ 0"""
    function make_PSD_con(d,t,Lx;noBlock=false)
        MMCoefᴿ,MMexᴿ = mom.get_ℂ_block_diag(d,t;noBlock=noBlock)
        return  Utils_cons.idx2var_dic(Lx,MMexᴿ,MMCoefᴿ)
    end

    """L(xx*⊗yy*) = ρ  ⟺ Lᴿ_ℝxx̄ᵀ⨂yȳᵀ = ℜe(ρ), Lᴿ_𝕀xx̄ᵀ⨂yȳᵀ = ℑm(ρ)  """
    function make_ord4_con(d,Lx)
        γδ_dict = mom.get_γδ_dict(d,(2,2))
        xx̄ᵀ⨂yȳᵀ= mom.make_xx̄ᵀ⨂yȳᵀ(d,γδ_dict)
        (cℝxx̄ᵀ⨂yȳᵀ,eℝxx̄ᵀ⨂yȳᵀ),(c𝕀xx̄ᵀ⨂yȳᵀ,e𝕀xx̄ᵀ⨂yȳᵀ) = mom.get_ℝℂ_coef_expo(d,xx̄ᵀ⨂yȳᵀ,γδ_dict)
        Lᴿ_ℝxx̄ᵀ⨂yȳᵀ = uc.idx2var_a(Lx,eℝxx̄ᵀ⨂yȳᵀ,cℝxx̄ᵀ⨂yȳᵀ)
        Lᴿ_𝕀xx̄ᵀ⨂yȳᵀ = uc.idx2var_a(Lx,e𝕀xx̄ᵀ⨂yȳᵀ,c𝕀xx̄ᵀ⨂yȳᵀ)
        return Lᴿ_ℝxx̄ᵀ⨂yȳᵀ,Lᴿ_𝕀xx̄ᵀ⨂yȳᵀ
    end

    """ Lᴿ(gᴿ⋅[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁ᵀ) ⪰ 0
        gᴿ = √(ρₘₐₓ) - ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ²) , i ∈ [d_1],
        or   √(ρₘₐₓ) - ((yᵣₑ)ⱼ² + (yᵢₘ)ⱼ²) , j ∈ [d_2] """
    make_loc_con(Lx,sqρ,eᵤ,eᵥ,B,C) = sqρ*uc.idx2var_a(Lx,B,C) - (uc.idx2var_a(Lx,B,C,[2*eᵤ]) + uc.idx2var_a(Lx,B,C,[2*eᵥ]))  # L((√ρₘₐₓ - ((xᵣₑ)ᵢ²+(xᵢₘ)ᵢ²))⋅ηₜ₋₁)
    function make_loc_cons_S1(ρ,d,t,Lx;noBlock=false)
        d₁,d₂ = d ; n = sum(2 .*d) ; sqρ   = sqrt(maximum(norm.(ρ)))
        MMCoefᴿ,MMexᴿ = mom.get_ℂ_block_diag(d,t .- 1;noBlock=noBlock)
        loc_con = Dict()
        for b in keys(MMCoefᴿ)
            for k in 1:d₁ # Constraint:  Lᴿ( (√ρₘₐₓ-((xᵣₑ)ᵢ²-(xᵢₘ)ᵢ²))⋅ηₜ₋₁)) ⪰ 0 for k ∈ [d₁]
                loc_con[(b,"x²_$k")] = make_loc_con(Lx,sqρ,uc.eᵢ(n,k),uc.eᵢ(n,k+d₁),MMexᴿ[b],MMCoefᴿ[b])
            end
            for k in 1:d₂ # Constraint:  Lᴿ( (√ρₘₐₓ-((yᵣₑ)ᵢ²+(yᵢₘ)ᵢ²))⋅ηₜ₋₁) ⪰ 0 for k ∈ [d₂]
                loc_con[(b,"y²_$k")] = make_loc_con(Lx,sqρ,uc.eᵢ(n,k+2*d₁),uc.eᵢ(n,k+2*d₁+d₂),MMexᴿ[b],MMCoefᴿ[b])
            end
        end
        return loc_con
    end

    """Lᴿ(gᴿ⋅[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁ᵀ) ⪰ 0
        gᴿ   = √Tr(ρ) - ∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ²),
                ∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ²) - ∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ²)"""
    tmp(Lx,ind,coef,n,s,l) = sum([uc.idx2var_a(Lx,ind,coef,[2* uc.eᵢ(n,k)]) for k in s:l]) #### SOmething is wrong
    function make_loc_cons_S2(ρ,d,t,Lx;noBlock=false)
        d₁,d₂ = d ; n = sum(2 .*d) ; stρ = sqrt(real(tr(ρ)))
        MMCoefᴿ,MMexᴿ = mom.get_ℂ_block_diag(d, t.- 1,noBlock=noBlock)

        loc_con = Dict(); g₂ = Dict()
        for b in keys(MMCoefᴿ)
            e_η = MMexᴿ[b] ; c_η = MMCoefᴿ[b]
            xRterm = tmp(Lx,e_η,c_η,n,1,d₁)    + tmp(Lx,e_η,c_η,n,d₁+1,2*d₁)       #  Lᴿ(∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            yRterm = tmp(Lx,e_η,c_η,n,2*d₁+1,2*d₁+d₂) + tmp(Lx,e_η,c_η,n,2*d₁+d₂+1,2*d₁+2*d₂)  #  Lᴿ(∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            loc_con[b,"x"] = stρ*uc.idx2var_a(Lx,e_η,c_η) - xRterm #  √Tr(ρ) ⋅ Lᴿ(ηₜ₋₁) - Lᴿ(∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            #loc_con[b,"y"] = stρ*uc.idx2var_a(Lx,MMexᴿ[b],MMCoefᴿ[b]) - yRterm #  √Tr(ρ) ⋅ Lᴿ(ηₜ₋₁) - Lᴿ(∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            g₂[b] = xRterm - yRterm #
        end
        return loc_con, g₂
    end

    """Lᴿ(gᴿ⋅[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁[xᵣₑ,xᵢₘ,yᵣₑ,yᵢₘ]ₜ₋₁ᵀ) ⪰ 0
        gᴿ = Tr(ρ) - ∑ᵈᵢ(uₓᵢ² + vₓᵢ²),
                ∑ᵈᵢ(u_yⱼ² + v_yⱼ²)  - 1"""
    function make_loc_cons_S3(ρ,d,t,Lx;noBlock=false)
        d₁,d₂ = d ; n = sum(2 .*d) ; tr_ρ  = real(tr(ρ))
        MMCoefᴿ,MMexᴿ = mom.get_ℂ_block_diag(d, t.- 1,noBlock=noBlock)
        loc_con    = Dict() ; loc_con_eq = Dict()
        for b in keys(MMCoefᴿ)
            e_η = MMexᴿ[b] ; c_η = MMCoefᴿ[b]
            xRterm = tmp(Lx,e_η,c_η,n,1,d₁)           + tmp(Lx,e_η,c_η,n,d₁+1,2*d₁) #  Lᴿ(∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            yRterm = tmp(Lx,e_η,c_η,n,2*d₁+1,2*d₁+d₂) + tmp(Lx,e_η,c_η,n,2*d₁+d₂+1,2*d₁+2*d₂) #  Lᴿ(∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ )
            loc_con[b]    = tr_ρ*uc.idx2var_a(Lx,MMexᴿ[b],MMCoefᴿ[b]) - xRterm #  Tr(ρ)⋅L(ηₜ₋₁) - L(∑ᵈᵢ((xᵣₑ)ᵢ² + (xᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ ) ⪰ 0
            loc_con_eq[b] = uc.idx2var_a(Lx, MMexᴿ[b],MMCoefᴿ[b]) - yRterm #  1⋅L(ηₜ₋₁) -  L(∑ᵈᵢ((yᵣₑ)ᵢ² + (yᵢₘ)ᵢ² ) ⋅ ηₜ₋₁ ) = 0
        end
        return loc_con, loc_con_eq
    end

    """ρ⊗L(η) - L( (xx*⊗yy*) ⊗ η )⪰ 0 ∀ η ∈ blocks of ([x,̄x,y,̄y]ₜ₋₂[x,̄x,y,̄y]*ₜ₋₂)
     L( ρℝ ⊗ η)  - L( Gℝ[k] ⊗ η)⪰ 0 ∀ η ∈ blocks of [uₓ,vₓ,u_y,v_y]ₜ₋₂[uₓ,vₓ,u_y,v_y]^Tₜ₋₂  """
    function make_Gᴿ_con(ρ,d,t,Lx;noBlock=false)
        ρᴿ = [real(ρ) -imag(ρ); imag(ρ) real(ρ)]
        γδ_dict = mom.get_γδ_dict(d,(2,2))
        xx̄ᵀ⨂yȳᵀ = Moments.make_xx̄ᵀ⨂yȳᵀ(d,γδ_dict)
        c_xx̄ᵀ⨂yȳᵀᴿ,e_xx̄ᵀ⨂yȳᵀᴿ = mom.get_coef_expoᴿ(d,xx̄ᵀ⨂yȳᵀ,γδ_dict)
        MMᴿcoef,MMᴿexᴿ = mom.get_ℂ_block_diag(d, t.- 2,noBlock=noBlock)

        spec_sum(𝐴,𝐵) = [a + b  for a in 𝐴 for b in 𝐵]
        spec_mult(𝐶,𝐷) = [c * d  for c in 𝐶 for d in 𝐷]
        LGℝη = Dict()
        for b in keys(MMᴿcoef)
            Ex_temp   = mom.var_kron_C(e_xx̄ᵀ⨂yȳᵀᴿ,MMᴿexᴿ[b],spec_sum)
            Coef_temp = mom.var_kron_C(c_xx̄ᵀ⨂yȳᵀᴿ,MMᴿcoef[b],spec_mult)

            TEMP1 = uc.idx2var_a(Lx,Ex_temp,Coef_temp)  # L(Gᴿ[k]⋅ηₜ₋₂)
            TEMP2 = la.kron(ρᴿ,uc.idx2var_a(Lx,MMᴿexᴿ[b],MMᴿcoef[b])) # L(ρᴿ ⊗ ηₜ₋₂)
            LGℝη[b]   = TEMP2 - TEMP1
        end
        return  LGℝη
    end
end

#
# for b in keys(MMᴿcoef)
#     TEMP1 = sum([tmp2(MMexᴿ[b],MMCoefᴿ[b],k) for k in 1:8]) # ∑ₖL(Gᴿ[k]⋅ηₜ₋₂)
#     TEMP2 = la.kron(ρᴿ,uc.idx2var_arr(Lx,MMexᴿ[b],MMCoefᴿ[b])) # L(ρᴿ ⊗ ηₜ₋₂)
#     LGℝη[b]   = TEMP2 - TEMP1 # L(ρᴿ ⊗ ηₜ₋₂)  - ∑ₖL(Gᴿ[k] ⊗ ηₜ₋₂)
# end

        # LGℝη         = Dict()
        # tmp2(B,C,k) = la.kron(sm[k],ones(D .*size(B))) .*
        #     uc.idx2var_arr(
        #         Lx,
        #         mom.var_kron_C(Gᴿ[k],B),
        #         mom.var_kron_C(fill(0.0,size(Gᴿ[k])...),C))

# make_ord4_con(d,Lx) = uc.idx2varxx̄ᵀtyȳᵀ(Lx,mom.make_xx̄ᵀ⨂yȳᵀ(d))

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
#=
circular_linked_list:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-01
=#

module CircularLinkedList

    # node class
    mutable struct Node
        data
        next
    end

    # head and tail pointer as null
    head = Node(nothing, nothing)
    tail = Node(nothing, nothing)
    head.next = tail
    tail.next = head
    size = 0

    function add(data)
        new_node = Node(data, nothing)

        # check if list is empty
        if head.data == nothing
            global head = new_node
            global tail = new_node
            new_node.next = head
        else
            tail.next = new_node
            tail = new_node
            tail.next = head
        end
        global size += 1
    end

    function delete_middle()
        global head, tail, size
        if head == nothing
            return
        else
            # mid position of list
            if size % 2 == 0
                count = size ÷ 2
            else
                count = (size + 1) ÷ 2
            end
            # check whether head is equal to tail or not
            # if yes then list has only one node
            if head != tail
                tmp = head
                current = nothing
                # current will point to node previous to tmp
                for i in 1:count-1
                    current = tmp
                    tmp = tmp.next
                end
                if current != nothing
                    # tmp is mid node that need to be removed
                    # skip removed node
                    current.next = tmp.next
                    # delete tmp
                    tmp = nothing
                else
                    head = tail = tmp.next
                    tail.next = head
                    tmp = nothing
                end
            else
                head = tail = nothing
            end
        end
        size -= 1
    end

    function show_list()
        current = head
        if head == nothing
            println("List is empty")
            return
        else
            print("$(current.data) ")
            while current.next != head
                current = current.next
                print("$(current.data) ")
            end
            println("")
        end
    end

    function main()
        # build circular linked list
        add(1)
        add(2)
        add(3)
        add(4)

        println("Original Linked List")
        show_list()

        while head != nothing
            delete_middle()
            println("Delete Middle from List")
            show_list()
            println("")
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .CircularLinkedList
    CircularLinkedList.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
###############################################################################
#
# Surface roughness
#
###############################################################################
"""
    aerodynamic_roughness(h_can::FT) where {FT<:AbstractFloat}

Calculate the aerodynamic roughness for wind, given
- `h_can` Average canopy height
"""
function aerodynamic_roughness(h_can::FT) where {FT<:AbstractFloat}
    return h_can * FT(0.07)
end

###############################################################################
#
# Zero plane displacement
#
###############################################################################
"""
    zero_plane_displacement(h_can::FT) where {FT<:AbstractFloat}

Calculate the zero plane displacement height, given
- `h_can` Average canopy height
"""
function zero_plane_displacement(h_can::FT) where {FT<:AbstractFloat}
    return h_can * 2/3
end

###############################################################################
#
# Calculate wind speed at different height
#
###############################################################################
"""
    wind_speed(us::FT, z0::FT, d::FT, z::FT) where {FT<:AbstractFloat}
    wind_speed(wr::FT, zr::FT, z0::FT, d::FT, z::FT) where {FT<:AbstractFloat}

Calculate the wind speed at height `z`, given
- `us` Friction velocity `[m s⁻¹]`
- `d` Zero plane displacement `[m]`
- `z0` Aerodynamic roughnes `[m]`
- `z` Height `[m]`
- `wr` Wind speed at reference height `[m s⁻¹]`
- `zr` Reference height `[m s⁻¹]`

The equation used is
```math
u(z) = \\dfrac{us}{κ} ⋅ \\log \\left( \\dfrac{z - d}{z0} \\right)
```
where `κ` is the Von Kármán constant.
"""
function wind_speed(
            us::FT,
            z0::FT,
            d::FT,
            z::FT
) where {FT<:AbstractFloat}
    return us / K_VON_KARMAN(FT) * log((z-d)/z0)
end

function wind_speed(
            wr::FT,
            zr::FT,
            z0::FT,
            d::FT,
            z::FT
) where {FT<:AbstractFloat}
    return wr * log((z-d)/z0) / log((zr-d)/z0)
end

###############################################################################
#
# Update wind speed at different canopy layer
#
###############################################################################
"""
    wind_speed!(spac::SPACMono{FT}, us::FT) where {FT<:AbstractFloat}
    wind_speed!(spac::SPACMono{FT}, wr::FT, zr::FT) where {FT<:AbstractFloat}

Update wind speed profile in the canopy, given
- `spac` [`SPACMono`](@ref) type struct
- `us` Friction velocity `[m s⁻¹]`
- `wr` Wind speed at reference height `[m s⁻¹]`
- `zr` Reference height `[m s⁻¹]`
"""
function wind_speed!(
            spac::SPACMono{FT},
            us::FT
) where {FT<:AbstractFloat}
    @unpack air_bounds, n_canopy, wind_d, wind_z0, wind_zs = spac;
    spac.winds .= wind_speed.(us, wind_z0, wind_d, wind_zs);

    return nothing
end

function wind_speed!(
            spac::SPACMono{FT},
            wr::FT,
            zr::FT
) where {FT<:AbstractFloat}
    @unpack air_bounds, n_canopy, wind_d, wind_z0, wind_zs = spac;
    spac.winds .= wind_speed.(wr, zr, wind_z0, wind_d, wind_zs);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
__precompile__()

module QMTK
using Compat
using Compat.SparseArrays
using Compat.Random

"""
    __version__

QMTK version number
"""
const __version__ = v"0.1.0"
# NOTE: remember to replace `include` when issue #4600
# is solved

# orthogonal part
include("Consts/Consts.jl")
# include("Utils.jl")
include("Base/Base.jl")

# Depends on Base.jl
include("Lattice/Lattice.jl")
include("Sample/Sample.jl")

# Hamiltonian (depends on Lattice.jl, Basis.jl)
include("Hamiltonian/Hamiltonian.jl")

# Model
include("Model/Model.jl")

# Kronnecker Sugar
include("Kron.jl")

end # module


text:
#module curve_utils

#export snake_normals, eliminate_loopy_stuff, get_outline, closed_curve_to_binary_mat, curve_lengths

using Dierckx
using PolygonOps
using StaticArrays
using IterTools

function normalize!(n)
    lengths = map(i-> sum(n[i,:].^2)^0.5, 1:size(n)[1])
    lindexes = findall(l -> l > 0.0, lengths)
    n[lindexes,:] = n[lindexes,:]./lengths[lindexes]
end

function snake_normals(curve::AbstractArray{T,2}) where T<:AbstractFloat
    X = cat(dims = 1, curve[end,:]', curve, curve[1,:]'); # extended S
    dX = X[1:end-1,:]-X[2:end,:]; # dX
    normals_edges = hcat(dX[:,2], -dX[:,1])
    normalize!(normals_edges)
    normals_vertices = 0.5*(normals_edges[1:end-1,:]+normals_edges[2:end,:])
    normalize!(normals_vertices)
    return normals_vertices
end

function segment_length(a::AbstractArray{T,1},b::AbstractArray{T,1}) where T<:AbstractFloat
    value = sqrt((a[1]-b[1])^2+(a[2]-b[2])^2)
    return value
end

function curve_lengths(arr::AbstractArray{T,2}) where T<:AbstractFloat
    result = Array{Float64,1}(undef,size(arr)[1])
    result[1]=0.0
    for i in 1:size(arr)[1]-1
        result[i+1] = segment_length(arr[i,:],arr[i+1,:])
    end
    csum = cumsum(result)
    return csum
end

function eliminate_loopy_stuff(curve, limit)
    L = size(curve,1)
    tokeep = collect(1:L)
    for i in subsets(1:L, 2)
        #if the two points are not neighbours, and the distance between them is less than the limit and the loop is longer than 4*limit
        if abs(i[1]-i[2]) > 2 && segment_length(curve[i[1],:], curve[i[2],:]) < limit && curve_lengths(curve[i[1]:i[2],:])[end] > 4*limit#limit
            filter!(keeper -> keeper < (i[1]+1) || keeper > (i[2]-1), tokeep)
        end
    end
    curve = curve[tokeep,:]
    return curve
end

function get_xy(cks)
    Re = real(cks)
    Im = imag(cks)
    y = cumsum(Im);
    x = cumsum(Re);
    prepend!(y,[0.0]);
    prepend!(x,[0.0]);
    return cat(x,y, dims =2)
end

function get_ck(xy)
    xy_1 = xy[1:end-1,:]
    xy_2 = xy[2:end,:]
    vecs = xy_2.-xy_1
    return vecs[:,1]+im*vecs[:,2]
end

function get_outline(centerline_points, radius_func)
    L = size(centerline_points,1)
    t = curve_lengths(centerline_points)

    spl = ParametricSpline(t,centerline_points',k=1)
    tspl = range(0, t[end], length=L)

    derr = derivative(spl,tspl)'
    normal = hcat(-derr[:,2], derr[:,1])
    radii = radius_func.(tspl)
    ronsplinetop = spl(tspl)'.+(radii.*normal)
    ronsplinebot = (spl(tspl)'.-(radii.*normal))[end:-1:1,:]

    head = (ronsplinetop[1,:].+ronsplinebot[end,:])./2
    tail = (ronsplinetop[end,:].+ronsplinebot[1,:])./2

    ronsplinetop = cat(head', ronsplinetop, dims = 1)
    ronsplinebot = cat(tail', ronsplinebot, dims = 1)

    outline_xy = cat(ronsplinetop, ronsplinebot, dims=1)

    normal = snake_normals(outline_xy)

    return (outline_xy, normal)
end

#Makes a matrix where the matrix entry is true iff the center of corresponding pixel is not outside the closed curve
#curve is closed curve where first and last point should be the same.
#xs and ys denote the center of each pixel column/rownorm
function closed_curve_to_binary_mat(curve::Array{T1,2}, xs::Array{T2,1}, ys::Array{T2,1}) where {T1,T2 <: Number}
    N = length(curve[:,1])
    poly = map(i -> SVector{2,Float64}(curve[i,1],curve[i,2]), 1:N)
    result = zeros(Int64, length(xs), length(ys))
    i = 0
    for x in xs
        i += 1
        j = 0
        for y in ys
            j +=1
            result[j,i] = inpolygon(SVector(x,y), poly)
        end
    end
    return result
end

#end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
##########################################################
"""
$(TYPEDEF)

Abstract supertype of quantities
"""
abstract type AbstractQuantity{Ti<:Integer} end

##########################################################
"""
$(TYPEDEF)

A continuous quantity is represented by exactly one species number

$(TYPEDFIELDS)
"""
struct ContinuousQuantity{Ti} <: AbstractQuantity{Ti}
    """
    Species number representing the quantity
    """
    ispec::Ti

    """
    Quantity identifier allowing to use the quantity as index
    in parameter fields
    """
    id::Ti
end

"""
     ContinuousQuantity(system,regions; ispec=0, id=0)

Add continuous quantity to the regions listed in `regions`.

Unless specified in `ispec`, the species number is generated
automtically.

Unless specified by `id`, the quantity ID is generated automatically.
"""
function ContinuousQuantity(sys::AbstractSystem{Tv,Ti,Tm},regions; ispec=0, id=0) where {Tv,Ti,Tm}
    if ispec==0
        nspec=num_species(sys)
        nspec=nspec+1
    else
        nspec=ispec
    end
    sys.num_quantities+=1
    if id==0
        id=sys.num_quantities
    end
    enable_species!(sys,nspec,regions)
    ContinuousQuantity{Ti}(nspec,id)
end

##########################################################
"""
$(TYPEDEF)

An interface quantity is represented by exactly one species number

$(TYPEDFIELDS)
"""
struct InterfaceQuantity{Ti}<: AbstractQuantity{Ti}
    """
    Species number representing the quantity
    """
    ispec::Ti

    """
    Quantity identifier allowing to use the quantity as index
    in parameter fields
    """
    id::Ti
end

"""
     InterfaceQuantity(system,regions; ispec=0, id=0)

Add interface quantity to the boundary regions listed in `regions`.

Unless specified in `ispec`, the species number is generated
automtically.

Unless specified by `id`, the quantity ID is generated automatically.
"""
function InterfaceQuantity(sys::AbstractSystem{Tv,Ti,Tm},regions;ispec=1,id=0) where {Tv,Ti,Tm}
    if ispec==0
        nspec=num_species(sys)
        nspec=nspec+1
    else
        nspec=ispec
    end
    enable_boundary_species!(sys,nspec,regions)
    sys.num_quanties+=1
    if id==0
        id=sys.num_quantities
    end
    InterfaceQuantity{Ti}(nspec,id)
end

###########################################################
"""
$(TYPEDEF)

A discontinuous quantity is represented by different species
in neigboring regions.

$(TYPEDFIELDS)
"""
struct DiscontinuousQuantity{Ti} <: AbstractQuantity{Ti}
    """
    Species numbers representing the quantity in each region
    """
    regionspec::Vector{Ti}

    """
    Quantity identifier allowing to use the quantity as index
    in parameter fields
    """
    id::Ti
end

"""
     DiscontinuousQuantity(system,regions; regionspec=nothing, id=0)

Add discontinuous quantity to the regions listed in `regions`.

Unless specified in `regionspec`, the species numbers for each region
are generated automatically.

Unless specified by `id`, the quantity ID is generated automatically.
"""
function DiscontinuousQuantity(sys::AbstractSystem{Tv,Ti,Tm},regions; regionspec=nothing, id=0) where {Tv,Ti,Tm}
    rspec=zeros(Ti,num_cellregions(sys.grid))
    if regionspec==nothing
        nspec=num_species(sys)
        for ireg ∈ regions
            nspec=nspec+1
            enable_species!(sys,nspec,[ireg])
            rspec[ireg]=nspec
        end
    else
        for ireg ∈ regions
            enable_species!(sys,regionspec[ireg],[ireg])
            rspec[ireg]=regionspec[ireg]
        end
    end
    sys.num_quantities+=1
    if id==0
        id=sys.num_quantities
    end
    quantity=DiscontinuousQuantity{Ti}(rspec,id)
    quantity
end

"""
    $(SIGNATURES)

Number of quantities defined for system
"""
num_quantities(system::AbstractSystem) = system.num_quantities

"""
    subgrids(quantity, system)

Return a vector of subgrids containing a subgrid for each
region where discontinuous quantitiy is defined.
"""
function subgrids(quantity::DiscontinuousQuantity, sys)
    grid=sys.grid
    subgrids=Vector{ExtendableGrid}(undef,0)
    for ireg=1:num_cellregions(grid)
        ispec=quantity.regionspec[ireg]
        if ispec>0
            push!(subgrids,subgrid(grid,[ireg]))
        end
    end
    subgrids
end

"""
    views(quantity, subgrids,system)

Return a vector of subgrids containing a subgrid for each
region where discontinuous quantitiy is defined.
"""
function views(U, quantity::DiscontinuousQuantity, subgrids,sys)
    grid=sys.grid
    projections=Vector[]
    j=1
    for ireg=1:num_cellregions(grid)
        ispec=quantity.regionspec[ireg]
        if ispec>0
            push!(projections,view(U[ispec,:],subgrids[j]))
            j=j+1
        end
    end
    projections
end

function views(U, quantity::ContinuousQuantity, subgrids,sys)
    grid=sys.grid
    projections=Vector[]
    j=1
    for ireg=1:num_cellregions(grid)
        push!(projections,view(U[quantity.ispec,:],subgrids[j]))
        j=j+1
    end
    projections
end

# just return the first which comes into mind.
# we need to ensure homogeneity of bregion-region structure.
# if that is true, this works.
function get_cellregion(sys,ibc)
    grid=sys.grid
    bfregions=grid[BFaceRegions]
    cregions=grid[CellRegions]
    for ibface=1:num_bfaces(sys.grid)
        if bfregions[ibface]==ibc
            bfcells=grid[BFaceCells]
            return cregions[bfcells[ibface,1]]
        end
    end
    return 0
end

"""
    boundary_dirichlet(system, quantity, ibc, value)

Set Dirichlet boundary `value` for `quantity` at boundary `ibc`.
"""
function boundary_dirichlet!(sys::AbstractSystem, q::DiscontinuousQuantity, ibc, v)
    boundary_dirichlet!(sys,
                        q.regionspec[get_cellregion(sys,ibc)],
                        ibc,
                        v)
end

function boundary_dirichlet!(sys::AbstractSystem, q::ContinuousQuantity, ibc, v)
    boundary_dirichlet!(sys,
                        q.ispec,
                        ibc,
                        v)
end

"""
    node[quantity]
    edge[quantity]
Return species number on [`AbstractNode`](@ref) or [`AbstractEdge`](@ref)
"""
Base.getindex(q::ContinuousQuantity,node::AbstractNode)=q.ispec
Base.getindex(q::AbstractQuantity,edge::AbstractEdge)=q.ispec

Base.getindex(q::DiscontinuousQuantity,edge::Edge)=@inbounds q.regionspec[edge.region]
Base.getindex(q::DiscontinuousQuantity,edge::BEdge)=nothing
Base.getindex(q::DiscontinuousQuantity,node::Node)=@inbounds q.regionspec[node.region]

"""
    bnode[quantity]
Return species number of discontinuous quantity region `ireg`  adjacent
to  [`BoundaryNode`](@ref) for outer boundary nodes.
"""
Base.getindex(q::DiscontinuousQuantity,bnode::BNode)=@inbounds q.regionspec[bnode.cellregions[1]]

"""
    bnode[quantity,ireg]
Return species number of discontinuous quantity region `ireg`  adjacent
to  [`BoundaryNode`](@ref).
"""
Base.getindex(q::DiscontinuousQuantity,bnode::BNode,j)=@inbounds q.regionspec[bnode.cellregions[j]]

"""
    u[q,j]
Return value of quantity in unknowns on edge in flux callbacks.
"""
Base.getindex(u::AbstractEdgeData,q::AbstractQuantity,j)=@inbounds u[q[u.geom],j]

"""
    u[q]
Return value of quantity in unknowns on node in  node callbacks.
"""
Base.getindex(u::AbstractNodeData,q::AbstractQuantity)=@inbounds u[q[u.geom]]

"""
    f[q]=value
Set rhs value for quantity in callbacks
"""
Base.setindex!(f::AbstractNodeData,v,q::AbstractQuantity)=@inbounds f[q[f.geom]]=v

"""
    u[q,ireg]
Return value of discontinuous quantity in unknowns adjacent to unknowns on boundary node.
"""
Base.getindex(u::BNodeUnknowns,q::DiscontinuousQuantity,j)=@inbounds u[q[u.geom,j]]
Base.getindex(f::BNodeRHS,q::DiscontinuousQuantity,j)=@inbounds f[q[f.geom,j]]

"""
    f[q,ireg]=v
Set rhs value for discontinuous quantity in adjacent regions of  boundary node.
"""
Base.setindex!(f::BNodeRHS,v,q::DiscontinuousQuantity,j)=@inbounds f[q[f.geom,j]]=v

"""
    M[q,i]

Access columns  `M` using id of quantity `q`
"""
Base.getindex(m::AbstractMatrix,q::AbstractQuantity,j)= m[q.id,j]

"""
    M[q,i]

Set element of `M` using id of quantity `q`
"""
Base.setindex!(m::AbstractMatrix,v,q::AbstractQuantity,j)= m[q.id,j]=v

"""
    A[q]

Access columns  of Array `A` using id of quantity `q`
"""
Base.getindex(A::AbstractArray,q::AbstractQuantity)= A[q.id]

"""
    A[q]

Set element of `A` using id of quantity `q`
"""
Base.setindex!(A::AbstractMatrix,v,q::AbstractQuantity)= A[q.id]=v