text:
#############################################################################
#Copyright (c) 2016 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
#############################################################################

function rameau(paths::Tuple,switches::Tuple;input_properties=('\t',0),prediction_coef=[0.0059;0.0005],temperature=23.0,laser=532.0,lb_break=1600.,hb_start=2600.,basetype="gcvspline")

	# some function definition
	calibration_model(x, p) = p[1].*x
	baseline_lfmodel(x, p) = p[1].*x + p[2].*x.^2 + p[3].*x.^3

	if switches[1] == "internal"

		scale = 1

		liste=readcsv(paths[1], skipstart=1)
		names = liste[:,1]
		water = Array{Float64}(liste[:,3])
		smo_lf = Array{Float64}(liste[:,4])
		smo_hf = Array{Float64}(liste[:,5])
		roi = ones(Int((size(liste,2)-5)/2),2,size(liste,1))
		for j = 1:size(liste,1)
			roi[:,1,j] = (liste[j,6:2:end])
			roi[:,2,j] = (liste[j,7:2:end])
		end

		rws = ones(size(liste,1),3)

		for i = 1:size(liste,1)
			spectra = readdlm(string(paths[2],names[i]),input_properties[1],skipstart=input_properties[2])

			if spectra[end,1] < spectra[1,1]
				spectra = flipdim(spectra,1)
			end

			x = spectra[:,1]
			spectra[:,2] = spectra[:,2] - minimum(spectra[:,2])
			y = spectra[:,2]./maximum(spectra[:,2]).*scale + 0.1# area normalisation, and we avoid any 0
			ratio_bkg = minimum(y)/maximum(y) # to keep a record of the ratio of maximum signal intensity over minimum background intensity

			#### PRELIMINARY STEP: FIRST WE GRAB THE GOOD SIGNAL IN THE ROI
			interest_index::Array{Int64} = find(roi[1,1,i] .<= x[:,1] .<= roi[1,2,i])
			if size(roi)[1] > 1
				for j = 2:size(roi)[1]
					interest_index = vcat(interest_index,  find(roi[j,1,i] .<= x[:,1] .<= roi[j,2,i]))
				end
			end
			interest_x = x[interest_index,1]
			interest_y = y[interest_index,1]

			if switches[3] == "yes" # if user asks for the experimental baseline

				if switches[4] == "yes"
					x, y_long, ~ = tlcorrection([x[:] y[:]],temperature,laser) # Long correction
					y_long=y_long./maximum(y_long).*scale

					### THE FOLLOWING IS PURELY EXPERIMENTAL

					lb = roi[:,1,i]
					hb = roi[:,2,i]

					roi_hf =  [lb[hb.>hb_start] hb[hb.>hb_start]]
					roi_lf =  [lb[hb.<lb_break] hb[hb.<lb_break]]

					y_hf = y_long[x.>hb_start]
					y_lf = y_long[x.<lb_break]

					# A first linear baseline in signals between ~1300 and ~ 2000 cm-1
					#y_calc_hf, baseline_hf = baseline(x,y_long,[2600. 3100.;3775. 3900.],"gcvspline",p=smo_hf[i])
					#y_calc_hf, baseline_hf = baseline(x,y_long,[roi_hf[end-1,1] roi_hf[end-1,2];roi_hf[end,1] roi_hf[end,2]],"poly",p=3.0)
					y_calc_hf, baseline_hf = baseline(x,y_long,[3000. 3100;3750. 3800.],"poly",p=1.0)

					x_fit = [x[0 .< x .< 150];x[1250. .< x .< 1500]]
					y_fit = [y_long[0 .< x .< 150];y_long[1250. .< x .< 1500.]]

					#sss = size(x_fit,1)
					#mod = Model(solver=IpoptSolver())
					#@variable(mod, k1 >= 0.0,start=1.0)
					#@variable(mod, k2 >= 0.0, start = 1.0)

					#@NLexpression(mod, base[k=1:sss],k1*x_fit[k] + k2*x_fit[k]*x_fit[k])
					#@NLobjective(mod,Min,sum{(base[k] - y_fit[k])^2,k=sss})
					#status = solve(mod)

					#p1 = getvalue(k1)
					#p2 = getvalue(k2)

					#fit_bas = curve_fit(baseline_lfmodel, , [y_long[0 .< x .< 100];y_long[1250. .< x .< 1500.]], [0.5,0.5,0.5])
					#coef_bas = fit_bas.param
					#baseline_lf = baseline_lfmodel(x,[p1;p2])
					#y_calc_lf = y_long - baseline_lf

					y_calc_lf, baseline_lf = baseline(x,y_long,[0. 150;1250. 1600.],"poly",p=3.0)
					#y_calc_lf, baseline_lf = baseline(x,y_long,roi[:,:,i],basetype,p=smo_lf[i])

					y_calc2 = [y_calc_lf[x.<=1600.];y_calc_hf[x.>1600.]]
					bas2 = [baseline_lf[x.<=1600.];baseline_hf[x.>1600.]]
					#y_calc2 = y_calc2[:,1]./trapz(x[:],y_calc2[:,1]).*scale # area normalisation

					figure(figsize=(20,20))
					suptitle("Baseline sub. with Long corr., sp. $(names[i])")

					subplot(3,2,(1,2))
					plot(x,y,"black",label="Raw sp.")
					plot(x,y_long,"blue",label="Long corr. sp.")
					scatter(interest_x,y_long[interest_index,1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x,bas2,"red",label="Baseline")
					legend(loc="best",fancybox="true")

					subplot(3,2,3)
					plot(x[x.<1500],y_long[x.<1500],"blue",label="Long corr. sp.")
					scatter(interest_x[interest_x.<1500],y_long[interest_index[interest_x.<1500],1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x[x.<1500],bas2[x.<1500],"red",label="Baseline")
					xlim(0,1500)
					legend(loc="best",fancybox="true")

					subplot(3,2,4)
					plot(x[x.>2800],y_long[x.>2800],"blue",label="Long corr. sp.")
					scatter(interest_x[interest_x.>2800],y_long[interest_index[interest_x.>2800],1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x[x.>2800],bas2[x.>2800],"red",label="Baseline")
					xlim(2800,4000)
					legend(loc="best",fancybox="true")

					subplot(3,2,(5,6))
					plot(x,y_calc2,"cyan",label="Final sp.")
					xlabel(L"Raman shift, cm$^{-1}$")
					ylabel("Intensity, area normalised")
					legend(loc="best",fancybox="true")

				else
					# ALSO FULLY EXPERIMENTAL!!!! IT DOES NOT WORK!
					y_hf = y[x.>hb_start]
					y_lf = y[x.<lb_break]

					# A first linear baseline in signals between ~1300 and ~ 2000 cm-1
					y_calc_hf, baseline_hf = baseline(x,y,[2600. 2900.;3775. 3900.],"poly",p=3.0)

				    #A first linear baseline in signals between ~1300 and ~ 2000 cm-1
				   	y_calc1_part1, baseline1 = baseline(x,y,[1300 1800.],"poly",p=1.0)

					# we get a mean signal for background below 1300 cm-1
					X_low_lf = x[y .== minimum(y[100 .<x.<1400])]
					b_low_lf = mean(y[X_low_lf[1]-10 .< x .< X_low_lf[1]+10])

					baseline1[x.<1250.] = b_low_lf
					y_calc1_part1 = y - baseline1

					bas2 = [baseline1[x.<lb_break];baseline_hf[x.>lb_break]]
					y_calc2 = [y_calc1_part1[x.<lb_break];y_calc_hf[x.>lb_break]]

					figure()
					plot(x,y,color="black")
					plot(x,baseline_hf,color="red")
					plot(x,y_calc2,color="blue")
				end

			else # only a single baseline treatment is asked

				# we test if thre Long correction is asked, and react accordingly for baseline substraction
				if switches[4] == "yes"
					x, y_long, ~ = tlcorrection([x[:] y[:]],temperature,laser) # Long correction
					y_long=y_long./maximum(y_long).*scale

					y_calc2, bas2 = baseline(x,y_long,roi[:,:,i],basetype,p=smo_lf[i])
					y_calc2 = y_calc2[:,1]./trapz(x[:],y_calc2[:,1]).*(scale./10) # area normalisation

					figure(figsize=(20,20))
					suptitle("LL2012 method, sp. $(names[i])")

					subplot(3,2,(1,2))
					plot(x,y,"black",label="Raw sp.")
					plot(x,y_long,"blue",label="Long corr. sp.")
					scatter(interest_x,y_long[interest_index,1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x,bas2,"red",label="Baseline")
					legend(loc="best",fancybox="true")

					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,3)
					plot(x[x.<1500],y_long[x.<1500],"blue",label="Long corr. sp.")
					scatter(interest_x[interest_x.<1500],y_long[interest_index[interest_x.<1500],1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x[x.<1500],bas2[x.<1500],"red",label="Baseline")
					xlim(0,1500)
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,4)
					plot(x[x.>2800],y_long[x.>2800],"blue",label="Long corr. sp.")
					scatter(interest_x[interest_x.>2800],y_long[interest_index[interest_x.>2800],1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x[x.>2800],bas2[x.>2800],"red",label="Baseline")
					xlim(2600,4000)
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,(5,6))
					plot(x,y_calc2,"cyan",label="Final sp.")
					xlabel(L"Raman shift, cm$^{-1}$")
					ylabel("Intensity, area normalised")
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

				else
					y_calc2, bas2 = baseline(x,y,roi[:,:,i],basetype,p=smo_lf[i])
					y_calc2 = y_calc2[:,1]./trapz(x[:],y_calc2[:,1]).*scale # area normalisation

					figure()
					plot(x,y.*10,"black",label="Raw sp.")
					scatter(interest_x,interest_y.*10,s=10.,color="red",label="ROI")
					plot(x,bas2.*10,"red",label="Baseline")
					plot(x,y_calc2,"cyan",label="Final sp.")
					xlabel(L"Raman shift, cm$^{-1}$")
					ylabel("Intensity, area normalised")
					title("Baseline sub., NO Long corr.")
					legend(loc="best",fancybox="true")
				end

			end

			# saving the corrected spectrum
			writecsv(string(paths[3],names[i]),[x y_calc2 bas2])

			#Saving the generated figure (specified directory should exist)
			savefig(string(paths[4],"baseline_",names[i],".pdf"))
    		close()
			if switches[3] == "yes"
				# calc=ulating the areas under silicate and water bands

				if switches[4] == "yes"
					As = (trapz(x[150 .<x.<600],y_calc2[150 .<x.<600])+trapz(x[1000 .<x.<1250],y_calc2[1000 .<x.<1250]))

					Aw = trapz(x[3100 .< x .<3750],y_calc2[3100 .< x .<3750])
				else

					#idx1a = find(y[400 .<x.<600] .== maximum(y[400 .<x.<600]))
					#idx1b = find(y[600 .<x.<1000] .== minimum(y[600 .<x.<1000]))

					#idx2a = find( y[900 .<x.<1250] .== maximum(y[900 .<x.<1250]))
					#idx2b = find( y[1000 .<x.<1250] .== minimum(y[1000 .<x.<1250]))

					#idx3a = find( y[3000 .<x.<3700] .== maximum(y[3000 .<x.<3700]))
					#idx3b = find( y[3700 .<x.<3800] .== minimum(y[3700 .<x.<3800]))

					#As = mean(y[(idx1a[1]-4):(idx1a[1]+4),1]) - mean(y[(idx1b[1]-4):(idx1b[1]+4),1]) + mean(y[(idx2a[1]-4):(idx2a[1]+4),1]) - mean(y[(idx2b[1]-4):(idx2b[1]+4),1])
					#Aw = mean(y[idx3a[1]-4 : idx3a[1]+4,1]) - mean(y[idx3b[1]-4 : idx3b[1]+4,1])

					As = maximum(y[400 .<x.<600]) - minimum(y[600 .<x.<1000]) + maximum(y[900 .<x.<1250]) - minimum(y[1000 .<x.<1250])
					Aw = maximum(y[3000 .<x.<3700]) - minimum(y[3700 .<x.<3800])
				end
			else
				# calc=ulating the areas under silicate and water bands
				As = trapz(x[150 .< x .<1300],y_calc2[150 .< x .<1300])
				Aw = trapz(x[2800 .< x .<3750],y_calc2[2800 .< x .<3750])
			end

			# recording them
			rws[i,1] = As;
			rws[i,2] = Aw;
			rws[i,3] = Aw./(As);
			println("Spectrum $(names[i]), the rws is $(rws[i,3])")
		end

		if switches[2] == "yes"
			fit = curve_fit(calibration_model, rws[:,3], water[:]./(100.-water[:]), [0.5])
			coef = fit.param
			sigma = estimate_errors(fit, 0.95)

			rws_calibration = collect(0:0.01:round(maximum(rws[:,3]),2))
			water_ratio_calibration = calibration_model(rws_calibration,coef)
			water_compare =  100.*calibration_model(rws[:,3],coef)./(calibration_model(rws[:,3],coef)+1) # eq. 3 Le Losq et al. (2012)

			rmse_calibration = sqrt(1./(size(rws,1)-1).*sum((water_compare-water).^2))
			rmse_ratio_ws = rmse_calibration/100*(rmse_calibration/(100-rmse_calibration)+1)

			figure()
			scatter(rws[:,3],water[:]./(100.-water[:]))
			plot(rws_calibration,water_ratio_calibration,color="red",linewidth=2.0)
			plot(rws_calibration,water_ratio_calibration+rmse_ratio_ws,color="red",linestyle="--",linewidth=1.0)
			plot(rws_calibration,water_ratio_calibration-rmse_ratio_ws,color="red",linestyle="--",linewidth=1.0)
			plot(rws_calibration,water_ratio_calibration,color="red")
			xlim(0,maximum(rws_calibration)+1./4*maximum(rws_calibration))
			ylim(0,maximum(water_ratio_calibration)+1./4*maximum(water_ratio_calibration))

			xlabel(L"A$_{water}$/A$_{silicates}$, area ratio",fontsize=18,fontname="Arial")
			ylabel("Water/Glass, weight ratio",fontsize=18,fontname="Arial")
			annotate(L"R$_{ws}$=",xy=(0.3,0.9),xycoords="axes fraction",fontsize=18,fontname="Arial",horizontalalignment="center")
			annotate("$(round(coef[1],5)) +/- $(round(sigma[1],5))",xy=(0.3,0.8),xycoords="axes fraction",fontsize=18,fontname="Arial",horizontalalignment="center")
			annotate("Standard deviation\n= $(round(rmse_calibration,2)) wt%",xy=(0.7,0.3),xycoords="axes fraction",horizontalalignment="center",fontsize=18,fontname="Arial")
			savefig(paths[6])

			writecsv(string(paths[5]), ["spectrum" "product" "Rws" "Water input" "Water Raman";liste[:,1] liste[:,2] rws[:,3] water water_compare])
			#return [rws[:,3] water water_compare]
		else
			water_predicted = 100.*(rws[:,3].*prediction_coef[1]./(rws[:,3].*prediction_coef[1]+1))
			water_predicted_high = 100.*(rws[:,3].*sum(prediction_coef)./(rws[:,3].*sum(prediction_coef)+1))
			writecsv(string(paths[5]), ["spectrum" "product" "Rws" "Water Raman" "Error";liste[:,1] liste[:,2] rws[:,3] water_predicted water_predicted_high-water_predicted])
			#return [rws[:,3] water_predicted]
		end

	elseif switches[1] == "external"
		liste=readcsv(paths[1], skipstart=1)

		reference = liste[:,1]
		# reference product names are in first column for record
		reference_density = liste[:,3]
		reference_water = liste[:,4]

		sample = liste[:,5]
		# sample product names are in column 6 for record
		sample_density = liste[:,7]

		water = ones(size(liste,1),1) # the array to record water content

		for i = 1:size(liste,1)

			#importing data
			reference_sp = readdlm(string(paths[2],reference[i]),input_properties[1],skipstart=input_properties[2])
			sample_sp = readdlm(string(paths[2],sample[i]),input_properties[1],skipstart=input_properties[2])

			if reference_sp[end,1] < reference_sp[1,1]
				reference_sp = flipdim(reference_sp,1)
			end

			if sample_sp[end,1] < sample_sp[1,1]
				sample_sp = flipdim(sample_sp,1)
			end

			# temperature - excitation line correction with the Hehlen formula
			x_ref, reference_sp_long, ~ = tlcorrection(reference_sp,temperature,laser,correction="hehlen",normalisation="no",density=reference_density[i])
			x_sp, sample_sp_long, ~ = tlcorrection(sample_sp,temperature,laser,correction="hehlen",normalisation="no",density=sample_density[i])

			# Linear baseline subtraction
			roi_hf = [3000. 3100.; 3800. 3900.] # the region where the basline is constrained
			reference_sp_corr2, reference_baseline = baseline(x_ref,reference_sp_long,roi_hf,"poly",p=1.0)
			sample_sp_corr2, sample_baseline = baseline(x_sp,sample_sp_long,roi_hf,"poly",p=1.0)

			# Area calculation between 3100 and 3800 cm-1
			Area_reference = trapz(x_ref[3100. .<x_ref[:,1].<3800.],reference_sp_corr2[3100. .<x_ref[:,1].<3800.])
			Area_sample = trapz(x_sp[3100. .<x_sp[:,1].<3800.],sample_sp_corr2[3100. .<x_sp[:,1].<3800.])

			# With the intensity
			Intensity_reference = maximum(reference_sp_corr2[3100. .<x_ref[:,1].<3800])
			Intensity_sample = maximum(sample_sp_corr2[3100. .<x_ref[:,1].<3800])

			# Water calculation
			water[i] = (reference_water[i]./1.8*reference_density[i]) .*Area_sample./Area_reference # water is in mol/L
			#water[i] = (reference_water[i]./1.8*reference_density[i]) .*Intensity_sample./Intensity_reference # water is in mol/L

			water[i] = water[i] *1.8./sample_density[i] # water converted in wt%

			figure(figsize=(7.5,15))

			subplot(311)
			title("External calibration: $(round(water[i],2)) wt% water,\n sp. $(sample[i]) with ref. $(reference[i]) ")
			plot(reference_sp[:,1],reference_sp[:,2],color="black",label="Reference")
			plot(sample_sp[:,1],sample_sp[:,2],color="blue",label="Sample")
			xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
			legend(loc="best",frameon=false)

			subplot(312)
			plot(x_ref[:],reference_sp_long[:],color="black",label=L"T-$\nu$ corrected reference")
			plot(x_ref[:],reference_baseline,color="grey",label="Reference baseline")
			plot(x_sp[:],sample_sp_long[:],color="blue",label=L"T-$\nu$ corrected sample")
			plot(x_sp[:],sample_baseline,color="cyan",label="Sample baseline")
			xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
			ylabel("Normalized intensity, a. u.",fontsize=18,fontname="Arial")
			legend(loc="best",frameon=false)

			subplot(313)
			plot(x_ref,reference_sp_corr2,color="black",label="Ref. spectrum")
			plot(x_sp,sample_sp_corr2,color="blue",label="Sample spectrum")
			xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
			legend(loc="best",frameon=false)
			tight_layout()
			# saving the water content
			writecsv(string(paths[5]), ["Spectrum" "Sample Name" "Estimated water content";liste[:,5] liste[:,6] water])

			#Saving the generated figure (specified directory should exist)
			savefig(string(paths[4],"Internal_Std_",sample[i],".pdf"))

		end

	else # or error message
		error("the first switch indicates if you want to use an internal or external calibration mode. Please choose between both.")
	end
end


text:
mutable struct Tree
    ptr::Ptr{TSTree}
    root_node::TSNode
    language::Ptr{TSLanguage}
    Tree() = new()
end

getroot(🌲::Tree) = 🌲.root_node
getlanguage(🌲::Tree) = 🌲.language

ext: jl


text:
using OmegaCore

function test()
  x = 1 ~ StdNormal{Float64}()
  y = 2 ~ StdUniform{Float64}()
  xy = @joint x y
  ω = OmegaCore.Space.LinearΩ([0.3, 0.5], OmegaCore.Tags(), Dict(x => 1, y => 2), [])
  xy(ω)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Autogenerated wrapper script for Modflow6_jll for armv7l-linux-gnueabihf-libgfortran3
export libmf6, mf5to6, mf6, zbud6

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Modflow6")
JLLWrappers.@declare_library_product(libmf6, "libmf6.so")
JLLWrappers.@declare_executable_product(mf5to6)
JLLWrappers.@declare_executable_product(mf6)
JLLWrappers.@declare_executable_product(zbud6)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libmf6,
        "lib/libmf6.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mf5to6,
        "bin/mf5to6",
    )

    JLLWrappers.@init_executable_product(
        mf6,
        "bin/mf6",
    )

    JLLWrappers.@init_executable_product(
        zbud6,
        "bin/zbud6",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ABCDistances
using Base.Test

include("abc.jl")


text:
module MultiScaleTreeGraph
using AbstractTrees
using Printf
using DataFrames
import MutableNamedTuples: MutableNamedTuple
using DelimitedFiles
using OrderedCollections
import XLSX: readxlsx, sheetnames
import SHA: sha1 # for naming the cache variable
import Base.setindex!
import DataFrames: transform!, transform # We define our own version for transforming the MTG
import DataFrames: select!, select # We define our own version for transforming the MTG
import DataFrames: rename!, names # We define our own version for renaming node attributes
import MetaGraphsNext: MetaGraph, code_for, add_edge! # Transform to MetaGraph
import Graphs.DiGraph
using RecipesBase # For plotting with Plots.jl backends

include("types/AbstractNodeMTG.jl")
include("types/Node.jl")
include("read_MTG/read_MTG.jl")
include("read_MTG/strip_comments.jl")
include("read_MTG/utils-string.jl")
include("read_MTG/parse_section.jl")
include("read_MTG/parse_mtg.jl")
include("read_MTG/expand_node.jl")
include("compute_MTG/node_funs.jl")
include("print_MTG/print.jl")
include("compute_MTG/caching.jl")
include("compute_MTG/equality.jl")
include("compute_MTG/indexing.jl")
include("compute_MTG/descendants.jl")
include("compute_MTG/ancestors.jl")
include("compute_MTG/check_filters.jl")
include("compute_MTG/mutation.jl")
include("compute_MTG/append_attributes.jl")
include("compute_MTG/traverse.jl")
include("compute_MTG/transform.jl")
include("compute_MTG/select.jl")
include("compute_MTG/delete_nodes.jl")
include("compute_MTG/filter/filter-funs.jl")
include("write_mtg/update_sections.jl")
include("write_mtg/write_mtg.jl")
include("compute_MTG/insert_nodes.jl")
include("compute_MTG/mutation_helpers.jl")
include("compute_MTG/nleaves.jl")
include("compute_MTG/pipe_model.jl")
include("compute_MTG/get_node.jl")
include("conversion/DataFrame.jl")
include("conversion/MetaGraph.jl")
include("plot_MTG/coordinates.jl")
include("plot_MTG/plot.jl")

export read_mtg
export isleaf
export isroot
export children
export addchild!
export traverse!
export traverse
export transform!, transform, select!, select
export get_root
export nextsibling
export print
export show
export length
export DataFrame
export MetaGraph
export iterate
export siblings
export append!
export @mutate_node!
export @mutate_mtg!
export is_filtered
export delete_nodes!
export delete_node!
export insert_parents!, insert_generations!, insert_children!, insert_siblings!
export insert_parent!, insert_generation!, insert_child!, insert_sibling!
export write_mtg
export is_segment!
export descendants, ancestors, descendants!
export Node
export AbstractNodeMTG
export NodeMTG
export MutableNodeMTG
export (==)
export check_filters
export get_features, get_attributes
export list_nodes
export get_classes
export get_description
export clean_cache!
export branching_order!
export nleaves!, nleaves, nleaves_siblings!
export pipe_model!
export get_node
export new_child_link
export new_id, max_id

end


text:
## High level

# Here is where we can add a default algorithm for computing sensitivities
# Based on problem information!
function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::Nothing,u0,p,args...;kwargs...)
  default_sensealg = if isgpu(u0) && !DiffEqBase.isinplace(prob)
    # only Zygote is GPU compatible and fast
    # so if in-place, try Zygote
    if p === nothing || p === DiffEqBase.NullParameters()
        # QuadratureAdjoint skips all p calculations until the end
        # So it's the fastest when there are no parameters
        QuadratureAdjoint(autojacvec=ZygoteVJP())
      else
        InterpolatingAdjoint(autojacvec=ZygoteVJP())
      end
    else
      if p === nothing || p === DiffEqBase.NullParameters()
        QuadratureAdjoint()
      else
        InterpolatingAdjoint()
      end
    end
  DiffEqBase._concrete_solve_adjoint(prob,alg,default_sensealg,u0,p,args...;kwargs...)
end

function DiffEqBase._concrete_solve_adjoint(prob::SteadyStateProblem,alg,sensealg::Nothing,u0,p,args...;kwargs...)
  default_sensealg = SteadyStateAdjoint()
  DiffEqBase._concrete_solve_adjoint(prob,alg,default_sensealg,u0,p,args...;kwargs...)
end

function DiffEqBase._concrete_solve_adjoint(prob,alg,
                                 sensealg::AbstractAdjointSensitivityAlgorithm,
                                 u0,p,args...;save_start=true,save_end=true,
                                 saveat = eltype(prob.tspan)[],
                                 save_idxs = nothing,
                                 kwargs...)
  _prob = remake(prob,u0=u0,p=p)

  # Force `save_start` and `save_end` in the forward pass This forces the
  # solver to do the backsolve all the way back to `u0` Since the start aliases
  # `_prob.u0`, this doesn't actually use more memory But it cleans up the
  # implementation and makes `save_start` and `save_end` arg safe.
  kwargs_fwd = NamedTuple{Base.diff_names(Base._nt_names(
  values(kwargs)), (:callback_adj,))}(values(kwargs))
  kwargs_adj = NamedTuple{Base.diff_names(Base._nt_names(values(kwargs)), (:callback_adj,:callback))}(values(kwargs))
  if haskey(kwargs, :callback_adj)
    kwargs_adj = merge(kwargs_adj, NamedTuple{(:callback,)}( [get(kwargs, :callback_adj, nothing)] ))
  end

  if typeof(sensealg) <: BacksolveAdjoint
    sol = solve(_prob,alg,args...;save_noise=true,save_start=save_start,save_end=save_end,saveat=saveat,kwargs...)
  elseif ischeckpointing(sensealg)
    sol = solve(_prob,alg,args...;save_noise=true,save_start=true,save_end=true,saveat=saveat,kwargs...)
  else
    sol = solve(_prob,alg,args...;save_noise=true,save_start=true,save_end=true,kwargs...)
  end

  if typeof(sensealg) <: BacksolveAdjoint
    # Saving behavior unchanged
    ts = sol.t
    only_end = length(ts) == 1 && ts[1] == _prob.tspan[2]
    out = DiffEqBase.sensitivity_solution(sol,sol.u,ts)
  elseif saveat isa Number
    if _prob.tspan[2] > _prob.tspan[1]
      ts = _prob.tspan[1]:convert(typeof(_prob.tspan[2]),abs(saveat)):_prob.tspan[2]
    else
      ts = _prob.tspan[2]:convert(typeof(_prob.tspan[2]),abs(saveat)):_prob.tspan[1]
    end
    _out = sol(ts)
    out = if save_idxs === nothing
      out = DiffEqBase.sensitivity_solution(sol,_out.u,sol.t)
    else
      out = DiffEqBase.sensitivity_solution(sol,[_out[i][save_idxs] for i in 1:length(_out)],ts)
    end
    only_end = length(ts) == 1 && ts[1] == _prob.tspan[2]
  elseif isempty(saveat)
    no_start = !save_start
    no_end = !save_end
    sol_idxs = 1:length(sol)
    no_start && (sol_idxs = sol_idxs[2:end])
    no_end && (sol_idxs = sol_idxs[1:end-1])
    only_end = length(sol_idxs) <= 1
    _u = sol.u[sol_idxs]
    u = save_idxs === nothing ? _u : [x[save_idxs] for x in _u]
    ts = sol.t[sol_idxs]
    out = DiffEqBase.sensitivity_solution(sol,u,ts)
  else
    _saveat = saveat isa Array ? sort(saveat) : saveat # for minibatching
    _saveat = eltype(_saveat) <: typeof(prob.tspan[2]) ? convert.(typeof(_prob.tspan[2]),_saveat) : _saveat
    ts = _saveat
    _out = sol(ts)
    out = if save_idxs === nothing
      out = DiffEqBase.sensitivity_solution(sol,_out.u,ts)
    else
      out = DiffEqBase.sensitivity_solution(sol,[_out[i][save_idxs] for i in 1:length(_out)],ts)
    end
    only_end = length(ts) == 1 && ts[1] == _prob.tspan[2]
  end

  _save_idxs = save_idxs === nothing ? Colon() : save_idxs

  function adjoint_sensitivity_backpass(Δ)
    function df(_out, u, p, t, i)
      if only_end
        if typeof(Δ) <: AbstractArray{<:AbstractArray} && length(Δ) == 1 && i == 1
          # user did sol[end] on only_end
          if typeof(_save_idxs) <: Number
            _out[_save_idxs] = -vec(Δ[1])[_save_idxs]
          elseif _save_idxs isa Colon
            vec(_out) .= -vec(Δ[1])
          else
            vec(@view(_out[_save_idxs])) .= -vec(Δ[1])[_save_idxs]
        end
        else
          if typeof(_save_idxs) <: Number
            _out[_save_idxs] = -vec(Δ)[_save_idxs]
          elseif _save_idxs isa Colon
            vec(_out) .= -vec(Δ)
          else
            vec(@view(_out[_save_idxs])) .= -vec(Δ)[_save_idxs]
          end
        end
      else
        if typeof(Δ) <: AbstractArray{<:AbstractArray} || typeof(Δ) <: DESolution
          if typeof(_save_idxs) <: Number
            _out[_save_idxs] = -Δ[i][_save_idxs]
          elseif _save_idxs isa Colon
            vec(_out) .= -vec(Δ[i])
          else
            vec(@view(_out[_save_idxs])) .= -vec(Δ[i][_save_idxs])
          end
        else
          if typeof(_save_idxs) <: Number
            _out[_save_idxs] = -adapt(typeof(u0),reshape(Δ, prod(size(Δ)[1:end-1]), size(Δ)[end])[_save_idxs, i])
          elseif _save_idxs isa Colon
            vec(_out) .= -vec(adapt(typeof(u0),reshape(Δ, prod(size(Δ)[1:end-1]), size(Δ)[end])[:, i]))
          else
            vec(@view(_out[_save_idxs])) .= -vec(adapt(typeof(u0),reshape(Δ, prod(size(Δ)[1:end-1]), size(Δ)[end])[_save_idxs, i]))
          end
        end
      end
    end

    du0, dp = adjoint_sensitivities(sol,alg,args...,df,ts; sensealg=sensealg,
                                    kwargs_adj...)

    du0 = reshape(du0,size(u0))
    dp = p === nothing || p === DiffEqBase.NullParameters() ? nothing : reshape(dp',size(p))

    (nothing,nothing,du0,dp,nothing,ntuple(_->nothing, length(args))...)
  end
  out, adjoint_sensitivity_backpass
end

# Prefer this route since it works better with callback AD
function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::AbstractForwardSensitivityAlgorithm,
                                 u0,p,args...;save_idxs = nothing,
                                 kwargs...)
   _save_idxs = save_idxs === nothing ? (1:length(u0)) : save_idxs
   _prob = ODEForwardSensitivityProblem(prob.f,u0,prob.tspan,p,sensealg)
   sol = solve(_prob,alg,args...;kwargs...)
   _,du = extract_local_sensitivities(sol, sensealg, Val(true))
   out = DiffEqBase.sensitivity_solution(sol,[sol[i][_save_idxs] for i in 1:length(sol)],sol.t)
   function forward_sensitivity_backpass(Δ)
     adj = sum(eachindex(du)) do i
       J = du[i]
       if Δ isa AbstractVector
         v = Δ[i]
       else
         v = @view Δ[:, i]
       end
       J'v
     end
     (nothing,nothing,nothing,adj,nothing,ntuple(_->nothing, length(args))...)
   end
   out,forward_sensitivity_backpass
end

function DiffEqBase._concrete_solve_forward(prob,alg,
                                 sensealg::AbstractForwardSensitivityAlgorithm,
                                 u0,p,args...;save_idxs = nothing,
                                 kwargs...)
   _prob = ODEForwardSensitivityProblem(prob.f,u0,prob.tspan,p,sensealg)
   sol = solve(_prob,args...;kwargs...)
   u,du = extract_local_sensitivities(sol,Val(true))
   _save_idxs = save_idxs === nothing ? (1:length(u0)) : save_idxs
   out = DiffEqBase.sensitivity_solution(sol,[ForwardDiff.value.(sol[i][_save_idxs]) for i in 1:length(sol)],sol.t)
   function _concrete_solve_pushforward(Δself, ::Nothing, ::Nothing, x3, Δp, args...)
     x3 !== nothing && error("Pushforward currently requires no u0 derivatives")
     du * Δp
   end
   out,_concrete_solve_pushforward
end

# Generic Fallback for ForwardDiff
function DiffEqBase._concrete_solve_adjoint(prob,alg,
                                 sensealg::ForwardDiffSensitivity,
                                 u0,p,args...;saveat=eltype(prob.tspan)[],
                                 save_idxs = nothing,
                                 kwargs...)
  _save_idxs = save_idxs === nothing ? (1:length(u0)) : save_idxs
  MyTag = typeof(prob.f)
  pdual = seed_duals(p,MyTag)
  u0dual = convert.(eltype(pdual),u0)

  if (convert_tspan(sensealg) === nothing && (
        (haskey(kwargs,:callback) && has_continuous_callback(kwargs[:callback])) ||
        (haskey(prob.kwargs,:callback) && has_continuous_callback(prob.kwargs[:callback]))
        )) || (convert_tspan(sensealg) !== nothing && convert_tspan(sensealg))

    tspandual = convert.(eltype(pdual),prob.tspan)
  else
    tspandual = prob.tspan
  end

  _prob = remake(prob,u0=u0dual,p=pdual,tspan=tspandual)

  if saveat isa Number
    _saveat = prob.tspan[1]:saveat:prob.tspan[2]
  else
    _saveat = saveat
  end

  sol = solve(_prob,alg,args...;saveat=_saveat,kwargs...)
  _,du = extract_local_sensitivities(sol, sensealg, Val(true))
  out = DiffEqBase.sensitivity_solution(sol,[ForwardDiff.value.(sol[i][_save_idxs]) for i in 1:length(sol)],ForwardDiff.value.(sol.t))
  function forward_sensitivity_backpass(Δ)
    adj = sum(eachindex(du)) do i
      J = du[i]
      if Δ isa AbstractVector
        v = Δ[i]
      else
        v = @view Δ[:, i]
      end
      ForwardDiff.value.(J'v)
    end
    (nothing,nothing,nothing,adj,nothing,ntuple(_->nothing, length(args))...)
  end
  out,forward_sensitivity_backpass
end

function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::ZygoteAdjoint,
                                 u0,p,args...;kwargs...)
    Zygote.pullback((u0,p)->solve(prob,alg,args...;u0=u0,p=p,kwargs...),u0,p)
end

function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::TrackerAdjoint,
                                            u0,p,args...;
                                            kwargs...)

  local sol
  function tracker_adjoint_forwardpass(_u0,_p)
    if DiffEqBase.isinplace(prob)
      # use Array{TrackedReal} for mutation to work
      # Recurse to all Array{TrackedArray}
      _prob = remake(prob,u0=map(identity,_u0),p=_p)
    else
      # use TrackedArray for efficiency of the tape
      _f(args...) = Tracker.collect(prob.f(args...))
      if prob isa SDEProblem
        _g(args...) = Tracker.collect(prob.g(args...))
        _prob = remake(prob,f=DiffEqBase.parameterless_type(prob.f)(_f,_g),u0=_u0,p=_p)
      else
        _prob = remake(prob,f=DiffEqBase.parameterless_type(prob.f)(_f),u0=_u0,p=_p)
      end
    end
    sol = solve(_prob,alg,args...;sensealg=DiffEqBase.SensitivityADPassThrough(),kwargs...)

    if typeof(sol.u[1]) <: Array
      return Array(sol)
    else
      tmp = vec(sol.u[1])
      for i in 2:length(sol.u)
        tmp = hcat(tmp,vec(sol.u[i]))
      end
      return reshape(tmp,size(sol.u[1])...,length(sol.u))
    end
    #adapt(typeof(u0),arr)
    sol
  end

  out,pullback = Tracker.forward(tracker_adjoint_forwardpass,u0,p)
  function tracker_adjoint_backpass(ybar)
    u0bar, pbar = pullback(ybar)
    _u0bar = u0bar isa Tracker.TrackedArray ? Tracker.data(u0bar) : Tracker.data.(u0bar)
    (nothing,nothing,_u0bar,Tracker.data(pbar),nothing,ntuple(_->nothing, length(args))...)
  end

  u = u0 isa Tracker.TrackedArray ? Tracker.data.(sol.u) : Tracker.data.(Tracker.data.(sol.u))
  DiffEqBase.sensitivity_solution(sol,u,Tracker.data.(sol.t)),tracker_adjoint_backpass
end

function DiffEqBase._concrete_solve_adjoint(prob,alg,sensealg::ReverseDiffAdjoint,
                                            u0,p,args...;kwargs...)

  t = eltype(prob.tspan)[]
  u = typeof(u0)[]

  local sol

  function reversediff_adjoint_forwardpass(_u0,_p)
    if DiffEqBase.isinplace(prob)
      # use Array{TrackedReal} for mutation to work
      # Recurse to all Array{TrackedArray}
      _prob = remake(prob,u0=map(identity,_u0),p=_p)
    else
      # use TrackedArray for efficiency of the tape
      _f(args...) = reduce(vcat,prob.f(args...))
      if prob isa SDEProblem
        _g(args...) = reduce(vcat,prob.g(args...))
        _prob = remake(prob,f=DiffEqBase.parameterless_type(prob.f)(_f,_g),u0=_u0,p=_p)
      else
        _prob = remake(prob,f=DiffEqBase.parameterless_type(prob.f)(_f),u0=_u0,p=_p)
      end
    end

    sol = solve(_prob,alg,args...;sensealg=DiffEqBase.SensitivityADPassThrough(),kwargs...)
    t = sol.t
    if DiffEqBase.isinplace(prob)
      u = map.(ReverseDiff.value,sol.u)
    else
      u = map(ReverseDiff.value,sol.u)
    end

    Array(sol)
  end

  tape = ReverseDiff.GradientTape(reversediff_adjoint_forwardpass,(u0, p))
  tu, tp = ReverseDiff.input_hook(tape)
  output = ReverseDiff.output_hook(tape)
  ReverseDiff.value!(tu, u0)
  typeof(p) <: DiffEqBase.NullParameters || ReverseDiff.value!(tp, p)
  ReverseDiff.forward_pass!(tape)
  function reversediff_adjoint_backpass(ybar)
    ReverseDiff.increment_deriv!(output, ybar)
    ReverseDiff.reverse_pass!(tape)
    (nothing,nothing,ReverseDiff.deriv(tu),ReverseDiff.deriv(tp),nothing,ntuple(_->nothing, length(args))...)
  end
  DiffEqBase.sensitivity_solution(sol,u,t),reversediff_adjoint_backpass
end

function DiffEqBase._concrete_solve_adjoint(prob::SteadyStateProblem,alg,sensealg::SteadyStateAdjoint,
                                 u0,p,args...;save_idxs = nothing, kwargs...)

    _prob = remake(prob,u0=u0,p=p)
    sol = solve(_prob,alg,args...;kwargs...)
    _save_idxs = save_idxs === nothing ? Colon() : save_idxs

    if save_idxs === nothing
      out = sol
    else
      out = DiffEqBase.sensitivity_solution(sol,sol[_save_idxs])
    end

    function steadystatebackpass(Δ)
      # Δ = dg/dx or diffcache.dg_val
      # del g/del p = 0
      dp = adjoint_sensitivities(sol,alg;sensealg=sensealg,g=nothing,dg=Δ,save_idxs=save_idxs)
      dp
      (nothing,nothing,nothing,dp,nothing,ntuple(_->nothing, length(args))...)
    end
    out, steadystatebackpass
end


text:
# - https://thenewphalls.wordpress.com/2014/02/19/understanding-object-oriented-programming-in-julia-part-1/
# - http://stackoverflow.com/questions/26800811/julia-constructor-embedding-a-function-inside-a-type
type Visualizer
  ### variables
  pomdp::BotWorld
  CELL_SIZE::Int64
  fig
  ax
  placedTargets::Vector{Tuple{Int64, Int64}}

  #### functions
  init::Function
  #draw::Function
  #update::Function

  function Visualizer(pomdp::BotWorld; CELL_SIZE::Int64=1)
    self = new()
    self.pomdp = pomdp
    self.CELL_SIZE = CELL_SIZE

    self.init = function()

      # plot setup
      self.fig = figure()
      PyPlot.hold(true)
      PyPlot.ion()

      # general setup
      self.placedTargets = Vector{Tuple{Int64, Int64}}()

      # draw initial grid
      #state = getInitialBotState(self.pomdp)
      #self.draw()

    end

    # return self
    return self
  end
end

# - http://stackoverflow.com/questions/30351546/using-matplotlibs-patches-in-julia
# - http://matthiaseisen.com/pp/patterns/p0203/
# - http://matplotlib.org/api/patches_api.html
# - http://matplotlib.org/examples/shapes_and_collections/artist_reference.html
# - https://groups.google.com/forum/#!msg/julia-users/94P2OGaHVUA/oW2thTo4D1cJ
# - http://stackoverflow.com/questions/4098131/how-to-update-a-plot-in-matplotlib
# - http://stackoverflow.com/questions/12998430/remove-xticks-in-a-matplot-lib-plot
function visualizerDraw(self::Visualizer, prevState::BotState=getInitialBotState(self.pomdp), action::BotAction=POMDPs.create_action(self.pomdp), newState::BotState=getInitialBotState(self.pomdp))

  # clear
  PyPlot.clf()
  self.ax = self.fig[:add_subplot](1,1,1)

  # draw grid
  numCellsX = pomdp.WORLD_SIZE_X
  numCellsY = pomdp.WORLD_SIZE_Y
  PyPlot.xlim([-0.5,numCellsX*self.CELL_SIZE+0.5])
  PyPlot.ylim([-0.5,numCellsY*self.CELL_SIZE+0.5])
  numVerticalLines = numCellsX + 1
  numHorizontalLines = numCellsY + 1
  for i=1:numVerticalLines
    PyPlot.plot([self.CELL_SIZE*(i-1), self.CELL_SIZE*(i-1)], [0, self.CELL_SIZE*numCellsY], color="grey", alpha=0.2)
  end
  for i=1:numHorizontalLines
    PyPlot.plot([0, self.CELL_SIZE*numCellsX], [self.CELL_SIZE*(i-1), self.CELL_SIZE*(i-1)], color="grey", alpha=0.2)
  end

  # draw base
  #base = patch.RegularPolygon([pomdp.BASE_LOCATION[1]-0.5, pomdp.BASE_LOCATION[2]-0.5], 4, 0.4)
  base = patch.Rectangle([self.pomdp.BASE_LOCATION[1]-0.9, self.pomdp.BASE_LOCATION[2]-0.9], 0.8, 0.8, alpha=0.5, facecolor="blue")
  self.ax[:add_artist](base)

  # check if placed target this round
  if (prevState.isHoldingTarget == true) && (newState.isHoldingTarget == false) && (newState.numTargetsRemaining == prevState.numTargetsRemaining - 1)
    for targetSpecification in pomdp.TARGET_SPECIFICATIONS
      if (prevState.x == targetSpecification.x) && (prevState.y == targetSpecification.y)
        push!(self.placedTargets, (targetSpecification.x, targetSpecification.y))
      end
    end
  end

  # draw target positions
  for targetSpecification in pomdp.TARGET_SPECIFICATIONS
    # if not already placed
    if !( (targetSpecification.x, targetSpecification.y) in self.placedTargets )
      color = "green"
      fill = false
      linestyle = "dashed"
    # if already placed
    else
      color = "green"
      fill = true
      linestyle = "solid"
    end
    target = patch.RegularPolygon([targetSpecification.x-0.5, targetSpecification.y-0.5], 4, 0.3, alpha=0.8, facecolor=color, fill=fill, linestyle=linestyle)
    self.ax[:add_artist](target)
  end

  # draw obstacles
  for obstacleSpecification in pomdp.OBSTACLE_SPECIFICATIONS
    obstacle = patch.Circle([obstacleSpecification.x-0.5, obstacleSpecification.y-0.5], 0.3, alpha=0.8, facecolor="red", fill=true)
    self.ax[:add_artist](obstacle)
  end

  # draw robot
  if !isEqualStates(newState, getIntendedNextState(pomdp, prevState, action)) # something went unexpected
    color = "red"
  elseif newState.isHoldingTarget == true
    color = "green"
  else
    color = "black"
  end
  if newState.facingDirection == :up
    orientation = 0
  elseif newState.facingDirection == :left
    orientation = 90*pi/180
  elseif newState.facingDirection == :down
    orientation = 180*pi/180
  elseif newState.facingDirection == :right
    orientation = 270*pi/180
  end
  robot = patch.RegularPolygon([newState.x-0.5, newState.y-0.5], 3, 0.25, orientation=orientation, facecolor=color)
  self.ax[:add_artist](robot)

  # show
  #PyPlot.show()
  PyPlot.draw()

end

max_stars_repo_name: abrahambotros/MarsHabBot

max_issues_repo_name: abrahambotros/MarsHabBot

max_forks_repo_name: abrahambotros/MarsHabBot


text:
## Conversion to and from Matrix objects ##

function Graph(m::Matrix{T}; kwargs...) where {T<:Number}
    n, n_ = size(m)
    @assert(n == n_, "adjacency matrix should be square")
    g = Graph(n; kwargs...)
    for i = 1:n, j = 1:n
        if m[i, j] > zero(T)
            addedge!(g, i, j)
        end
    end
    return g
end

function Matrix(g::Graph)
    n = nodecount(g)
    m = zeros(n, n)
    for i = nodes(g), j = outneighbors(g, i)
        m[i, j] += 1
    end
    return m
end


text:
module TestTreeDemo3x7

using Test
import Pkg: TOML
import GraphIdx

include("data.jl")

const fname = data_dir("test", "tree_demo3x7.toml")

@testset "Demo 3x7                       " begin
    tree = TOML.parsefile(fname)["tree"][1]
    pi = tree["parent"] .+ 1
    cidx = GraphIdx.Tree.ChildrenIndex(pi)
    @test GraphIdx.Tree.hierarchy_string(cidx) == tree["str"]
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
struct Rectangle
    top::Int
    left::Int
    bottom::Int
    right::Int
    # return two upper-left and lower-right points of the Rectangle
    Rectangle(p::Position, s::Size) = new(p.y + s.height, p.x, p.y, p.x + s.width)
end

# check if the two rectangles overlap
function overlap(A::Rectangle, B::Rectangle)
    return A.left < B.right && A.right > B.left && A.top > B.bottom && A.bottom < B.top
end

function collide(A::Thing, B::Thing)
    println("Checking collision of thing vs. thing")
    rectA = Rectangle(position(A), size(A))
    rectB = Rectangle(position(B), size(B))
    return overlap(rectA, rectB)
end

function collide(A::Spaceship, B::Spaceship)
    println("Checking collision of spaceship vs. spaceship")
    return true
end

function collide(A::Asteroid, B::Thing)
    println("Checking collision of asteroid vs. thing")
    return true
end

function collide(A::Thing, B::Asteroid)
    println("Checking collision of thing vs. asteroid")
    return true
end

function collide(A::Asteroid, B::Asteroid)
    println("Checking collision of asteroid vs. asteroid")
    return true
end

# ramdomly pick two things and check if they collide
function check_randomly(things)
    for i in 1:5
        two = rand(things, 2)
        collide(two...)
    end
end


text:
"""
timestats(cube, )

Compute the multi temporal statistics of a cube with temporal axis and return a cube with an `Stats` axis.
"""
function timestats(cube;kwargs...)

    indims = InDims("Time")
    funcs = ["Mean", "5th Quantile", "25th Quantile", "Median", "75th Quantile", "95th Quantile",
            "Standard Deviation", "Minimum", "Maximum",
            "Skewness", "Kurtosis", "Median Absolute Deviation"]

    stataxis = CategoricalAxis("Stats", funcs)
    od = OutDims(stataxis)
    stats = mapCube(ctimestats!, cube, indims=indims, outdims=od, kwargs...)
end

function ctimestats!(xout, xin)
    ts = collect(skipmissing(xin))
    m = mean(ts)
    T = eltype(m)
    stats = Vector{T}(undef,12)
    stats[1] = m
    stats[2:6] = quantile(ts, [0.05,0.25,0.5, 0.75,0.95])
    stats[7] = std(ts)
    stats[8] = minimum(ts)
    stats[9] = maximum(ts)
    stats[10] = skewness(ts)
    stats[11] = kurtosis(ts)
    stats[12] = mad(ts, normalize=true)
    xout .=stats
end

function decompose(cube, num_imfs=6)
    indims = InDims("Time")
    imfax = CategoricalAxis("IntrinsicModeFunctions", [(Ref("IMF ") .* string.(1:num_imfs))..., "Residual"])
    @show imfax
    @show length(imfax)
    timeax = ESDL.getAxis("Time", cube)
    od = OutDims(imfax, timeax)

    dateint = intify(cube)
    mapCube(cubeemd!, cube, dateint, num_imfs,  indims=indims, outdims=od)
end

function rqats(cube::YAXArrays.Cubes.YAXArray, dist=1)
    indims = InDims("Time")
    r = rqa(RecurrenceMatrix(rand(10), 1))
    rqanames = string.(keys(r))
    @show rqanames
    rqaaxis = CategoricalAxis("RQA Metrics", collect(rqanames))
    @show rqaaxis
    od = OutDims(rqaaxis)
    @show dist
    @show od
    mapCube(crqa!, cube, indims=indims, outdims=od)
end

function crqa!(xout, xin)
    ts = collect(skipmissing(xin))
    rp = RecurrenceMatrix(ts, 1)
    xout .= values(rqa(rp))
end

function cubeemd!(xout, xin, times, num_imfs)
    ts = collect(skipmissing(xin))
    ind = .!ismissing.(xin)
    dateint = times[ind]
    @show size(xout)
    imfs = EMD.ceemd(ts, dateint, num_imfs=num_imfs)
    @show length(imfs)
    for i in 1:length(imfs)

        xout[i,:] .= imfs[i]
    end
end

function lombscargle(cube, kwargs...)
    indims = InDims("Time")
    lombax = CategoricalAxis("LombScargle", ["Number of Frequencies", "Period with maximal power", "Maximal Power"])
    @show cube
    dateint = intify(cube)
    od = OutDims(lombax)
    mapCube(clombscargle, cube, dateint, indims=indims, outdims=od)
end

function clombscargle(xout, xin, times)
    ind = .!ismissing.(xin)
    ts = collect(nonmissingtype(eltype(xin)), xin[ind])
    dateint = times[ind]
    if length(ts) < 10
        @show length(ts)
        xout .= missing
        return
    end
    pl = LombScargle.plan(dateint, ts)
    #@show pl
    pgram = LombScargle.lombscargle(pl)
    lsperiod= findmaxperiod(pgram)
    lspower = findmaxpower(pgram)
    lsnum = LombScargle.M(pgram)
    #@show lsperiod, lspower
    #@show findmaxfreq(pgram), findmaxpower(pgram)
    xout .= [lsnum, lsperiod[1], lspower]
end

function ctslength(xin)
    ind = .!ismissing.(xin)
    return count(ind)
end

tslength(cube) = mapslices(ctslength, cube, dims="Time")

function intify(cube)
    timeax = ESDL.getAxis("Time", cube)
    x = timeax.values
    datediff = Date.(x) .- Date(x[1])
    dateint = getproperty.(datediff, :value)
end


text:
import LinearAlgebra

################################################################################
#                          Negative log-likelihood                             #
################################################################################

"""
Negative log-likelihood. ŷ should be a vector of normalized log probabilities.
"""
function nllloss(ŷ::AbstractMatrix, c::AbstractVector{<:Real})
    # loss = @tullio (+) r = ŷ[c[k], k]
    loss = LinearAlgebra.tr(NNlib.gather(ŷ', c))
    return -loss / length(c)
end

function ∇nllloss(Δ::Real, ŷ::AbstractMatrix, c::AbstractVector{<:Real})
    dŷ = zero(ŷ)
    val = -Δ / length(c)
    @tullio dŷ[c[k], k] = val
end

function ChainRulesCore.rrule(::typeof(nllloss), ŷ::AbstractMatrix, c::AbstractVector{<:Real})
    y = nllloss(ŷ, c)
    function nllloss_pullback(Δ)
        return NoTangent(), ∇nllloss(unthunk(Δ), ŷ, c), NoTangent()
    end
    return y, nllloss_pullback
end

################################################################################
#                               Cross Entropy                                  #
################################################################################

crossentropyloss(x::AbstractMatrix, c::Union{AbstractMatrix{<:Real}, AbstractVector{<:Real}}) =
    nllloss(logsoftmax(x), c)

################################################################################
#                                 Mean Squared Loss                            #
################################################################################

mseloss(inp::AbstractArray, target::AbstractArray) = mean((inp .- target) .^ 2)


text:
using memoJuliaClassJa
using Test

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


text:
# ---
# title: 931. Minimum Falling Path Sum
# id: problem931
# author: zhwang
# date: 2022-02-23
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-falling-path-sum/description/>
# hidden: true
# ---
#
# Given a **square** array of integers `A`, we want the **minimum** sum of a
# _falling path_ through `A`.
#
# A falling path starts at any element in the first row, and chooses one element
# from each row.  The next row's choice must be in a column that is different
# from the previous row's column by at most one.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[1,2,3],[4,5,6],[7,8,9]]
#     Output: 12
#     Explanation:
#     The possible falling paths are:
#
#
#   * `[1,4,7], [1,4,8], [1,5,7], [1,5,8], [1,5,9]`
#   * `[2,4,7], [2,4,8], [2,5,7], [2,5,8], [2,5,9], [2,6,8], [2,6,9]`
#   * `[3,5,7], [3,5,8], [3,5,9], [3,6,8], [3,6,9]`
#
# The falling path with the smallest sum is `[1,4,7]`, so the answer is `12`.
#
#
#
# **Constraints:**
#
#   * `1 <= A.length == A[0].length <= 100`
#   * `-100 <= A[i][j] <= 100`
#
#
## @lc code=start
using LeetCode

function min_falling_path_sum(matrix::Vector{Vector{Int}})::Int
    n = length(matrix)
    function f(l1, l2)
        return [
            num + minimum(l1[max(i - 1, 1):min(i + 1, n)]) for (i, num) in enumerate(l2)
        ]
    end
    return minimum(foldl(f, matrix))
end

## @lc code=end


text:
#=
 Explicit MDP model of an urban environment with one pedestrian
(extension of the single crosswalk model)
=#

# State type
struct PedMDPState
    crash::Bool
    ego::VehicleState
    ped::VehicleState
end

# copy b to a
function Base.copyto!(a::PedMDPState, b::PedMDPState)
    a.crash = b.crash
    a.ego = b.ego
    a.ped = b.ped
end

function Base.hash(s::PedMDPState, h::UInt64 = zero(UInt64))
    return hash(s.crash, hash(s.ego, hash(s.ped, h)))
end

function Base.:(==)(a::PedMDPState, b::PedMDPState)
    return a.crash == b.crash && a.ego == b.ego && a.ped == b.ped
end

# Action type
const PedMDPAction = UrbanAction

@with_kw mutable struct PedMDP <: MDP{PedMDPState, PedMDPAction}
    env::UrbanEnv = UrbanEnv(params=UrbanParams(nlanes_main=1,
                     crosswalk_pos =  [VecSE2(6, 0., pi/2), VecSE2(-6, 0., pi/2), VecSE2(0., -5., 0.)],
                     crosswalk_length =  [10.0, 10., 10.0],
                     crosswalk_width = [4.0, 4.0, 3.1],
                     stop_line = 22.0))
    ego_type::VehicleDef = VehicleDef()
    ped_type::VehicleDef = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)
    ped_model::DriverModel = ConstantPedestrian()
    max_acc::Float64 = 2.0
    pos_res::Float64 = 1.0
    vel_res::Float64 = 2.0
    vel_ped_res::Float64 = 1.0
    ego_start::Float64 = env.params.stop_line - ego_type.length/2
    ego_goal::LaneTag = LaneTag(2,1)
    off_grid::VecSE2 = VecSE2(UrbanEnv().params.x_min+VehicleDef().length/2, env.params.y_min+VehicleDef().width/2, 0)
    ΔT::Float64 = 0.5
    ped_birth::Float64 = 0.3
    collision_cost::Float64 = -1.0
    action_cost::Float64 = 0.
    goal_reward::Float64 = 1.
    γ::Float64 = 0.95
    _ped_grid::Dict{LaneTag, RectangleGrid{3}} = init_ped_grid(env, pos_res, vel_ped_res)
    _car_grid::Dict{LaneTag, RectangleGrid{2}} = init_car_grid(env, pos_res, vel_res)
    _l_grid::Dict{LaneTag, RectangleGrid{1}} = init_l_grid(env, pos_res)
    _v_grid::RectangleGrid{1} = init_v_grid(env, vel_res)
    _collision_checker::Dict{Tuple{Vararg{VehicleState, 3}}, Bool} = Dict{Tuple{Vararg{VehicleState, 3}}, Bool}()
end

### REWARD MODEL ##################################################################################

function POMDPs.reward(mdp::PedMDP, s::PedMDPState, a::PedMDPAction, sp::PedMDPState)
    r = mdp.action_cost
    ego = sp.ego
    if sp.crash
        r += mdp.collision_cost
    end
    if ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
       get_lane(mdp.env.roadway, ego).tag == mdp.ego_goal
        r += mdp.goal_reward
    end
    return r
end

function POMDPs.isterminal(mdp::PedMDP, s::PedMDPState)
    if s.crash
        return true
    elseif s.ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
       get_lane(mdp.env.roadway, s.ego).tag == mdp.ego_goal
       return true
   end
   return false
end

## Helpers

function POMDPs.discount(mdp::PedMDP)
    return mdp.γ
end

POMDPs.actions(mdp::PedMDP) = [PedMDPAction(-4.0), PedMDPAction(-2.0), PedMDPAction(0.0), PedMDPAction(2.0)]

function POMDPs.actionindex(mdp::PedMDP, action::PedMDPAction)
    if action.acc == -4.0
        return 1
    elseif action.acc == -2.0
        return 2
    elseif action.acc == 0.
        return 3
    else
        return 4
    end
end

function get_off_the_grid(mdp::PedMDP)
    return VehicleState(mdp.off_grid, mdp.env.roadway, 0.)
end

# create labels for model checking
function labeling(mdp::PedMDP)
    labels = Dict{PedMDPState, Vector{String}}()
    for s in ordered_states(mdp)
        if s.crash
            labels[s] = ["crash"]
        elseif s.ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
            get_lane(mdp.env.roadway, s.ego).tag == mdp.ego_goal
            labels[s] = ["goal"]
        end
    end
    return labels
end

function POMDPs.convert_s(::Type{Vector{Float64}}, s::PedMDPState, mdp::PedMDP)
    n_features = 9 # x_ego, y_ego, th_ego, v_ego, x_ped, y_ped, th_ped, v_ped, crash
    z = zeros(n_features)
    z[1] = s.ego.posG.x / abs(mdp.env.params.x_min)
    z[2] = s.ego.posG.y / abs(mdp.env.params.y_min)
    z[3] = s.ego.posG.θ / π
    z[4] = s.ego.v / abs(mdp.env.params.speed_limit)
    z[5] = s.ped.posG.x / abs(mdp.env.params.x_min)
    z[6] = s.ped.posG.y / abs(mdp.env.params.y_min)
    z[7] = s.ped.posG.θ / π
    z[8] = s.ped.v / abs(mdp.env.params.speed_limit)
    z[9] = float(s.crash)
    return z
end

function POMDPs.convert_s(::Type{PedMDPState}, z::V, mdp::PedMDP) where V<:AbstractArray{Float64}
    n_features = 4
    @assert length(z) == n_features*2 + 1
    ego_x = z[1]*abs(mdp.env.params.x_min)
    ego_y = z[2]*abs(mdp.env.params.y_min)
    ego_θ = z[3]*π
    ego_v = z[4]*abs(mdp.env.params.speed_limit)
    ego = VehicleState(VecSE2(ego_x, ego_y, ego_θ), mdp.env.roadway, ego_v)
    ped_x = z[5]*abs(mdp.env.params.x_min)
    ped_y = z[6]*abs(mdp.env.params.y_min)
    ped_θ = z[7]*π
    ped_v = z[8]*abs(mdp.env.params.speed_limit)
    ped = VehicleState(VecSE2(ped_x, ped_y, ped_θ), mdp.env.ped_roadway, ped_v)
    collision = Bool(z[9])
    return PedMDPState(collision, ego, ped)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
function astar(graph, graphWidth, startNode, goal)

    endNode = VertexStruct(goal, 0, goal)

    openList = PriorityQueue()
    closedList = []

    olLock = Threads.Mutex()
    clLock = Threads.Mutex()

    enqueue!(openList, VertexStruct(startNode, 0, startNode), 0)

    # distribute search across the (4) threads
    Threads.@threads for thread in 1:Threads.nthreads()

        # check if the next vertex is the end vertex
        while !isEnd(openList, goal, olLock)

            # get the next vertex off the queue
            lock(olLock)
            try
                current = peek(openList)
            catch
                unlock(olLock)
                continue
            end
            current = dequeue!(openList)
            unlock(olLock)

            # create vertex for the current to keep track of cost and parent as well
            currentNode = VertexStruct(current.key, current.cost, current.parent)

            # LOCK - can't compare and swap because it's an array of VertexStruct
            # add the element to closedList if it's not already there
            lock(clLock)
            exists = in(currentNode, closedList)
            if !exists
                push!(closedList, currentNode)
            end
            unlock(clLock)

            neighborVals = LightGraphs.neighbors(graph, currentNode.key)
            neighbors = []
            for neighborVal in neighborVals
                # create vertexStructs out of neighbors, and add them to an array
                push!(neighbors, VertexStruct(neighborVal, 0, currentNode.key))
            end
            for neighbor in neighbors
                cost = currentNode.cost + 1
                # lock the section to avoid race conditions
                lock(olLock)
                lock(clLock)
                if !in(neighbor, keys(openList)) && isnothing(findfirst((vertex -> vertex.key == neighbor.key), closedList))
                    neighbor.cost = cost
                    #  put the new value onto the openList to flag as "reachable"
                    enqueue!(openList, neighbor, cost+heuristic(endNode.key, neighbor.key, graphWidth))
                    neighbor.parent = currentNode.key
                end
                unlock(olLock)
                unlock(clLock)
            end
        end
    end

    finalVertex = dequeue!(openList)

    path = []
    currentVertex = finalVertex
    push!(path, currentVertex)

    while currentVertex.key != startNode
        parentVertex = closedList[findfirst((vertex -> vertex.key == currentVertex.parent), closedList)]
        push!(path, parentVertex)
        currentVertex = parentVertex
    end

    path = reverse(path)
    pathEdges = Array{Edge{Int64}}(undef, 0)
    length = size(path)[1]
    prev = path[1]

    for i = 2:length
        push!(pathEdges, Edge(prev.key, path[i].key))
        prev = path[i]
    end

    return pathEdges
end

# manhattan distance - no diagonal movement
function heuristic(goal, current, graphWidth)

    x1 = (current-1) % graphWidth
    x2 = (goal-1) % graphWidth
    y1 = div((current-1), graphWidth)
    y2 = div((goal-1), graphWidth)
    return abs(x1-x2) + abs(y1-y2)

end

function isEnd(vertexArray, goalVertex, arrayLock)
    lock(arrayLock)
    isEnd = ((peek(vertexArray)[1]).key == goalVertex)
    unlock(arrayLock)
    return isEnd
end

isnothing(::Any) = false
isnothing(::Nothing) = true

mutable struct VertexStruct
    key::Int64
    cost::Int64
    parent::Int64
end


text:
# https://arxiv.org/abs/2012.06678

using Flux
using Flux: @functor
using Parameters: @with_kw
using FastAI
import FastAI.Models.tabular_embedding_backbone
include("../layers/transformer_block.jl")
include("../layers/categorical_embeddings.jl")

# Struct to define hyperparameters
@with_kw mutable struct TabTransfortmerArgs
    embedding_dims::Vector{Tuple{Int64, Int64}}
    mha_input_dims::Int64
    cont_input_dim::Int64
    output_dim::Int64 = 2
    lr::Float64 = 3e-4		# learning rate
    epochs::Int64 = 10        # number of epochs
    use_cuda::Bool = true   # use gpu (if cuda available)
    dropout::Float64 = 0.10 # dropout from dense layers
    hidden_dims::Vector{Int64} = [128, 64] # Size of hidden layers
    embedding_dropout::Float64 = 0.10 # dropout for categorical embeddin
    batchnorm::Bool = true # batchnorm on dense layers
    linear_first::Bool = true # linear layer before or after batch norm
    activation = Flux.relu
    output_activation = sigmoid
    mha_heads::Int64 = 8
    mha_head_dims::Int64 = 32
    transformer_dropout::Float64 = 0.1
    transformer_dense_hidden_dim::Int64 = 64
    transformer_blocks::Int64 = 2
end

"""
Create a series of dense layers following the
embedding and continious data vcat

https://github.com/FluxML/FastAI.jl/blob/master/src/models/tabularmodel.jl
"""
function dense_layers(args)
    layers = []
    first_layer = linbndrop(args.mha_input_dims + args.cont_input_dim, first(args.hidden_dims);
                    use_bn=args.batchnorm, p=args.dropout, lin_first=args.linear_first,
                    act=args.activation)
    push!(layers, first_layer)

    for (isize, osize) in zip(args.hidden_dims[1:(end-1)], args.hidden_dims[2:end])
        layer = linbndrop(isize, osize; use_bn=args.batchnorm, p=args.dropout,
                        lin_first=args.linear_first, act=args.activation)
        push!(layers, layer)
    end

    return Chain(layers...)

end

struct TabTransformer
    cat_embeddings
    cat_backbone
    cont_backbone
    dense
    output
end
@functor TabTransformer

TabTransformer(args::TabTransfortmerArgs) = TabTransformer(
    tabular_embedding_backbone(args.embedding_dims, args.embedding_dropout),
    Chain([TransformerBlock(args) for i in 1:args.transformer_blocks]...,
            Flux.flatten),
    BatchNorm(args.cont_input_dim),
    dense_layers(args),
    Dense(args.hidden_dims[end], args.output_dim, args.output_activation)
)

function (t::TabTransformer)(X_cat, X_cont)
    h_cat = t.cat_embeddings(X_cat...)
    h_cat = t.cat_backbone(reshape(h_cat, (size(h_cat, 1), 1, size(h_cat, 2))))
    h_cont = t.cont_backbone(X_cont)
    h2 = t.dense(vcat(h_cat, h_cont))
    t.output(h2)
end


text:
using JLBoost
using Test
using DataFrames
using Random

X = rand(100, 50)

#y = 1.0*(rand(100) .> 0.5)
y = (sum(X', dims=[1]).>25)[1,:];

df = DataFrame(X, :auto);
df[!,"y"] = y;

xgb = jlboost(df, "y"; verbose=false);


text:
using Documenter
import Mosek

makedocs(
    format = :html,
    sitename = "Mosek",
    pages = [
        "Index"         => "index.md",
        "API Reference" => "Mosek-Functions.md",
        "General Convex API" => "geco.md",
        "Solver parameters"    => "parameters.md",
        "Symbolic values"    => "enums.md"
    ]
)

deploydocs(
    repo   = "github.com/JuliaOpt/Mosek.jl.git",
    target = "build",
    osname = "linux",
    julia  = "0.6",
    deps   = nothing,
    make   = nothing
)


text:
using Test

include("high-score.jl")

@testset "Adding players" begin
    @testset "Add player with score to empty score dict" begin
        scores = Dict{String, Int}()
        add_player!(scores, "Kalei", 233_436)
        @test scores == Dict(
            "Kalei" => 233_436,
        )
    end

    @testset "Add two players with scores to empty score dict" begin
        scores = Dict{String, Int}()
        add_player!(scores, "Léna", 901_806)
        add_player!(scores, "Ljuba", 302_799)
        @test scores == Dict(
            "Léna" => 901_806,
            "Ljuba" => 302_799,
        )
    end

    @testset "Add player without score to empty score dict" begin
        scores = Dict{String, Int}()
        add_player!(scores, "Nicolet")
        @test scores == Dict(
            "Nicolet" => 0,
        )
    end

    @testset "Add two players without score to empty dict" begin
        scores = Dict{String, Int}()
        add_player!(scores, "Dina")
        add_player!(scores, "Yasu")
        @test scores == Dict(
            "Dina" => 0,
            "Yasu" => 0,
        )
    end
end

@testset "Removing players" begin
    @testset "Removing a player from empty score dict results in empty score dict" begin
        scores = Dict{String, Int}()
        remove_player!(scores, "Waldemar")
        @test isempty(scores)
    end

    @testset "Removing a player after adding results in empty score dict" begin
        scores = Dict{String, Int}()
        add_player!(scores, "Arisha")
        remove_player!(scores, "Arisha")
        @test isempty(scores)
    end

    @testset "Removing the first player after adding two results in dict with remaining player" begin
        scores = Dict{String, Int}()
        add_player!(scores, "Fritjof")
        add_player!(scores, "Valèria")
        remove_player!(scores, "Fritjof")
        @test scores == Dict(
            "Valèria" => 0,
        )
    end

    @testset "Removing the second player after adding two results in dict with remaining player" begin
        scores = Dict{String, Int}()
        add_player!(scores, "Fritjof")
        add_player!(scores, "Valèria")
        remove_player!(scores, "Valèria")
        @test scores == Dict(
            "Fritjof" => 0,
        )
    end
end

@testset "Updating scores" begin
    @testset "Updating score for non-existent player initializes value" begin
        scores = Dict{String, Int}()
        update_score!(scores, "Dov", 721_957)
        @test scores == Dict(
            "Dov" => 721_957,
        )
    end

    @testset "Updating score for existing player without score replaces previous score" begin
        scores = Dict(
            "Sabriyya" => 0,
        )
        update_score!(scores, "Sabriyya", 560_806)
        @test scores == Dict(
            "Sabriyya" => 560_806,
        )
    end

    @testset "Updating score for existing player with score replaces previous score" begin
        scores = Dict(
            "Edvard" => 151_216,
        )
        update_score!(scores, "Edvard", 600_259)
        @test scores == Dict(
            "Edvard" => 600_259,
        )
    end

    @testset "Updating score for existing player in a dict with several players with score replaces player's previous score" begin
        scores = Dict(
            "Əli" => 977_076,
            "Lisbet" => 288_100,
        )
        update_score!(scores, "Lisbet", 989_376)
        @test scores == Dict(
            "Əli" => 977_076,
            "Lisbet" => 989_376,
        )
    end
end

@testset "Resetting scores" begin
    @testset "Resetting score for existing player sets score to 0" begin
        scores = Dict(
            "Thea" => 307_261,
        )
        reset_score!(scores, "Thea")
        @test scores == Dict(
            "Thea" => 0,
        )
    end

    @testset "Resetting score for existing player in a dict with several players sets player's score to 0" begin
        scores = Dict(
            "Anelia" => 307_261,
            "Robert" => 785_497,
        )
        reset_score!(scores, "Anelia")
        @test scores == Dict(
            "Anelia" => 0,
            "Robert" => 785_497,
        )
    end

    @testset "Resetting score for non-existent player sets player score to 0" begin
        scores = Dict{String, Int}()
        reset_score!(scores, "Erwin")
        @test scores == Dict(
            "Erwin" => 0,
        )
    end
end

@testset "List of players" begin
    @testset "Score board with one entry gives one result" begin
        scores = Dict(
            "Rodolfito" => 820_876
        )
        @test players(scores) == ["Rodolfito"]
    end

    @testset "Score board with multiple entries gives results in unknown order" begin
        scores = Dict(
            "Naheed" => 97_506,
            "Ivona" => 837_854,
            "Iztok" => 160_855,
            "Sara" => 936_371,
            "Udo" => 25,
        )
        @test sort(players(scores)) == ["Ivona", "Iztok", "Naheed", "Sara", "Udo"]
    end

    @testset "Empty score board results in empty player list" begin
        scores = Dict()
        @test isempty(players(scores))
    end
end


text:
using LearningStrategies
import LearningStrategies: update!

mutable struct GPModel{N, C}
    fitness::Function
    population::Vector{DecisionTree{N, C}}
    cache::Vector{DecisionTree{N, C}}
    population_fitnesses::Vector{Float64}
    cache_fitnesses::Vector{Float64}
end

function GPModel(fitness, population)
    fitnesses = fitness.(population)
    GPModel(fitness, population, similar(population), fitnesses, similar(fitnesses))
end

mutable struct SSGPModel{N, C}
    fitness::Function
    population::Vector{DecisionTree{N, C}}
    population_fitnesses::Vector{Float64}
end

function SSGPModel(fitness, population)
    fitnesses = fitness.(population)
    SSGPModel(fitness, population, fitnesses)
end

struct GPModelSolver{N, C} <: LearningStrategy
    min_depth::Int
    max_depth::Int
    tournament_size::Int
    crossover_probability::Float64
    mutation_probability::Float64
    tree_sampler::TreeSampler{N, C}
    rng::AbstractRNG
end

function selection(parents, parent_fitnesses, s)
    candidates = rand(s.rng, eachindex(parents, parent_fitnesses), s.tournament_size)
    parents[indmax(parent_fitnesses[candidates])]
end

function mate(p, q, s)
    # Choose a random subtree from `p` and replace it by a random child of `q`
    newtree, _ = randsplit((subtree) -> randchild(s.rng, q), s.rng, p)
    newtree
end

function crossover(parent₁, parent₂, s)
    p = s.crossover_probability
    r = rand(s.rng)
    if r ≤ p
        candidate = mate(parent₁, parent₂, s)
        while !(s.min_depth ≤ treedepth(candidate) ≤ s.max_depth)
            candidate = mate(parent₁, parent₂, s)
        end

        return candidate
    elseif p < r ≤ (1 - p) / 2
        parent₁
    else # (1 - p) / 2 < r ≤ 1
        parent₂
    end
end

function splice(i, s)
    # Choose a random subtree from `i` and replace it by a random tree
    newtree, _ = randsplit((subtree) -> rand(s.rng, s.tree_sampler), i)
    newtree
end

function mutate(individual, s)
    if rand(s.rng) ≤ s.mutation_probability
        candidate = splice(individual, s)
        while !(s.min_depth ≤ treedepth(candidate) ≤ s.max_depth)
            candidate = splice(individual, s)
        end

        return candidate
    else
        return individual
    end
end

function update!(model::GPModel, s::GPModelSolver)
    parents = model.population
    children = model.cache
    parent_fitnesses = model.population_fitnesses
    child_fitnesses = model.cache_fitnesses
    fitness = model.fitness

    # breed
    for i in eachindex(children, child_fitnesses)
        p₁ = selection(parents, parent_fitnesses, s)
        p₂ = selection(parents, parent_fitnesses, s)
        child = crossover(p₁, p₂, s)
        children[i] = mutate(child, s)
        child_fitnesses[i] = fitness(child)
    end

    # update: swap actual values with caches
    model.population = children
    model.population_fitnesses = child_fitnesses
    model.cache = parents
    model.cache_fitnesses = parent_fitnesses
end

function update!(model::SSGPModel, s::GPModelSolver)
    population = model.population
    fitnesses = model.population_fitnesses
    fitness = model.fitness

    for _ = eachindex(population)
        # breed
        p₁ = selection(population, fitnesses, s)
        p₂ = selection(population, fitnesses, s)
        child = crossover(p₁, p₂, s)
        offspring = mutate(child, s)
        offspring_fitness = fitness(offspring)

        # replace worst individual
        candidate = indmin(fitnesses)
        population[candidate] = offspring
        fitnesses[candidate] = offspring_fitness
    end
end

function rungp{N, C}(fitness, psize::Int, sampler::TreeSampler{N, C}, maxiter::Int;
                     tracer = Tracer(Void, (m, i) -> nothing, typemax(Int)),
                     breaker = Breaker((m, i) -> false),
                     min_depth = 1, max_depth = 20,
                     tournament_size = 7,
                     mutation_rate = 0.5, crossover_rate = 0.5,
                     depth_penalty = 2.0, size_penalty = 0.5,
                     rng = Base.GLOBAL_RNG, verbose = true)
    initial_population = rand(rng, sampler, psize)
    model = GPModel(float ∘ fitness, initial_population)
    solver = GPModelSolver(min_depth, max_depth, tournament_size, crossover_rate, mutation_rate,
                           sampler, rng)
    iteration_control = verbose ? Verbose(MaxIter(maxiter)) : MaxIter(maxiter)
    learn!(model, strategy(solver, iteration_control, tracer, breaker))

    return model.population, tracer
end

function runssgp{N, C}(fitness, psize::Int, sampler::TreeSampler{N, C}, maxiter::Int;
                       tracer = Tracer(Void, (m, i) -> nothing, typemax(Int)),
                       breaker = Breaker((m, i) -> false),
                       min_depth = 1, max_depth = 20,
                       tournament_size = 7,
                       mutation_rate = 0.5, crossover_rate = 0.5,
                       depth_penalty = 2.0, size_penalty = 0.5,
                       rng = Base.GLOBAL_RNG,
                       verbose = true)
    initial_population = rand(rng, sampler, psize)
    model = SSGPModel(float ∘ fitness, initial_population)
    solver = GPModelSolver(min_depth, max_depth, tournament_size, crossover_rate, mutation_rate,
                           sampler, rng)
    iteration_control = verbose ? Verbose(MaxIter(maxiter)) : MaxIter(maxiter)
    learn!(model, strategy(solver, iteration_control, tracer, breaker))

    return model.population, tracer
end

max_stars_repo_name: phipsgabler/soft-computing-final

max_issues_repo_name: phipsgabler/soft-computing-final

max_forks_repo_name: phipsgabler/soft-computing-final


text:
using AbstractFFTs
using FFTW
using LinearAlgebra
using SubpixelRegistration: fourier_shift!

"""
    fourier_lucky_image(cube; dims, q, register=:dft, kwargs...)

Perform Fourier lucky imaging following the algorithm in Garrel, Guyon, and Baudoz (2012).[^1] This technique provides superior sharpness compared to classic lucky imaging *for the same selection percentage*. This means for a desired angular resolution (or Strehl ratio) a higher signal-to-noise ratio can be achieved.

[^1]: Vincent Garrel et al 2012 ["A Highly Efficient Lucky Imaging Algorithm: Image Synthesis Based on Fourier Amplitude Selection"](https://iopscience.iop.org/article/10.1086/667399) PASP 124 861

# Registration

Before processing the Fourier transform of the cube, the frames need coregistered to avoid phase ramps in the Fourier transform modulus. The following algorithms are available for registering the cube, set with the `register` keyword argument. Extra options can be provided to `kwargs...`

* `:dft` - use `SubpixelRegistration.jl` to register the frames using their phase offsets. The reference frame will be the one with the highest metric, and keyword arguments like `upsample_factor` and `refidx` can be passed directly. By default the `refidx` will be set to the frame with the highest peak flux.
* `:peak` - register to maximum value
* `:com` - register to center of mass

# Keyword arguments

* `dims` - the dimension along which to perform lucky imaging (required)
* `q` - the selection quantile (required)
* `register` - the method used for registration. Default is `:dft`
* `maxfreq` - if provided, will exclude frequencies higher than `maxfreq * maximum(fftfreqs)`. In other words, a value of 1 uses the full Fourier transform, but a value of 0.5 will low-pass filter
* `kwargs...` - additional keyword arguments will be passed to the `register` method (e.g., `upsample_factor`)

# Examples

```julia
julia> cube = # load data ...

julia> res = fourier_lucky_image(cube; dims=3, q=0.5, upsample_factor=10);
```
"""
function fourier_lucky_image(cube::AbstractArray{T,3}; dims, q, register=:dft, maxfreq=1, metric=nothing, kwargs...) where T
    # register cube
    if register === :dft
        # to help coregistration, choose refidx from frame with highest flux
        refidx = @compat findmax(maximum, eachslice(cube; dims))[2]
        refframe = selectdim(cube, dims, refidx)
        # get shift to center refframe
        refshift = center_of_mass(refframe) .- center(refframe)
        reffreq = fft(refframe)
        # go ahead and shift it now
        fourier_shift!(reffreq, refshift)
        # don't actually register yet, do this during loop
        registered = cube
    else
        registered = copy(cube)
        @inbounds for didx in axes(registered, dims)
            frame = selectdim(registered, dims, didx)
            if register === :peak
                index = argmax(frame).I
            elseif register === :com
                index = center_of_mass(frame)
            end
            shift = index .- center(frame)
            # update frame in-place (it is a view)
            frame .= fourier_shift(frame, shift)
        end
    end

    # plan fft using first slice
    first_frame = selectdim(registered, dims, firstindex(registered, dims))
    plan = plan_fft(first_frame)
    # set up work arrays to avoid allocating inside loop
    mean_freq = zeros(Complex{T}, size(first_frame))
    # tmp_mod = similar(first_frame)
    norm_value = size(cube, dims)
    # figure out frequency mask ahead of time

    freqs = map(fftfreq, size(mean_freq))
    _maxfreq = norm(maxfreq .* maximum.(freqs))
    mask = @. hypot(freqs[1], freqs[2]')  > _maxfreq

    if !isnothing(metric)
        frame_gen = eachslice(cube; dims=dims)
        if metric === :peak
            _metric = map(maximum, frame_gen)
        elseif metric === :mean
            _metric = map(mean, frame_gen)
        else
            _metric = map(metric, frame_gen)
        end
        cutoff = quantile(_metric, q)
    end

    # create frequency cube
    freq_cube = similar(registered, complex(T))
    @inbounds for didx in axes(registered, dims)
        frame = selectdim(registered, dims, didx)
        freq_frame = selectdim(freq_cube, dims, didx)
        freq_frame .= plan * frame
        # if we are doing DFT registration we can just do that during this step to
        # save of time spent doing FFTs
        if register === :dft && didx != refidx
            result = phase_offset(plan, reffreq, freq_frame; kwargs...)
            fourier_shift!(freq_frame, result.shift, result.phasediff)
        end
    end

    mod_cube = map(abs2, freq_cube)
    mod_cutoff = mapslices(s -> quantile(s, q), mod_cube; dims)
    freq_cube[mod_cube .< mod_cutoff] .= zero(eltype(freq_cube))
    mean_freq = mean(freq_cube; dims)

    # mean_freq = mapslices(freq_cube; dims) do slice
    #     # get modulus cutoff from quantile
    #     @. mod_vec = abs2(slice)
    #     cutoff = quantile(mod_vec, q)
    #     mask = mod_vec .≥ cutoff
    #     return mean(slice[mask])
    # end
    # @inbounds for didx in axes(freq_cube, dims)
    #     frame_freq = selectdim(freq_cube, dims, didx)

    #     # low-pass filter
    #     if isnothing(metric) || _metric[didx] < cutoff
    #         frame_freq[mask] .= zero(eltype(frame_freq))
    #     end

    #     # get selection cutoff from remaining frequencies
    #     @. tmp_mod = abs(frame_freq)
    #     freq_cutoff = quantile(vec(tmp_mod), q)
    #     # update running mean of frequencies
    #     @. mean_freq[tmp_mod ≥ freq_cutoff] += frame_freq[tmp_mod ≥ freq_cutoff] / norm_value
    # end
    # ifft
    return real(plan \ view(mean_freq, :, :, 1))
end


text:
@testset "Misc" begin
  include("Misc/psibound.jl")
  include("Misc/Primes.jl")
  include("Misc/Places.jl")
  include("Misc/meataxe.jl")
  include("Misc/stable_subgroups.jl")
  include("Misc/NumberField.jl")
end


text:
#=
Functions to compute quadratic sums

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date : 2014-08-19
=#

@doc doc"""
Computes the expected discounted quadratic sum

```math
    q(x_0) = \mathbb{E} \sum_{t=0}^{\infty} \beta^t x_t' H x_t
```

Here ``{x_t}`` is the VAR process ``x_{t+1} = A x_t + C w_t`` with ``{w_t}``
standard normal and ``x_0`` the initial condition.

##### Arguments
- `A::Union{Float64, Matrix{Float64}}` The `n x n` matrix described above (scalar)
  if `n = 1`
- `C::Union{Float64, Matrix{Float64}}` The `n x n` matrix described above (scalar)
  if `n = 1`
- `H::Union{Float64, Matrix{Float64}}` The `n x n` matrix described above (scalar)
  if `n = 1`
- `beta::Float64`: Discount factor in `(0, 1)`
- `x_0::Union{Float64, Vector{Float64}}` The initial condtion. A conformable
  array (of length `n`) or a scalar if `n = 1`

##### Returns

- `q0::Float64` : Represents the value ``q(x_0)``

##### Notes

The formula for computing ``q(x_0)`` is ``q(x_0) = x_0' Q x_0 + v`` where

- ``Q`` is the solution to ``Q = H + \beta A' Q A`` and
- ``v = \frac{trace(C' Q C) \beta}{1 - \beta}``

"""
function var_quadratic_sum(A::ScalarOrArray, C::ScalarOrArray, H::ScalarOrArray,
                           bet::Real, x0::ScalarOrArray)
    n = size(A, 1)

    # coerce shapes
    A = reshape([A;], n, n)
    C = reshape([C;], n, n)
    H = reshape([H;], n, n)
    x0 = reshape([x0;], n)

    # solve system
    Q = solve_discrete_lyapunov(sqrt(bet) .* A', H)
    cq = C'*Q*C
    v = tr(cq) * bet / (1 - bet)
    q0 = x0'*Q*x0 + v
    return q0[1]
end

@doc doc"""
Computes the quadratic sum

```math
    V = \sum_{j=0}^{\infty} A^j B A^{j'}
```

``V`` is computed by solving the corresponding discrete lyapunov equation using the
doubling algorithm.  See the documentation of `solve_discrete_lyapunov` for
more information.

##### Arguments

- `A::Matrix{Float64}` : An `n x n` matrix as described above.  We assume in order
  for convergence that the eigenvalues of ``A`` have moduli bounded by unity
- `B::Matrix{Float64}` : An `n x n` matrix as described above.  We assume in order
  for convergence that the eigenvalues of ``B`` have moduli bounded by unity
- `max_it::Int(50)` : Maximum number of iterations

##### Returns

- `gamma1::Matrix{Float64}` : Represents the value ``V``

"""
function m_quadratic_sum(A::Matrix, B::Matrix; max_it=50)
    solve_discrete_lyapunov(A, B, max_it)
end


text:
using DifferentialEquations
# using Plots
include("YfindNP.jl")
using PyCall
pygui(:qt)
using PyPlot
pygui(true)

H=.25
h=1/(2*H)
max_time=50
a=(2+h-2*sqrt(h+1))/h; # paramter a (ratio) in terms of h

#exit
# Q=0.0015000000999000002
# P=0.0411
# Y=0.5780004334287135
#before
# Q=0.7040082229574488
# P=-1.798419378432707
# Y=1.523090142608121
#nowhere
Q=-1.0010302468906507
P=1.7385260494565795
Y=YfindNP(Q,P,H)

u0=[0;P;Q;Y[1];0;0;0]
# u0= [0;0;pert;-Y_0;0;0;0]
tspan = (0.0,max_time)
#10.531140269496431
#4.75752093999673
ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
Centroid_ODE(z,w)=conj( 1 ./(1+z.^2)+1 ./(1+w.^2) )
# refernce
# z=X+iP W=Q+iY . Plane Q=0, P=0
# z=u[1]+i u[2] ; w=u[3]+i u[4]

function Eq_of_M(du,u,p,t)
du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
du[5]=real(Centroid_ODE(u[1]+im*u[2],u[3]+im*u[4]));
du[6]=imag(Centroid_ODE(u[1]+im*u[2],u[3]+im*u[4]))
du[7]=0
return
end
condition(u,t,integrator)= u[7]==20
affect!(integrator) = terminate!(integrator)
condition2(u,t,integrator)=  u[1]
function affect2!(integrator)
integrator.u[7]=integrator.u[7]+1
end
cb1 = DiscreteCallback(condition,affect!)
cb2 =ContinuousCallback(condition2,affect2!,rootfind=true)
cb=CallbackSet(cb2,cb1)

prob = ODEProblem(Eq_of_M,u0,tspan)
t,A=solve(prob, Vern9(),reltol=1e-13,abstol=1e-20,maxiters=1e15,callback=cb)
N=length(t.t)
Z=A[1,:]+im *A[2,:];
W=A[3,:]+im *A[4,:];
Energy=zeros(N)
for k=1:N
    Energy[k]=abs(1/(1+Z[k]^2)-1/(1+W[k]^2))
end

Centeroid=zeros(N);
Centeroid=A[5,:]+im *A[6,:];
lin_impulse=im*(1+a )*ones(N);
z_1pos=.5*(Centeroid+Z+lin_impulse+W);
z_1neg=.5*(Centeroid+Z-lin_impulse-W);
z_2pos=.5*(Centeroid-Z+lin_impulse-W);
z_2neg=.5*(Centeroid-Z-lin_impulse+W);
DDD=(z_2pos+z_1pos-z_2neg-z_1neg)/2;

z_1posX=real(z_1pos)
z_1posY=imag(z_1pos)

z_2posX=real(z_2pos)
z_2posY=imag(z_2pos)

z_1negX=real(z_1neg)
z_1negY=imag(z_1neg)

z_2negX=real(z_2neg)
z_2negY=imag(z_2neg)

figure()
plot(t.t,A[1,:])
uf=zeros(7)
uf[1]=A[1,end]#X
uf[2]=A[2,end] #P
uf[3]=A[3,end] #Q
uf[4]=A[4,end] #Y
uf[5]=A[5,end] #Y
uf[6]=A[6,end]
#
# figure()
# tspan = (0.0,-max_time)
# prob = ODEProblem(Eq_of_M,uf,tspan)
# conditionM(u,t,integrator)= u[7]==20
# cb1 = DiscreteCallback(conditionM,affect!)
# cb2 =ContinuousCallback(condition2, affect2!, rootfind=true)
# cb=CallbackSet(cb2,cb1)
# t,B=solve(prob, Vern9(),reltol=1e-13,abstol=1e-20,maxiters=1e15,callback=cb)
# plot(t.t,B[1,:])
#
#
#
# figure()
# tspan = (0.0,-max_time)
# prob = ODEProblem(Eq_of_M,u0,tspan)
# conditionM(u,t,integrator)= u[7]==20
#
# cb1 = DiscreteCallback(conditionM,affect!)
# cb2 =ContinuousCallback(condition2, affect2!, rootfind=true)
# cb=CallbackSet(cb2,cb1)
# t,B=solve(prob, Vern9(),reltol=1e-13,abstol=1e-20,maxiters=1e15,callback=cb)
# plot(t.t,B[1,:])

# # col=rand(1)
# plot!(z_1pos,linewidth=1,seriescolor=:red)
# plot!(z_1neg, linewidth=1,seriescolor=:blue)
# plot!(z_2pos,linewidth=1)
# plot!(z_2neg,linewidth=1)


text:
module Soss

import Base.rand
using Random
using Reexport: @reexport

@reexport using Distributions
@reexport using StatsFuns
using NamedTupleTools

import MacroTools: prewalk, postwalk, @q, striplines, replace, flatten, @capture
import MLStyle
@reexport using MonteCarloMeasurements

using LazyArrays
using FillArrays

include("core/statement.jl")
include("core/model.jl")
include("core/jointdistribution.jl")
include("core/canonical.jl")
include("core/dependencies.jl")
include("core/toposort.jl")
include("core/weighted.jl")
include("core/utils.jl")

include("distributions/dist.jl")
include("distributions/for.jl")
include("distributions/iid.jl")
include("distributions/mix.jl")
# include("distributions/flat.jl")

include("primitives/rand.jl")
include("primitives/logpdf.jl")
include("primitives/xform.jl")
include("primitives/likelihood-weighting.jl")
# include("primitives/bijectors.jl")

include("transforms/predictive.jl")
include("transforms/markovblanket.jl")
include("transforms/do.jl")

include("symbolic/symbolic.jl")
include("symbolic/codegen.jl")

include("particles.jl")
include("plots.jl")

include("inference/rejection.jl")
include("inference/dynamicHMC.jl")
include("inference/advancedhmc.jl")

# include("weighted.jl")
#
# # include("graph.jl")
# # # include("optim.jl")
# include("importance.jl")
#
# # # include("sobols.jl")
# # # include("fromcube.jl")
# # # include("tocube.jl")
end # module