text:
# using CSV
# using PyPlot
# using Dates

path_45 = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp45/rcp45/"
Name_Projections_45 = readdir(path_45)
path_85 = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp85/rcp85/"
Name_Projections_85 = readdir(path_85)

include("loadfunctions.jl")
include("compare_Present_Future.jl")
include("compare_Present_Future_High_Flows.jl")
include("compare_Present_Future_snow.jl")
include("compare_Present_Future_low_flows.jl")

Area_Catchment_Gailtal = sum([98227533.0, 184294158.0, 83478138.0, 220613195.0])
Area_Catchment_Palten = sum([198175943.0, 56544073.0, 115284451.3])
Area_Catchment_Pitten = 115496400.
Area_Catchment_Silbertal = 100139168.
Area_Catchment_Defreggental = sum([235811198.0, 31497403.0])
Area_Catchment_Pitztal = sum([20651736.0, 145191864.0])

Discharge_Feistritz = CSV.read("/home/sarah/HBVModel/Feistritz/Q-Tagesmittel-214353.csv", header= false, skipto=388, decimal=',', delim = ';', types=[String, Float64])
Discharge_Palten = CSV.read("/home/sarah/HBVModel/Palten/Q-Tagesmittel-210815.csv", header= false, skipto=21, decimal=',', delim = ';', types=[String, Float64])
Discharge_Gailtal = CSV.read("/home/sarah/HBVModel/Gailtal/Q-Tagesmittel-212670.csv", header= false, skipto=23, decimal=',', delim = ';', types=[String, Float64])
Discharge_Silbertal = CSV.read("/home/sarah/HBVModel/Silbertal/Q-Tagesmittel-200048.csv", header= false, skipto=24, decimal=',', delim = ';', types=[String, Float64])
Discharge_Defreggental = CSV.read("/home/sarah/HBVModel/Defreggental/Q-Tagesmittel-212100.csv", header= false, skipto=26, decimal=',', delim = ';', types=[String, Float64])
Discharge_Pitztal = CSV.read("/home/sarah/HBVModel/Pitztal/Q-Tagesmittel-201335.csv", header= false, skipto=23, decimal=',', delim = ';', types=[String, Float64])

Catchment_Name = "Pitztal"
Area_Catchment = Area_Catchment_Pitztal
Discharge_Catchment = Discharge_Pitztal
nr_runs = 300
scaling_factor_discharge = 1
# make folder Projections/Catchment_Name/PastvsFuture/ - (Inputs, Monthly_Discharge, Annual_Discharge, FDC, Budyko, Annual_Max_Discharge/Circular, LowFlows, Snow_Cover, Drought (/none, /summer /winter)

# make plots of tmeperature and precipitation in past and future
#Prec_past, Prec_Future,Temp_past, Temp_Future, Months = plot_Monthly_Temperature_Precipitation(path_45, Catchment_Name)
#Prec_past, Prec_Future,Temp_past, Temp_Future, Months = plot_Monthly_Temperature_Precipitation(path_85, Catchment_Name)

# Prec_statistics_past_45, Prec_statistics_future_45 = monthly_prec_statistics(path_45, Catchment_Name)
# Prec_statistics_past_85, Prec_statistics_future_85 = monthly_prec_statistics(path_85, Catchment_Name)

# plot HYDROGRAPHS ------------
# @time begin
# plot_hydrographs_proj(path_45, Catchment_Name, Area_Catchment, nr_runs)
# plot_hydrographs_proj(path_85, Catchment_Name, Area_Catchment, nr_runs)
# end
# plot_hydrographs_proj(path_45, Catchment_Name, Area_Catchment, 1987, 2087)
# plot_hydrographs_proj(path_85, Catchment_Name, Area_Catchment, 1987, 2087)
# plot_hydrographs_proj(path_45, Catchment_Name, Area_Catchment, 1990, 2090)
# plot_hydrographs_proj(path_85, Catchment_Name, Area_Catchment, 1990, 2090)
# plot_hydrographs_proj(path_45, Catchment_Name, Area_Catchment, 1995, 2095)
# plot_hydrographs_proj(path_85, Catchment_Name, Area_Catchment, 1995, 2095)
# plot_hydrographs_proj(path_45, Catchment_Name, Area_Catchment, 1983, 2083)
# plot_hydrographs_proj(path_85, Catchment_Name, Area_Catchment, 1983, 2083)
# plot_hydrographs_proj(path_45, Catchment_Name, Area_Catchment, 1987, 2087)
# plot_hydrographs_proj(path_85, Catchment_Name, Area_Catchment, 1987, 2087)
# plot_hydrographs_proj(path_45, Catchment_Name, Area_Catchment, 1992, 2092)
# plot_hydrographs_proj(path_85, Catchment_Name, Area_Catchment, 1992, 2092)
# # # ----------------- MONTHLY DISCHARGES ------------------------
# # saves the results to a textfile
# @time begin
# Monthly_Discharge_Change_85, monthly_Discharge_past_85, monthly_Discharge_future_85, months_85 = change_monthly_Discharge(path_85, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_8.5.txt",hcat(months_85, monthly_Discharge_past_85, monthly_Discharge_future_85, Monthly_Discharge_Change_85),',')
# end
# #
# @time begin
# Monthly_Discharge_Change_45, monthly_Discharge_past_45, monthly_Discharge_future_45, months_45 = change_monthly_Discharge(path_45, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_4.5.txt",hcat(months_45, monthly_Discharge_past_45, monthly_Discharge_future_45, Monthly_Discharge_Change_45),',')
# end
# #
# monthly_changes_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_8.5.txt", ',')
# months_85 = monthly_changes_85[:,1]
# monthly_Discharge_past_85 = monthly_changes_85[:,2]
# monthly_Discharge_future_85  = monthly_changes_85[:,3]
# Monthly_Discharge_Change_85  = monthly_changes_85[:,4]
# monthly_changes_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_4.5.txt", ',')
# months_45 = monthly_changes_45[:,1]
# monthly_Discharge_past_45 = monthly_changes_45[:,2]
# monthly_Discharge_future_45  = monthly_changes_45[:,3]
# Monthly_Discharge_Change_45  = monthly_changes_45[:,4]
# #
# plot_changes_monthly_discharge(monthly_Discharge_past_45, monthly_Discharge_future_45, monthly_Discharge_past_85, monthly_Discharge_future_85, months_45, Catchment_Name, nr_runs)
# #
# plot_changes_monthly_discharge_mm(monthly_Discharge_past_45, monthly_Discharge_future_45, monthly_Discharge_past_85, monthly_Discharge_future_85, months_45, Area_Catchment, Catchment_Name)

## -------------- PLOT MONTHLY DISCHARGE PAST REAL DATA ---------------------------------
# Timeseries_Past = collect(Date(1981,1,1):Day(1):Date(2010,12,31))
# Discharge = convert(Matrix, Discharge_Catchment)
# startindex = findfirst(isequal("01.01."*string(1981)*" 00:00:00"), Discharge)
# endindex = findfirst(isequal("31.12."*string(2010)*" 00:00:00"), Discharge)
# Observed_Discharge = Array{Float64,1}[]
# push!(Observed_Discharge, Discharge[startindex[1]:endindex[1],2])
# Observed_Discharge_m3s = Observed_Discharge[1]
# Observed_Discharge = Observed_Discharge_m3s * 1000 / Area_Catchment * (3600 * 24)
# Monthly_Discharge_Observed, Months_Past = monthly_discharge(Observed_Discharge, Timeseries_Past)
# plot()
# for month in 1:12
#     boxplot!(Monthly_Discharge_Observed[findall(x-> x == month, Months_Past)], size=(2000,800), leg=false, color=["red"], left_margin = [5mm 0mm])
# end
# ylabel!("Averaged Monthly Discharge [mm/d]")
# title!("Averaged Measured Monthly Discharge 1981-2010")
# #ylims!((0,40))
# #hline!([0], color=["grey"], linestyle = :dash)
# xticks!([1:12;], ["Jan", "Feb", "Mar", "Apr", "May","Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])
# #xticks!([2:2:24;], ["Jan 8.5", "Feb 8.5", "Mar 8.5", "Apr 8.5", "May 8.5","Jun 8.5", "Jul 8.5", "Aug 8.5", "Sep 8.5", "Oct 8.5", "Nov 8.5", "Dec 8.5"])
# savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/"*Catchment_Name*"_measured_monthly_discharge_mm.png")
#
# # ---------------- ANNUAL DISCHARGE --------------------
#
# relative_change_45, Total_Discharge_Past_45, Total_Discharge_Future_45 = annual_discharge_new(monthly_Discharge_past_45, monthly_Discharge_future_45, Area_Catchment, nr_runs)
# relative_change_85, Total_Discharge_Past_85, Total_Discharge_Future_85 = annual_discharge_new(monthly_Discharge_past_85, monthly_Discharge_future_85, Area_Catchment, nr_runs)

# #can also be calculate directly from discharge data but takes more time
#relative_change_45, Total_Discharge_Past_45, Total_Discharge_Future_45, relative_change_85, Total_Discharge_Past_85, Total_Discharge_Future_85 = change_annual_discharge(path_45, path_85, Catchment_Name, Area_Catchment)
# plot_change_annual_discharge(relative_change_45, Total_Discharge_Past_45, Total_Discharge_Future_45, relative_change_85, Total_Discharge_Past_85, Total_Discharge_Future_85, Area_Catchment, Catchment_Name)

# ----------------------- FDC ------------------------

## the change in discharge of a certian percentile can be calculated
# percentile = 0.1
#discharge_past_45, discharge_future_45= FDC_compare_percentile(path_45, percentile, Catchment_Name)
#discharge_past_85, discharge_future_85= FDC_compare_percentile(path_85, percentile, Catchment_Name)
## or directly plotted
#plot_FDC_Percentile(path_45, path_85, Catchment_Name, percentile)

## or the relative change for a range of percentiles can be plotted
#plot_FDC_Percentile_relative_change(path_45, path_85, Catchment_Name, collect(0.1:0.1:0.9))

# --------------------------- BUDYKO -------------------------

## calculates the evaporative and aridity index
# aridity_past45, aridity_future_45, evaporative_past_45, evaporative_future_45 = aridity_evaporative_index(path_45, Area_Catchment, Catchment_Name)
# aridity_past85, aridity_future_85, evaporative_past_85, evaporative_future_85,  prec_past85, prec_future_85 = aridity_evaporative_index(path_85, Area_Catchment, Catchment_Name)
# plot_Budyko(aridity_past45, aridity_future_45, evaporative_past_45, evaporative_future_45, aridity_past85, aridity_future_85, evaporative_past_85, evaporative_future_85, path_45, Catchment_Name, nr_runs)
# plot_changes_Budyko(aridity_past45, aridity_future_45, evaporative_past_45, evaporative_future_45, aridity_past85, aridity_future_85, evaporative_past_85, evaporative_future_85, Area_Catchment, Catchment_Name, nr_runs)
#

# ----------- RUNOFF COEFFICIENT ----------------
# @time begin
#monthly_runoff_coef_past, monthly_runoff_coef_future, months, monthly_discharge_past, monthly_discharge_future = monthly_runoff_coefficient(path_45, Catchment_Name, Area_Catchment, nr_runs)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Budyko/monthly_runoff_coefficient_4.5_new.txt",hcat(monthly_runoff_coef_past, monthly_runoff_coef_future, months),',')
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/monthly_discharge_all_years_past_4.5.txt", monthly_discharge_past)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/monthly_discharge_all_years_future_4.5.txt", monthly_discharge_future)
# end
#
# @time begin
# monthly_runoff_coef_past, monthly_runoff_coef_future, months, monthly_discharge_past, monthly_discharge_future = monthly_runoff_coefficient(path_85, Catchment_Name, Area_Catchment, nr_runs)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Budyko/monthly_runoff_coefficient_8.5_new.txt",hcat(monthly_runoff_coef_past, monthly_runoff_coef_future, months),',')
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/monthly_discharge_all_years_past_8.5.txt", monthly_discharge_past)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/monthly_discharge_all_years_future_8.5.txt", monthly_discharge_future)
# end

# @time begin
#seasonal_runoff_coef_past, seasonal_runoff_coef_future, all_seasons = seasonal_runoff_coefficient(path_45, Catchment_Name, nr_runs)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Budyko/seasonal_runoff_coefficient_4.5.txt",hcat(seasonal_runoff_coef_past, seasonal_runoff_coef_future, all_seasons),',')
# end

# @time begin
# seasonal_runoff_coef_past_85, seasonal_runoff_coef_future_85, all_seasons_85 = seasonal_runoff_coefficient(path_85, Catchment_Name, nr_runs)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Budyko/seasonal_runoff_coefficient_8.5.txt",hcat(seasonal_runoff_coef_past_85, seasonal_runoff_coef_future_85, all_seasons_85 ),',')
# end
# # # --------------------------- Maximum Annual Discharge -------------------
# #
# @time begin
# average_max_Discharge_past_45, average_max_Discharge_future_45, Timing_max_Discharge_past_45, Timing_max_Discharge_future_45, All_Concentration_past_45, All_Concentration_future_45 =  change_max_Annual_Discharge(path_45, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_4.5.txt",hcat(average_max_Discharge_past_45, average_max_Discharge_future_45, Timing_max_Discharge_past_45, Timing_max_Discharge_future_45, All_Concentration_past_45, All_Concentration_future_45),',')
# end
#
# @time begin
# average_max_Discharge_past_85, average_max_Discharge_future_85, Timing_max_Discharge_past_85, Timing_max_Discharge_future_85, All_Concentration_past_85, All_Concentration_future_85 =  change_max_Annual_Discharge(path_85, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_8.5.txt",hcat(average_max_Discharge_past_85, average_max_Discharge_future_85, Timing_max_Discharge_past_85, Timing_max_Discharge_future_85, All_Concentration_past_85, All_Concentration_future_85),',')
# end
#
# annual_max_flow_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_4.5.txt",',')
# annual_max_flow_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_8.5.txt",',')
# average_max_Discharge_past_45 = annual_max_flow_45[:,1]
# average_max_Discharge_future_45 = annual_max_flow_45[:,2]
# Timing_max_Discharge_past_45 = annual_max_flow_45[:,3]
# Timing_max_Discharge_future_45 = annual_max_flow_45[:,4]
# All_Concentration_past_45 = annual_max_flow_45[:,5]
# All_Concentration_future_45 = annual_max_flow_45[:,6]
# average_max_Discharge_past_85 = annual_max_flow_85[:,1]
# average_max_Discharge_future_85 = annual_max_flow_85[:,2]
# Timing_max_Discharge_past_85 = annual_max_flow_85[:,3]
# Timing_max_Discharge_future_85 = annual_max_flow_85[:,4]
# All_Concentration_past_85 = annual_max_flow_85[:,5]
# All_Concentration_future_85 = annual_max_flow_85[:,6]
# plot_Max_Flows(average_max_Discharge_past_45, average_max_Discharge_future_45, average_max_Discharge_past_85, average_max_Discharge_future_85, Timing_max_Discharge_past_45, Timing_max_Discharge_future_45, Timing_max_Discharge_past_85, Timing_max_Discharge_future_85, Catchment_Name, Area_Catchment, nr_runs)
#
# std_timing(All_Concentration_past_45, All_Concentration_future_45, All_Concentration_past_85, All_Concentration_future_85, "Annual_Max_Discharge")
# # # # #--------------using probability distribution
# @time begin
# max_Discharge_Past_45, Max_Discharge_Future_45, Exceedance_Probability_45, Date_Past_45, Date_Future_45 =  change_max_Annual_Discharge_Prob_Distribution(path_45, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_prob_distr_4.5.txt",hcat(max_Discharge_Past_45, Max_Discharge_Future_45, Exceedance_Probability_45, Date_Past_45, Date_Future_45), ',')
# end
# @time begin
# max_Discharge_Past_85, Max_Discharge_Future_85, Exceedance_Probability_85, Date_Past_85, Date_Future_85 =  change_max_Annual_Discharge_Prob_Distribution(path_85, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_prob_distr_8.5.txt",hcat(max_Discharge_Past_85, Max_Discharge_Future_85, Exceedance_Probability_85, Date_Past_85, Date_Future_85), ',')
# end
# max_discharge_prob_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_prob_distr_4.5.txt", ',')
# max_Discharge_Past_45 = max_discharge_prob_45[:,1]
# Max_Discharge_Future_45 = max_discharge_prob_45[:,2]
# Exceedance_Probability_45 = max_discharge_prob_45[:,3]
# Date_Past_45 = max_discharge_prob_45[:,4]
# Date_Future_45 = max_discharge_prob_45[:,5]
#
# max_discharge_prob_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_prob_distr_8.5.txt", ',')
# max_Discharge_Past_85 = max_discharge_prob_85[:,1]
# Max_Discharge_Future_85 = max_discharge_prob_85[:,2]
# Exceedance_Probability_85 = max_discharge_prob_85[:,3]
# Date_Past_85 = max_discharge_prob_85[:,4]
# Date_Future_85 = max_discharge_prob_85[:,5]
#
#
# plot_change_timing_AMF_over_year(Date_Past_45, Date_Future_45, Catchment_Name, nr_runs, "45")
# plot_change_timing_AMF_over_year(Date_Past_85, Date_Future_85, Catchment_Name, nr_runs, "85")
# mean_change, min_change, max_change, std_change = plot_Max_Flows_Prob_Distribution(max_Discharge_Past_45, Max_Discharge_Future_45, max_Discharge_Past_85, Max_Discharge_Future_85, Exceedance_Probability_45, Catchment_Name)
# #
# # # ---------------- SUMMER LOW FLOWS 7 days ----------------------
# #
# @time begin
# Summer_Low_Flows_past45, Summer_Low_Flows_future45, Timing_Summer_Low_Flows_Past_45, Timing_Summer_Low_Flows_Future_45, concentration_past_45, concentration_future_45 = analyse_low_flows(path_45, Catchment_Name, "summer")
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_summer4.5.txt",hcat(Summer_Low_Flows_past45, Summer_Low_Flows_future45, Timing_Summer_Low_Flows_Past_45, Timing_Summer_Low_Flows_Future_45, concentration_past_45, concentration_future_45),',')
# end
# @time begin
# Summer_Low_Flows_past85, Summer_Low_Flows_future85, Timing_Summer_Low_Flows_Past_85, Timing_Summer_Low_Flows_Future_85, concentration_past_85, concentration_future_85  = analyse_low_flows(path_85, Catchment_Name, "summer")
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_summer8.5.txt",hcat(Summer_Low_Flows_past85, Summer_Low_Flows_future85, Timing_Summer_Low_Flows_Past_85, Timing_Summer_Low_Flows_Future_85, concentration_past_85, concentration_future_85),',')
# end
#
# low_flows_summer_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_summer8.5.txt", ',')
# Summer_Low_Flows_past85 = low_flows_summer_85[:,1]
# Summer_Low_Flows_future85 = low_flows_summer_85[:,2]
# Timing_Summer_Low_Flows_Past_85 = low_flows_summer_85[:,3]
# Timing_Summer_Low_Flows_Future_85 = low_flows_summer_85[:,4]
# concentration_past_85 = low_flows_summer_85[:,5]
# concentration_future_85 = low_flows_summer_85[:,6]
# low_flows_summer_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_summer4.5.txt", ',')
# Summer_Low_Flows_past45 = low_flows_summer_45[:,1]
# Summer_Low_Flows_future45 = low_flows_summer_45[:,2]
# Timing_Summer_Low_Flows_Past_45 = low_flows_summer_45[:,3]
# Timing_Summer_Low_Flows_Future_45 = low_flows_summer_45[:,4]
# concentration_past_45 = low_flows_summer_45[:,5]
# concentration_future_45 = low_flows_summer_45[:,6]
#
# plot_low_flows_summer(Summer_Low_Flows_past45, Summer_Low_Flows_future45, Summer_Low_Flows_past85, Summer_Low_Flows_future85,  Timing_Summer_Low_Flows_Past_45, Timing_Summer_Low_Flows_Future_45,  Timing_Summer_Low_Flows_Past_85, Timing_Summer_Low_Flows_Future_85, Catchment_Name, Area_Catchment, nr_runs)
#
# std_timing(concentration_past_45, concentration_future_45, concentration_past_85, concentration_future_85, "LowFlows/Summer")
# # # # # ----------------------- WINTER LOW FLOWS 7 DAYS ------------------------
# @time begin
# Winter_Low_Flows_past45, Winter_Low_Flows_future45, Timing_Winter_Low_Flows_Past_45, Timing_Winter_Low_Flows_Future_45, concentration_past_45, concentration_future_45  = analyse_low_flows(path_45, Catchment_Name, "winter")
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_winter4.5.txt",hcat(Winter_Low_Flows_past45, Winter_Low_Flows_future45, Timing_Winter_Low_Flows_Past_45, Timing_Winter_Low_Flows_Future_45, concentration_past_45, concentration_future_45),',')
# end
# @time begin
# Winter_Low_Flows_past85, Winter_Low_Flows_future85, Timing_Winter_Low_Flows_Past_85, Timing_Winter_Low_Flows_Future_85, concentration_past_85, concentration_future_85  = analyse_low_flows(path_85, Catchment_Name, "winter")
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_winter8.5.txt",hcat(Winter_Low_Flows_past85, Winter_Low_Flows_future85, Timing_Winter_Low_Flows_Past_85, Timing_Winter_Low_Flows_Future_85, concentration_past_85, concentration_future_85),',')
# end
#
# low_flows_winter_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_winter8.5.txt", ',')
# Winter_Low_Flows_past85 = low_flows_winter_85[:,1]
# Winter_Low_Flows_future85 = low_flows_winter_85[:,2]
# Timing_Winter_Low_Flows_Past_85 = low_flows_winter_85[:,3]
# Timing_Winter_Low_Flows_Future_85 = low_flows_winter_85[:,4]
# concentration_past_85 = low_flows_winter_85[:,5]
# concentration_future_85 = low_flows_winter_85[:,6]
# low_flows_winter_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_winter4.5.txt", ',')
# Winter_Low_Flows_past45 = low_flows_winter_45[:,1]
# Winter_Low_Flows_future45 = low_flows_winter_45[:,2]
# Timing_Winter_Low_Flows_Past_45 = low_flows_winter_45[:,3]
# Timing_Winter_Low_Flows_Future_45 = low_flows_winter_45[:,4]
# concentration_past_45 = low_flows_winter_45[:,5]
# concentration_future_45 = low_flows_winter_45[:,6]
#
# plot_low_flows_winter(Winter_Low_Flows_past45, Winter_Low_Flows_future45, Winter_Low_Flows_past85, Winter_Low_Flows_future85,  Timing_Winter_Low_Flows_Past_45, Timing_Winter_Low_Flows_Future_45,  Timing_Winter_Low_Flows_Past_85, Timing_Winter_Low_Flows_Future_85, Catchment_Name, Area_Catchment, nr_runs)
#
# std_timing(concentration_past_45, concentration_future_45, concentration_past_85, concentration_future_85, "LowFlows/Winter")
#
# # ------------------- LOW FLOWS ALL YEAR --------------------
# @time begin
# Summer_Low_Flows_past45, Summer_Low_Flows_future45, Timing_Summer_Low_Flows_Past_45, Timing_Summer_Low_Flows_Future_45, concentration_past_45, concentration_future_45 = analyse_low_flows(path_45, Catchment_Name, "none")
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_whole_year4.5.txt",hcat(Summer_Low_Flows_past45, Summer_Low_Flows_future45, Timing_Summer_Low_Flows_Past_45, Timing_Summer_Low_Flows_Future_45, concentration_past_45, concentration_future_45),',')
# end
# @time begin
# Summer_Low_Flows_past85, Summer_Low_Flows_future85, Timing_Summer_Low_Flows_Past_85, Timing_Summer_Low_Flows_Future_85, concentration_past_85, concentration_future_85  = analyse_low_flows(path_85, Catchment_Name, "none")
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_whole_year8.5.txt",hcat(Summer_Low_Flows_past85, Summer_Low_Flows_future85, Timing_Summer_Low_Flows_Past_85, Timing_Summer_Low_Flows_Future_85, concentration_past_85, concentration_future_85),',')
# end
#
# @time begin
# max_Discharge_Past_45, Max_Discharge_Future_45, Exceedance_Probability_45, Date_Past_45, Date_Future_45 =  analyse_low_flows_prob_distr(path_45, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_prob_distr_4.5.txt",hcat(max_Discharge_Past_45, Max_Discharge_Future_45, Exceedance_Probability_45, Date_Past_45, Date_Future_45), ',')
# end
# @time begin
# max_Discharge_Past_85, Max_Discharge_Future_85, Exceedance_Probability_85, Date_Past_85, Date_Future_85 =  analyse_low_flows_prob_distr(path_85, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_prob_distr_8.5.txt",hcat(max_Discharge_Past_85, Max_Discharge_Future_85, Exceedance_Probability_85, Date_Past_85, Date_Future_85), ',')
# end

## --------------- HYDROLOGICAL DROUGHT ----------------------------

#threshold = get_threshold_hydrological_drought(Discharge_Catchment, 1981,2010, 0.9, scaling_factor_discharge)
# "whole" means over whole timeseries
# season = "winter" # "summer", "winter", "none"
# @time begin
# Drought_45 = compare_hydrological_drought(path_45, threshold, season, Area_Catchment, Catchment_Name, "yearly")
# end
# @time begin
# Drought_85 = compare_hydrological_drought(path_85, threshold, season, Area_Catchment, Catchment_Name, "yearly")
# end
# plot_drought_total_deficit(Drought_45, Drought_85, threshold, Catchment_Name, Area_Catchment, season)
# plot_drought_statistics_change(Drought_45, Drought_85, threshold, Catchment_Name, Area_Catchment, season)
# plot_drought_statistics_rel_change(Drought_45, Drought_85, threshold, Catchment_Name, Area_Catchment, season)
# plot_drought_statistics(Drought_45, Drought_85, threshold, Catchment_Name, Area_Catchment, season)
# #
# season = "none" # "summer", "winter", "none"
# @time begin
# Drought_45 = compare_hydrological_drought(path_45, threshold, season, Area_Catchment, Catchment_Name, "yearly")
# end
# @time begin
# Drought_85 = compare_hydrological_drought(path_85, threshold, season, Area_Catchment, Catchment_Name, "yearly")
# end
# plot_drought_total_deficit(Drought_45, Drought_85, threshold, Catchment_Name, Area_Catchment, season)
# plot_drought_statistics_change(Drought_45, Drought_85, threshold, Catchment_Name, Area_Catchment, season)
# plot_drought_statistics_rel_change(Drought_45, Drought_85, threshold, Catchment_Name, Area_Catchment, season)
# plot_drought_statistics(Drought_45, Drought_85, threshold, Catchment_Name, Area_Catchment, season)
#
# @time begin
# Drought_Extremes_45 = compare_hydrological_drought_extremes(path_45, Area_Catchment, Catchment_Name)
# Drought_Extremes_85 = compare_hydrological_drought_extremes(path_85, Area_Catchment, Catchment_Name)
# # plot_drought_extremes_statistics(Drought_Extremes_45, Drought_Extremes_85, Catchment_Name)
# plot_drought_extremes_statistics_change(Drought_Extremes_45, Drought_Extremes_85, Catchment_Name)
# end

# @time begin
Nr_Days_Drought_monthly_past_45, Nr_Days_Drought_monthly_future_45, Deficit_monthly_past_45, Deficit_monthly_future_45, Threshold_45 = monthly_days_Q90(path_45, Area_Catchment, Catchment_Name)
Nr_Days_Drought_monthly_past_85, Nr_Days_Drought_monthly_future_85, Deficit_monthly_past_85, Deficit_monthly_future_85, Threshold_85 = monthly_days_Q90(path_85, Area_Catchment, Catchment_Name)
writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Drought/monthly_days_deficit_Q90_4.5.txt",hcat(Nr_Days_Drought_monthly_past_45, Nr_Days_Drought_monthly_future_45, Deficit_monthly_past_45, Deficit_monthly_future_45, Threshold_45),',')
writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Drought/monthly_days_deficit_Q90_8.5.txt",hcat(Nr_Days_Drought_monthly_past_85, Nr_Days_Drought_monthly_future_85, Deficit_monthly_past_85, Deficit_monthly_future_85, Threshold_85),',')
# plot_monthly_low_flows(Nr_Days_Drought_monthly_past_45, Nr_Days_Drought_monthly_future_45, Nr_Days_Drought_monthly_past_85, Nr_Days_Drought_monthly_future_85, Catchment_Name, nr_runs)
# plot_monthly_deficit(Deficit_monthly_past_45, Deficit_monthly_future_45, Deficit_monthly_past_85, Deficit_monthly_future_85, Catchment_Name, nr_runs)
# end
# # --------------- SNOW MELT ----------------
#
# @time begin
# Snow_Storage_Past85, Snow_Storage_Future85 = snow_contribution_yearly_projections(path_85, Catchment_Name)
# Snow_Storage_Past45, Snow_Storage_Future45 = snow_contribution_yearly_projections(path_45, Catchment_Name)
# end
#plot_snow_storage_contribution(Snow_Storage_Past45, Snow_Storage_Future45, Snow_Storage_Past85, Snow_Storage_Future85, Catchment_Name)
#
# @time begin
# Monthly_Discharge_Change_85, monthly_Discharge_past_85, monthly_Discharge_future_85, months_85 = change_monthly_snowmelt(path_85, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/snow_melt_months_8.5.txt",hcat(months_85, monthly_Discharge_past_85, monthly_Discharge_future_85, Monthly_Discharge_Change_85),',')
# end
# #
# @time begin
# Monthly_Discharge_Change_45, monthly_Discharge_past_45, monthly_Discharge_future_45, months_45 = change_monthly_snowmelt(path_45, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/snow_melt_months_4.5.txt",hcat(months_45, monthly_Discharge_past_45, monthly_Discharge_future_45, Monthly_Discharge_Change_45),',')
# end
#
# @time begin
# Monthly_Discharge_Change_85, monthly_Discharge_past_85, monthly_Discharge_future_85, months_85 = change_monthly_snowstorage(path_85, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/snow_storage_months_8.5.txt",hcat(months_85, monthly_Discharge_past_85, monthly_Discharge_future_85, Monthly_Discharge_Change_85),',')
# end
# #
# @time begin
# Monthly_Discharge_Change_45, monthly_Discharge_past_45, monthly_Discharge_future_45, months_45 = change_monthly_snowstorage(path_45, Catchment_Name)
# writedlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/snow_storage_months_4.5.txt",hcat(months_45, monthly_Discharge_past_45, monthly_Discharge_future_45, Monthly_Discharge_Change_45),',')
# end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
"""
    Percentile(x)

Indicate that `x` should be interpreted as a [percentile](https://en.wikipedia.org/wiki/Percentile) rather than an absolute value. For example,

- `canny(img, 1.4, (80, 20))` uses absolute thresholds on the edge magnitude image
- `canny(img, 1.4, (Percentile(80), Percentile(20)))` uses percentiles of the edge magnitude image as threshold
"""
struct Percentile{T} <: Real p::T end

"""
HomogeneousPoint(x::NTuple{N, T})

In projective geometry [homogeneous coordinates](https://en.wikipedia.org/wiki/Homogeneous_coordinates) are the
natural coordinates for describing points and lines.

For instance, the homogeneous coordinates for a planar point are a triplet of real numbers ``(u, v ,w)``, with ``w \\neq 0``.
This triple can be associated with a point ``P = (x,y)`` in Cartesian coordinates, where ``x = \\frac{u}{w}`` and ``y = \\frac{v}{w}``
[(more details)](http://www.geom.uiuc.edu/docs/reference/CRC-formulas/node6.html#SECTION01140000000000000000).

In particular, the `HomogeneousPoint((10.0,5.0,1.0))` is the standardised projective representation of the Cartesian
point `(10.0,5.0)`.
"""
struct HomogeneousPoint{T <: AbstractFloat,N}
    coords::NTuple{N, T}
end

# By overwriting Base.to_indices we can define how to index into an N-dimensional array
# given an (N+1)-dimensional [`HomogeneousPoint`](@ref) type.
# We do this by converting the homogeneous coordinates to Cartesian coordinates
# and rounding to nearest integer.
#
# For homogeneous coordinates of a planar point we return
# a tuple of permuted Cartesian coordinates, (y,x), since matrices
# are indexed  according to row and then column.
# For homogeneous coordinates of other dimensions we do not permute
# the corresponding Cartesian coordinates.
Base.to_indices(A::AbstractArray, p::Tuple{<: HomogeneousPoint}) = homogeneous_point_to_indices(p[1])

function homogeneous_point_to_indices(p::HomogeneousPoint{T,3}) where T
    if  p.coords[end] == 1
        return round(Int,  p.coords[2]), round(Int, p.coords[1])
    else
        return round(Int,  p.coords[2] / p.coords[end]), round(Int, p.coords[1] / p.coords[end])
    end
end

function homogeneous_point_to_indices(p::HomogeneousPoint)
    if  p.coords[end] == 1
        return round.(Int, p.coords)
    else
        return round.(Int, p.coords ./ p.coords[end])
    end
end

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
import Pkg.Types: VersionSpec

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.ZpN to X0Y.Z00.N00. So for example version 4.1.3p5 becomes 401.300.500.
#
# This reflects the fact that 4.2.0 will only be made if there is a change to
# Singular language itself, i.e., this is a pretty major change. The third digit
# is changed for regular interim releases, and corresponds roughly to a semver
# minor release; and the `p5` truly is a patch level.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor
# tweaks. If a rebuild of the JLL is needed which keeps the upstream version
# identical but breaks ABI compatibility for any reason, we can increment the
# minor version e.g. go from 401.300.500 to 401.301.500.
#
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 4.1.4 to 401.390.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.
#
name = "Singular"
upstream_version = v"4.3.0-2" # 4.3.0 plus some changes
version_offset = v"0.1.0"
version = VersionNumber(upstream_version.major * 100 + upstream_version.minor + version_offset.major,
                        upstream_version.patch * 100 + version_offset.minor,
                        Int(upstream_version.prerelease[1]) * 100 + version_offset.patch)

# Collection of sources required to build normaliz
sources = [
    GitSource("https://github.com/Singular/Singular.git", "a6029b833a41cc3e5095d46da10d0d366cf45c4f"),
    #ArchiveSource("https://www.mathematik.uni-kl.de/ftp/pub/Math/Singular/SOURCES/$(upstream_version.major)-$(upstream_version.minor)-$(upstream_version.patch)/singular-$(upstream_version).tar.gz",
    #              "5b0f6c036b4a6f58bf620204b004ec6ca3a5007acc8352fec55eade2fc9d63f6"),
    #DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd [Ss]ingular*

#for f in ${WORKSPACE}/srcdir/patches/*.patch; do
#    atomic_patch -p1 ${f}
#done

./autogen.sh
export CPPFLAGS="-I${prefix}/include"
./configure --prefix=$prefix --host=$target --build=${MACHTYPE} \
    --with-libparse \
    --enable-shared \
    --disable-static \
    --enable-p-procs-static \
    --disable-p-procs-dynamic \
    --enable-gfanlib \
    --with-readline=no \
    --with-gmp=$prefix \
    --with-flint=$prefix \
    --without-python \
    --with-builtinmodules=gfanlib,syzextra,customstd,interval,subsets,loctriv,gitfan,freealgebra \
    --disable-partialgb-module \
    --disable-polymake-module \
    --disable-pyobject-module \
    --disable-singmathic-module \
    --disable-systhreads-module \
    --disable-cohomo-module \
    --disable-machinelearning-module \
    --disable-sispasm-module

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
filter!(!Sys.iswindows, platforms)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpolys", :libpolys),
    LibraryProduct("libSingular", :libsingular),
    # LibraryProduct("customstd", :customstd),
    # LibraryProduct("subsets", :subsets),
    ExecutableProduct("Singular", :Singular),
    ExecutableProduct("libparse", :libparse),
    # LibraryProduct("syzextra", :syzextra),
    # LibraryProduct("interval", :interval),
    LibraryProduct("libfactory", :libfactory),
    LibraryProduct("libsingular_resources", :libsingular_resources),
    LibraryProduct("libomalloc", :libomalloc),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("cddlib_jll"),
    Dependency(PackageSpec(name="FLINT_jll"), compat = "~200.800.401"),
    Dependency("GMP_jll", v"6.2.0"),
    Dependency("MPFR_jll", v"4.1.1"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"6", julia_compat = "1.6")

max_stars_repo_name: cescalara/Yggdrasil

max_issues_repo_name: cescalara/Yggdrasil

max_forks_repo_name: cescalara/Yggdrasil


text:
#=
    Support for single variable functions i.e a <= b
=#

"""
prune_constraint!(com::CS.CoM, constraint::CS.SingleVariableConstraint, fct::MOI.ScalarAffineFunction{T}, set::LessThan{T}; logs = true) where T <: Real

Support for constraints of the form a <= b where a and b are single variables.
This function removes values which aren't possible based on this constraint.
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::CS.SingleVariableConstraint,
    fct::SAF{T},
    set::MOI.LessThan;
    logs = true,
) where {T<:Real}
    lhs = constraint.lhs
    rhs = constraint.rhs
    search_space = com.search_space
    !remove_above!(com, search_space[lhs], search_space[rhs].max) && return false
    !remove_below!(com, search_space[rhs], search_space[lhs].min) && return false
    return true
end

"""
    less_than(com::CoM, constraint::CS.SingleVariableConstraint, vidx::Int, val::Int)

Checks whether setting an `vidx` to `val` fulfills `constraint`
"""
function still_feasible(
    com::CoM,
    constraint::CS.SingleVariableConstraint,
    fct::SAF{T},
    set::MOI.LessThan,
    vidx::Int,
    val::Int,
) where {T<:Real}
    if constraint.lhs == vidx
        # if a > maximum possible value of rhs => Infeasible
        if val > com.search_space[constraint.rhs].max
            return false
        else
            return true
        end
    elseif constraint.rhs == vidx
        if val < com.search_space[constraint.lhs].min
            return false
        else
            return true
        end
    else
        error("This should not happen but if it does please open an issue with the information: SingleVariableConstraint index is neither lhs nor rhs and your model.")
    end
end

function is_constraint_solved(
    constraint::CS.SingleVariableConstraint,
    fct::SAF{T},
    set::MOI.LessThan,
    values::Vector{Int},
) where {T<:Real}
    return values[1] <= values[2]
end

"""
    is_constraint_violated(
        com::CoM,
        constraint::CS.SingleVariableConstraint,
        fct::SAF{T},
        set::Union{MOI.LessThan{T}, Union{MOI.LessThan{T},MOI.LessThan{T}}{T}},
    ) where {T<:Real}

Checks if the constraint is violated as it is currently set. This can happen inside an
inactive reified or indicator constraint.
"""
function is_constraint_violated(
    com::CoM,
    constraint::CS.SingleVariableConstraint,
    fct::SAF{T},
    set::MOI.LessThan,
) where {T<:Real}
    return com.search_space[constraint.indices[1]].min >
           com.search_space[constraint.indices[2]].max
end


text:
## Reactive signal update performance: o, h, l, c vectors from TimeArray values-component
using TimeSeries
using Reactive

nma = 50
ohlc_BA = TimeSeries.readtimearray("../../data/OHLC_BA.csv")

s_ohlc = Reactive.Input(ohlc_BA.values[1,:])
s_oc = lift(s -> (s[1], s[4]), s_ohlc, typ=(Float64, Float64))

vsig = [s_close.value - s_open.value]
function vsigupd!(vsig::Vector{Float64}, nma::Int64, nmax::Int64)
    for i in 2:nmax
        push!(s_ohlc, ohlc_BA.values[i,:])
        push!(vsig, s_oc.value[2] - s_oc.value[1])
    end
end

println("OHLC updates with 2-in-1-tuple lift-signal:")
println("reactive signal updates timing")
expr_run = parse("@time vsigupd!(vsig, nma, length(ohlc_BA))")
for i = 1:5
    eval(expr_run)
end

### results
# elapsed time: 0.21731065 seconds (23552912 bytes allocated, 17.07% gc time)
### more lift-signals really add time

Profile.clear()  # in case we have any previous profiling data
@profile vsigupd!(vsig, nma, length(ohlc_BA))
Profile.print()


text:
import OrdinaryDiffEq
import QuantumOpticsBase

"""
    liouvilleVonNeumann(rho0, tspan, H;
    \treltol=1.0e-12, abstol=1.0e-12, fout=nothing, alg=OrdinaryDiffEq.Tsit5())

Integrate Liouville-von Neumann equation to evolve states or compute propagators

``\\frac{d}{d t} \\rho(t) = - \\frac{i}{\\hbar} [ \\hat H, \\rho(t) ], \\quad \\hbar = 1``.

# Arguments
* `rho0`: Initial state vector (can be a bra or a ket) or initial propagator.
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `H`: Arbitrary operator specifying the Hamiltonian.
* `reltol`: Relative tolerance for OrdinaryDiffEq solver and its inner states.
* `abstol`: Absolute tolerance for OrdinaryDiffEq solver and its inner states.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The state `rho` is neither
        normalized nor permanent! It is still in use by the ode solver and
        therefore must not be changed.
* `alg`: Algorithm with which OrdinaryDiffEq will solve LvN equation.
"""
function LvN_sI(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::Float64 = 1.0e-12,
    abstol::Float64 = 1.0e-12,
    alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.Tsit5(),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    p = (agg.core, agg.tools, agg.operators, W0, eltype(W0))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    dLvN_(t, rho, drho, p) = dLvN_sI(t, rho, drho, p)

    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate(
        tspan_,
        dLvN_,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dLvN_sI(
    t::AbstractFloat,
    rho::T,
    drho::T,
    p
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, elementtype = p
    Ham = aggOperators.Ham

    Ham_0 = aggOperators.Ham_0
    Ham_I = aggOperators.Ham_I
    Ham_II_t = getInteractionHamIPicture(Ham_0, Ham_I, t)

    U_t = evolutionOperator(Ham, t)
    W_t = U_t * W0 * U_t'
    U_0_op = evolutionOperator(Ham_0, t)
    W_int_t = U_0_op' * W_t * U_0_op

    K = -elementtype(im) * (Ham_II_t.data * W_int_t.data - W_int_t.data * Ham_II_t.data)
    drho.data[:, :] = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
    return drho
end

function LvN_sS(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::Float64 = 1.0e-12,
    abstol::Float64 = 1.0e-12,
    alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.Tsit5(),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    p = (agg.core, agg.tools, agg.operators, W0, eltype(W0))
    rho0 = trace_bath(W0, agg.core, agg.tools; vib_basis=agg.operators.vib_basis)
    dLvN_(t, rho, drho, p) = dLvN_sS(t, rho, drho, p)

    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    OpenQuantumSystems.integrate(
        tspan_,
        dLvN_,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dLvN_sS(
    t::AbstractFloat,
    rho::T,
    drho::T,
    p
) where {B<:Basis,T<:Operator{B,B}}
    aggCore, aggTools, aggOperators, W0, elementtype = p
    Ham = aggOperators.Ham

    U_t = evolutionOperator(Ham, t)
    W_t = U_t * W0 * U_t'

    K = -elementtype(im) * (Ham.data * W_t.data - W_t.data * Ham.data)
    drho.data[:, :] = trace_bath(K, aggCore, aggTools; vib_basis=aggOperators.vib_basis)
    return drho
end

function LvN_SI(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::Float64 = 1.0e-12,
    abstol::Float64 = 1.0e-12,
    alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.Tsit5(),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    p = (agg.operators, eltype(W0))
    dLvN_(t, W::T, dW::T, p) = dLvN_SI(t, W, dW, p)

    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = W0.data
    state = T(W0.basis_l, W0.basis_r, W0.data)
    dstate = T(W0.basis_l, W0.basis_r, W0.data)

    OpenQuantumSystems.integrate(
        tspan_,
        dLvN_,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dLvN_SI(
    t::AbstractFloat,
    W::T,
    dW::T,
    p
) where {B<:Basis,T<:Operator{B,B}}
    aggOperators, elementtype = p

    Ham_0 = aggOperators.Ham_0.data
    Ham_I = aggOperators.Ham_I.data
    Ham_II_t = getInteractionHamIPicture(Ham_0, Ham_I, t)

    dW.data[:, :] = -elementtype(im) * (Ham_II_t * W.data - W.data * Ham_II_t)
    return dW
end

function LvN_SS(
    W0::T,
    tspan::Array,
    agg::Aggregate;
    reltol::Float64 = 1.0e-12,
    abstol::Float64 = 1.0e-12,
    alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.Tsit5(),
    fout::Union{Function,Nothing} = nothing,
    kwargs...,
) where {B<:Basis,T<:Operator{B,B}}
    p = (agg.operators, eltype(W0))
    dLvN_(t, W::T, dW::T, p) = dLvN_SS(t, W, dW, p)

    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = W0.data
    state = T(W0.basis_l, W0.basis_r, W0.data)
    dstate = T(W0.basis_l, W0.basis_r, W0.data)

    OpenQuantumSystems.integrate(
        tspan_,
        dLvN_,
        x0,
        state,
        dstate,
        fout;
        p = p,
        reltol = reltol,
        abstol = abstol,
        alg = alg,
        kwargs...,
    )
end

function dLvN_SS(
    t::AbstractFloat,
    W::T,
    dW::T,
    p
) where {B<:Basis,T<:Operator{B,B}}
    aggOperators, elementtype = p
    Ham = aggOperators.Ham.data
    dW.data[:, :] = -elementtype(im) * (Ham * W.data - W.data * Ham)
    return dW
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
"""
    rosenbrock!(x, g; b=2)

Return the value of Rosenbrock function in `N` dimensions.
"""
function rosenbrock!(x::AbstractVector, g::AbstractVector; b=2)
    f = zero(eltype(g))
    fill!(g, 0)
    inds = eachindex(x, g)
    for i in 2:last(inds)
        f += (1 - x[i-1])^2 + b * (x[i] - x[i-1]^2)^2
        g[i-1] += 2 * (x[i-1] - 1) + 4 * b * x[i-1] * (x[i-1]^2 - x[i])
        g[i] += 2 * b * (x[i] - x[i-1]^2)
    end
    return f, g
end

@testset "Convergence on the Rosenbrock function" begin
    dim = 2
    x0 = fill(-2, dim)
    @testset for method in OPT_TYPES
        opt = method(x0)
        optresult = optimize!(opt, rosenbrock!, x0; maxiter=1000)
        @test isapprox(optresult.argument, [1, 1], rtol=0.05)
    end
end


text:
using Serialization
include("CFRcs.jl")

function samp(probdict)
    s = rand()
    cumm = 0.0
    for (k, v) in probdict
        cumm += v
        if s <= cumm
            return k
        end
    end
end

function showpolicy(profile, hand, hist)
    hand = map(n -> UInt8(n), hand)
    hist = map(n -> UInt8(n), hist)
    k = key(hand, hist)
    policy = profile.policys[k].strategysum |> normalize
    for (k, v) in policy
        if v >= 0.01
            @printf("%2.0f => %0.2f\n", k, v)
        end
    end
end

function move(dict, hand, rawhist)
    hist = copy(rawhist)
    k = key(hand, hist)
    if !haskey(dict, k)
        while length(hist) > 0 && hist[1] >= 0 && hist[1] < 20
            popfirst!(hist)
        end
    end
    k = key(hand, hist)
    if !haskey(dict, k) && length(hist) > maxhist
        k = hist[end-maxhist+1] != 0 ? key(hand, hist[end-maxhist+1:end]) :
            key(hand, hist[end-maxhist+2:end])
    end
    d = dict[k].strategysum |> normalize
    samp(d)
end

function selfplay(dict, hands)
    hist = UInt8[]
    player = 1
    while true
        mv = move(dict, hands[player, :], hist)
        push!(hist, mv)
        if mv == 1
            w = win(hist, hands)
            x = (player == 1 && w) || (player == 2 && !w) ? 1 : -1
            return player, x, hist
        end
        player = 3 - player # change the player
    end
end

function headsup(dict1, dict2, hands)
    hist = UInt8[]
    player = rand(1:2)
    while true
        mv = player == 1 ? move(dict1, hands[1, :], hist) : move(dict2, hands[2, :], hist)
        push!(hist, mv)
        if mv == 1
            w = win(hist, hands)
            x = (player == 1 && w) || (player == 2 && !w) ? 1 : -1
            return player, x, hist
        end
        player = 3 - player
    end
end

function randplay(profile, dict1, dict2, n)
    total = 0
    hs = Matrix{UInt8}(undef, (2, quant))
    for _ = 1:n
        (_, x, _) = headsup(dict1, dict2, randhands!(profile, hs))
        total += x
    end
    println(total / n, " ", total, " ", n)
end

function play(profile)
    hist = UInt8[]
    hands = Matrix{UInt8}(undef, (2, profile.cards_per_hand))
    randhands!(profile, hands)
    you = rand(1:2)
    her = 3 - you
    player = 1
    println("Your hand: ", hands[you, :])
    while true
        if player == you
            print("Your move: ")
            t = readline()
            mv = tryparse(UInt8, t)
            if mv ∉ actions(profile.all_bids, hist)
                println("XXX - Illegal action, please try again - XXX")
                continue
            end
            player = her
        else
            mv = move(profile.policys, hands[player, :], hist)
            player = you
            if mv == 0
                println("Opponent : Challenge")
            elseif mv == 1
                println("Opponent : Counts")
            else
                @printf("Opponent : %.0f\n", mv)
            end
        end
        push!(hist, mv)
        if mv == 1
            w = win(hist, hands)
            x = (player == you && w) || (player == her && !w) ? -1 : 1
            println()
            println(x == 1 ? "++++++++++++++++ WIN +++++++++++++++++" :
                    "---------------- LOSE ----------------")
            println()
            println("Opponents hand: ", hands[her, :])
            return x
        end
    end
end

function games(profile)
    n = 0
    score = 0
    println("----------------------------")
    @printf("Liars Poker: %.0f Card, %.0f Ranks\n", profile.cards_per_hand, profile.num_ranks)
    println("0 - Challenge, 1 - Count")
    println("----------------------------")
    println()
    while true
        n += 1
        score += play(profile)
        @printf("Total: %.0f, Games: %.0f, Average: %.2f\n", score, n, score / n)
        println("======================================")
        println()
        @printf("Hand %.0f\n", n + 1)
    end
end

max_stars_repo_name: jeffreyrosenbluth/LP

max_issues_repo_name: jeffreyrosenbluth/LP

max_forks_repo_name: jeffreyrosenbluth/LP


text:
export Transformer, does_pcb_fit, core_name, ferrite_name

"""
    Transformer(ferrite, windings)

**Fields**
- `ferrite`   -- `FerriteProperties` object
- `windings`  -- `Windings` object
"""
struct Transformer
  ferrite :: FerriteProperties
  windings :: Windings
end

effective_volume(t::Transformer) = effective_volume(windings(t))
effective_area(t::Transformer) = effective_area(windings(t))
effective_length(t::Transformer) = effective_length(windings(t))
winding_aperature_height(t::Transformer) = winding_aperature_height(windings(t))
core_name(t::Transformer) = core_name(windings(t))
ferrite_name(t::Transformer) = name(ferrite(t))
"""
    does_pcb_fit(t)

`true` if PCB thickness is <= winding aperature height.
"""
does_pcb_fit(t) = pcb_thickness(t) <= winding_aperature_height(t)


text:
struct Zhuyin{T <: AbstractString}
    contents::T
end

struct Pinyin{T <: AbstractString}
    contents::T
end


text:
module Includeonce
# help include :
# include([mapexpr::Function,] path::AbstractString)
# help?> Base.include
#   Base.include([mapexpr::Function,] [m::Module,] path::AbstractString)

includeonce(path)=include(path)
end # module


text:
function DEFAULT_INITIALIZER end

"""
$(SIGNATURES)

Builds the callback to check positivity of the system density matrix during the evolution.
"""
function PositivityCheckCallback()
    affect! = function (u, t, integrator)
        if ndims(u) == 1
            dim = convert(Int, sqrt(length(u)))
            ρ = reshape(u, dim, dim)
        else
            ρ = u
        end
        if !check_positivity(ρ)
            @warn "The density matrix becomes negative at time $t."
            terminate!(integrator)
        end
        u_modified!(integrator, false)
    end
    FunctionCallingCallback(affect!, func_everystep=true, func_start=false)
end

"""
$(SIGNATURES)

Builds the callback to apply instantaneous pulses during the evolution.

...
# Arguments
- `tstops`: time points to apply pulses.
- `pulse_update`: udpate function for the state of the system. It taks two argument `pulse_update(c, i)` where `c` is the state vector or density matrix and `i` is the index of the pulse being applied.
...
"""
function InstPulseCallback(tstops, pulse_update)
    state = Ref{Int}(1)
    affect! = function (integrator)
        for c in full_cache(integrator)
            pulse_update(c, state[])
        end
        state[] += 1
    end

    initialize = function (c, u, t, integrator)
        state[] = 1
        u_modified!(integrator, false)
    end
    PresetTimeCallback(tstops, affect!, initialize=initialize)
end

function FluctuatorCallback(F::OpenQuantumBase.FluctuatorLiouvillian, initialize)
    time_choice = function (integrator)
        next_t = integrator.t + F.next_τ
        if next_t > integrator.sol.prob.tspan[2]
            return nothing
        else
            return next_t
        end
    end

    affect! = function (integrator)
        F.n = sum(F.b0, dims=1)[:]
        OpenQuantumBase.next_state!(F)
        u_modified!(integrator, false)
    end

    initialize_fluc = function (cb, u, t, integrator)
        OpenQuantumBase.reset!(F, initialize)
        u_modified!(integrator, false)
    end

    IterativeCallback(
        time_choice,
        affect!,
        initialize=initialize_fluc,
        save_positions=(false, false),
    )
end

function LindbladJumpCallback()
    r = Ref{Float64}(rand(Float64))
    condition = function (u, t, integrator)
        real(u' * u) - r[]
    end

    affect! = function (integrator)
        A = OpenQuantumBase.lindblad_jump(
            integrator.p.L,
            integrator.u,
            integrator.p,
            integrator.t,
        )
        r[] = rand()
        new_state = normalize(A * integrator.u)
        for c in full_cache(integrator)
            c .= new_state
        end
    end

    initialize = function (c, u, t, integrator)
        r[] = rand()
        u_modified!(integrator, false)
    end

    ContinuousCallback(condition, affect!, save_positions=(true, true), initialize=initialize)
end

"""
$(SIGNATURES)

Create callback to apply manifold retraction after each step of the ODE solver.
"""
function ManifoldRetractionCallback()
    affect! = function (integrator)
        rmul!(integrator.u, 1 / norm(integrator.u))
    end
    condition = (u, t, integrator) -> true
    DiscreteCallback(condition, affect!, save_positions=(false, save))
end

max_stars_repo_name: huangjingzheng/OpenQuantumTools

max_issues_repo_name: huangjingzheng/OpenQuantumTools

max_forks_repo_name: huangjingzheng/OpenQuantumTools


text:
"""
An `ArrayOperator` combines an `AbstractArray` with a source and destination
dictionary.

From the 'AbstractArray' we expect that `mul!`, `size`, `copy`, `diag`, `isdiag`, `inv`, `getindex`,
`adjoint`, `conj` are implemented where possible.
"""
abstract type ArrayOperator{T} <: DictionaryOperator{T} end

"Verify that size(A) == (length(dest),length(src)) ?"
verify_size(src, dest, A) = size(A) == (length(dest),length(src))

"Suggest a suitable src dictionary for a given matrix."
Asrc(A::AbstractMatrix{T}) where {T} = DiscreteVectorDictionary{T}(size(A,2))
"Suggest a suitable dest dictionary for a given matrix."
Adest(A::AbstractMatrix{T}) where {T} = DiscreteVectorDictionary{T}(size(A,1))

getindex(op::ArrayOperator, i::Int, j::Int) = op.A[i,j]

# Delegate to object
for f in (:size, :isefficient)
    @eval $f(op::ArrayOperator) = $f(op.A)
end#object related features

for f in (:isdiag, :diag, :eigvals, :Matrix, :svdvals, :norm, :rank, :sparse)
    @eval $f(op::ArrayOperator) = $f(unsafe_matrix(op))
end#matrix related features

for f in (:inv, :conj, :adjoint, :sqrt)
    @eval $f(op::ArrayOperator) = ArrayOperator($f(op.A), dest(op), src(op))
end#operator related features

unsafe_wrap_operator(src, dest, op::ArrayOperator) = similar_operator(op, src, dest)

similar_operator(op::ArrayOperator, src::Dictionary, dest::Dictionary) =
    ArrayOperator(op.A, src, dest)

apply_inplace!(op::ArrayOperator, coef_srcdest::AbstractVector) = _apply_inplace!(op, op.A, coef_srcdest)
_apply_inplace!(op::ArrayOperator, A::AbstractArray, x) = mul!(x, A, x)

apply!(op::ArrayOperator, coef_dest::AbstractVector, coef_src::AbstractVector) = _apply!(op, op.A, coef_dest, coef_src)
mul!(op::ArrayOperator, coef_dest::AbstractVector, coef_src::AbstractVector) = mul!(coef_dest, op.A, coef_src)

# Be forgiving for matrices: if the coefficients are multi-dimensional, reshape to a linear array first.
apply!(op::ArrayOperator{T}, coef_dest::AbstractArray{T,N1}, coef_src::AbstractArray{T,N2}) where {T,N1,N2} =
    apply!(op, reshape(coef_dest, length(coef_dest)), reshape(coef_src, length(coef_src)))

apply_inplace!(op::ArrayOperator{T}, coef_srcdest::AbstractArray{T,N}) where {T,N} =
    apply_inplace!(op, reshape(coef_srcdest, length(coef_srcdest)))

_apply!(op::ArrayOperator, A::AbstractArray, coef_dest, coef_src) = mul!(coef_dest, A, coef_src)

pinv(op::ArrayOperator, tolerance::Real = eps(real(eltype(op)))) = ArrayOperator(pinv(op.A, tolerance), src(op), dest(op))

matrix(op::ArrayOperator) = copy(op.A)
unsafe_matrix(op::ArrayOperator) = op.A

isefficient(::AbstractArray) = false
isefficient(::AbstractSparseArray) = true
isefficient(::UniformScaling) = true
isefficient(::ToeplitzMatrices.AbstractToeplitz) = true
isefficient(::Diagonal) = true

string(op::ArrayOperator) = string(op, op.A)
string(op::ArrayOperator,array) = "Multiplication by "*string(typeof(op.A))

"A `MatrixOperator` is a wrapper around an abstract matrix type."
struct MatrixOperator{T,ARRAY} <: ArrayOperator{T}
    A       ::  ARRAY
    src     ::  Dictionary
    dest    ::  Dictionary

    function MatrixOperator{T,ARRAY}(A::ARRAY, src::Dictionary, dest::Dictionary) where {T,ARRAY}
        @assert verify_size(src, dest, A)
        new(A, src, dest)
    end
end

"Find a suitable element type for a dictionary operator with the given construct arguments."
deduce_eltype(A::AbstractArray{T}, args...) where {T} = T
deduce_eltype(A::AbstractArray, src::Dictionary, args...) = promote_type(eltype(A), operatoreltype(src))
deduce_eltype(A::AbstractArray, src::Dictionary, dest::Dictionary, args...) = promote_type(eltype(A), operatoreltype(src,dest))
deduce_eltype(src::Dictionary, args...) = operatoreltype(src)
deduce_eltype(src::Dictionary, A::AbstractArray, args...) = promote_type(eltype(A),operatoreltype(src))
deduce_eltype(src::Dictionary, z::Number, args...) = promote_type(typeoef(z),operatoreltype(src))
deduce_eltype(src::Dictionary, dest::Dictionary, args...) = operatoreltype(src, dest)
deduce_eltype(src::Dictionary, dest::Dictionary, A::AbstractArray, args...) = promote_type(eltype(A),operatoreltype(src, dest))
deduce_eltype(src::Dictionary, dest::Dictionary, z::Number, args...) = promote_type(typeoef(z),operatoreltype(src, dest))

MatrixOperator(args...) = MatrixOperator{deduce_eltype(args...)}(args...)

MatrixOperator{T}(A::AbstractArray) where {T} = MatrixOperator{T}(A, Asrc(A))
MatrixOperator{T}(A::AbstractArray, src::Dictionary) where {T} = MatrixOperator{T}(A, src, src)

MatrixOperator{T}(A::AbstractArray{T}, src::Dictionary, dest::Dictionary) where {T} =
    MatrixOperator{T,typeof(A)}(A, src, dest)
MatrixOperator{T}(A::AbstractArray{S}, src::Dictionary, dest::Dictionary) where {S,T} =
    MatrixOperator{T}(convert(AbstractArray{T}, A), src, dest)

## Banded operators

"""
A banded operator of which every row contains equal elements.

The top row starts at index offset, the second row at step+offset.
"""
const HorizontalBandedOperator{T} = MatrixOperator{T,HorizontalBandedMatrix{T}}

HorizontalBandedOperator(args...) = HorizontalBandedOperator{deduce_eltype(args...)}(args...)

HorizontalBandedOperator{T}(src::Dictionary, dest::Dictionary, array::AbstractVector, step::Int=1, offset::Int=0) where {T} =
    HorizontalBandedOperator{T}(HorizontalBandedMatrix(length(dest), length(src), T.(array), step, offset), src, dest)

"""
A banded operator of which every column contains equal elements.

The top column starts at index offset, the second column at step+offset.
"""
const VerticalBandedOperator{T} = MatrixOperator{T,VerticalBandedMatrix{T}}

VerticalBandedOperator(args...) = VerticalBandedOperator{deduce_eltype(args...)}(args...)

VerticalBandedOperator{T}(src::Dictionary, dest::Dictionary, array::AbstractVector, step::Int=1, offset::Int=0) where {T} =
    VerticalBandedOperator{T}(VerticalBandedMatrix(length(dest), length(src), T.(array), step, offset), src, dest)

## Index extension and restriction

"An IndexRestriction selects a subset of coefficients based on their indices."
const IndexRestriction{T,N,I} = MatrixOperator{T,RestrictionIndexMatrix{T,N,I}}

IndexRestriction(src::Dictionary, args...) = IndexRestriction{operatoreltype(src)}(src, args...)
IndexRestriction(src::Dictionary, dest::Dictionary, args...) = IndexRestriction{operatoreltype(src,dest)}(src, dest, args...)

IndexRestriction{T}(A::RestrictionIndexMatrix{T,N,I}, src::Dictionary, dest::Dictionary) where {T,N,I} =
    IndexRestriction{T,N,I}(A, src, dest)

IndexRestriction{T}(src::Dictionary, subindices::AbstractVector) where {T} =
    IndexRestriction{T}(src, src[subindices], subindices)

function IndexRestriction{T}(src::Dictionary, dest::Dictionary, subindices::AbstractVector) where {T}
    @assert length(dest)==length(subindices) && length(dest)<=length(src)
    IndexRestriction{T}(RestrictionIndexMatrix{T}(size(src), subindices), src, dest)
end

subindices(op::IndexRestriction) = subindices(op.A)

show(io::IO, mime::MIME"text/plain", op::IndexRestriction) = composite_show(io, mime, op)
Display.object_parentheses(op::IndexRestriction) = false
Display.stencil_parentheses(op::IndexRestriction) = false
Display.displaystencil(op::IndexRestriction) = ["Restrict[𝕀 → ", subindices(op.A), "]"]

"An `IndexExtension` embeds coefficients in a larger set based on their indices."
const IndexExtension{T,N,I} = MatrixOperator{T,ExtensionIndexMatrix{T,N,I}}

IndexExtension(src::Dictionary, args...) = IndexExtension{operatoreltype(src)}(src, args...)
IndexExtension(src::Dictionary, dest::Dictionary, args...) = IndexExtension{operatoreltype(src,dest)}(src, dest, args...)

IndexExtension{T}(A::ExtensionIndexMatrix{T,N,I}, src::Dictionary, dest::Dictionary) where {T,N,I} =
    IndexExtension{T,N,I}(A, src, dest)
IndexExtension{T}(dest::Dictionary, subindices::AbstractVector) where {T} =
    IndexExtension{T}(dest[subindices], dest, subindices)
function IndexExtension{T}(src::Dictionary, dest::Dictionary, subindices::AbstractVector) where {T}
    @assert length(src)==length(subindices) && length(dest)>=length(src)
    IndexExtension{T}(ExtensionIndexMatrix{T}(size(dest), subindices), src, dest)
end

subindices(op::IndexExtension) = subindices(op.A)

show(io::IO, mime::MIME"text/plain", op::IndexExtension) = composite_show(io, mime, op)
Display.object_parentheses(op::IndexExtension) = false
Display.stencil_parentheses(op::IndexExtension) = false
Display.displaystencil(op::IndexExtension) = ["Extend[", subindices(op.A), " → 𝕀]"]

# Diagonal operator and friends

"A diagonal operator."
const DiagonalOperator{T,D} = MatrixOperator{T,Diagonal{T,D}}

# Convert various arguments to a diag matrix with a concrete diag vector
# with elements of type T.
to_diag(A::AbstractArray{T}) where {T} = to_diag(T, A)
to_diag(::Type{T}, A::AbstractVector) where {T} = to_diag(T, collect(A))
to_diag(::Type{T}, A::Vector{T}) where {T} = Diagonal(A)
to_diag(::Type{T}, A::FillArrays.AbstractFill{T,1}) where {T} = Diagonal(A)
to_diag(::Type{T}, A::Vector{S}) where {S,T} = to_diag(T, convert(Vector{T}, A))
to_diag(::Type{T}, A::Diagonal{T,Array{T,1}})  where {T} = A
to_diag(::Type{T}, A::Diagonal{T}) where {T} = to_diag(T, diag(A))
to_diag(::Type{T}, A::Diagonal{S}) where {S,T} = to_diag(T, convert(Diagonal{T}, A))

DiagonalOperator(args...) = DiagonalOperator{deduce_eltype(args...)}(args...)

DiagonalOperator{T}(A::AbstractMatrix) where {T} = DiagonalOperator{T}(A, Asrc(A), Adest(A))
DiagonalOperator{T}(A::AbstractMatrix, src::Dictionary) where {T} = DiagonalOperator{T}(A, src, src)

DiagonalOperator{T}(A::Diagonal{T,D}, src::Dictionary, dest::Dictionary) where {T,D} =
    DiagonalOperator{T,D}(A, src, dest)
DiagonalOperator{T}(A::AbstractArray, src::Dictionary, dest::Dictionary) where {T} =
    DiagonalOperator{T}(to_diag(T,A), src, dest)

# Support other order of arguments
DiagonalOperator{T}(src::Dictionary, A::AbstractArray) where {T} =
    DiagonalOperator{T}(A, src, src)
DiagonalOperator{T}(src::Dictionary, dest::Dictionary, A::AbstractArray) where {T} =
    DiagonalOperator{T}(A, src, dest)

DiagonalOperator(src::Dictionary, dest::Dictionary, A::AbstractOuterProductArray) =
    tensorproduct(map(DiagonalOperator, components(src), components(dest), components(A))...)

isinplace(op::DiagonalOperator) = true

_apply_inplace!(op::ArrayOperator, A::Diagonal, x) = mul!(x, A, x)
_apply!(op::ArrayOperator, A::Diagonal, coef_dest, coef_src) = mul!(coef_dest, A, coef_src)

# show(io::IO, ::MIME"text/plain", op::DiagonalOperator) = print(io, "DiagonalOperator with element type $(eltype(op))")
Display.displaystencil(op::DiagonalOperator) = ["DiagonalOperator(", diag(op.A), ')']

"The identity operator"
const IdentityOperator{T} = DiagonalOperator{T,<:Ones{T}}

IdentityOperator(src::Dictionary) = IdentityOperator{operatoreltype(src)}(src)
IdentityOperator(src::Dictionary, dest::Dictionary) = IdentityOperator{operatoreltype(src,dest)}(src, dest)

function IdentityOperator{T}(src::Dictionary, dest = src) where {T}
    diag = Ones{T}(length(src))
    DiagonalOperator{T,typeof(diag)}(Diagonal(diag), src, dest)
end

isidentity(op::IdentityOperator) = true

broadcasted(::Base.Broadcast.DefaultArrayStyle{N}, ::Type{T}, a::Ones{S,N}) where {T<:Number,S,N} =
    Ones{T}(axes(a))
broadcasted(::Base.Broadcast.DefaultArrayStyle{N}, ::Type{T}, a::Zeros{S,N}) where {T<:Number,S,N} =
    Zeros{T}(axes(a))

strings(op::IdentityOperator) = ("Identity operator of size $(size(op)) (T=$(eltype(op)))", strings(src(op)))
symbol(op::IdentityOperator) = "I"

"Scaling by a scalar value"
const ScalingOperator{T} = DiagonalOperator{T,<:Fill{T}}

ScalingOperator(args...) = ScalingOperator{deduce_eltype(args...)}(args...)
ScalingOperator(z::Number, src::Dictionary, dest = src) =
    ScalingOperator{promote_type(typeof(z),operatoreltype(src,dest))}(z, src, dest)
ScalingOperator(z::UniformScaling{S}, src::Dictionary, dest = src) where {S} =
    ScalingOperator{promote_type(S,operatoreltype(src,dest))}(z, src, dest)

# Allow other order of arguments
ScalingOperator(src::Dictionary, a) = ScalingOperator(a, src)
ScalingOperator(src::Dictionary, dest::Dictionary, a) = ScalingOperator(a, src, dest)
ScalingOperator{T}(src::Dictionary, a) where {T} = ScalingOperator{T}(a, src)
ScalingOperator{T}(src::Dictionary, dest::Dictionary, a) where {T} = ScalingOperator{T}(a, src, dest)

function ScalingOperator{T}(z::Number, src::Dictionary, dest::Dictionary = src) where {T}
    diag = Fill(T(z), length(src))
    DiagonalOperator{T}(diag, src, dest)
end
function ScalingOperator{T}(z::UniformScaling, src::Dictionary, dest::Dictionary = src) where {T}
    diag = Fill(T(z.λ), length(src))
    DiagonalOperator{T}(diag, src, dest)
end

scalar(op::ScalingOperator) = op.A.diag.value
strings(op::ScalingOperator) = ("Scaling by $(scalar(op))",)
symbol(op::ScalingOperator) = "α"

"The zero operator maps everything to zero."
const ZeroOperator{T} = MatrixOperator{T,<:Zeros{T,2}}

ZeroOperator(src::Dictionary) = ZeroOperator{operatoreltype(src)}(src)
ZeroOperator(src::Dictionary, dest::Dictionary) = ZeroOperator{operatoreltype(src,dest)}(src,dest)

function ZeroOperator{T}(src::Dictionary, dest::Dictionary = src) where {T}
    A = Zeros{T}(length(dest), length(src))
    MatrixOperator{T}(A, src, dest)
end

strings(op::ZeroOperator) = ("Zero operator of size $(size(op)) (T=$(eltype(op)))", strings(src(op)))
symbol(op::ZeroOperator) = "0"

"A circulant operator."
const CirculantOperator{T,S} = MatrixOperator{T,Circulant{T,S}}

CirculantOperator(args...) = CirculantOperator{deduce_eltype(args...)}(args...)

to_circulant(v::AbstractVector) = to_circulant(eltype(v), v)
to_circulant(::Type{T}, v::AbstractVector{S}) where {S,T} = to_circulant(T, T.(v))
to_circulant(::Type{T}, v::AbstractVector{T}) where {T} = Circulant(v)
to_circulant(::Type{T}, A::Circulant{T}) where {T} = A

CirculantOperator{T}(v::AbstractVector) where {T} = CirculantOperator{T}(to_circulant(T, v))
CirculantOperator{T}(A::AbstractMatrix) where {T} = CirculantOperator{T}(A, Asrc(A), Adest(A))
CirculantOperator{T}(A::AbstractArray, src::Dictionary) where {T} =
    CirculantOperator{T}(A, src, src)
CirculantOperator{T}(A::AbstractArray, src::Dictionary, dest::Dictionary) where {T} =
    CirculantOperator{T}(to_circulant(T, A), src, dest)
CirculantOperator{T}(A::Circulant{T,S}, src::Dictionary, dest::Dictionary) where {T,S} =
    CirculantOperator{T,S}(A, src, dest)

# Temporary fix for ToeplitzMatrices issue #47
inv(C::CirculantOperator{T}) where {T} = CirculantOperator{T}(Circulant(ifft(1 ./ C.A.vcvr_dft)), dest(C), src(C))
inv(C::CirculantOperator{T}) where {T<:Real} = CirculantOperator{T}(Circulant(real(ifft(1 ./ C.A.vcvr_dft))), dest(C), src(C))

# Support other order of arguments
CirculantOperator{T}(src::Dictionary, A::AbstractArray) where {T} =
    CirculantOperator{T}(A, src)
CirculantOperator{T}(src::Dictionary, dest::Dictionary, A::AbstractArray) where {T} =
    CirculantOperator{T}(A, src, dest)

# Convert any dictionary operator to a circulant operator
function CirculantOperator(op::DictionaryOperator)
    e = zeros(src(op))
    e[1] = 1
    C = CirculantOperator(op*e, src(op), dest(op), )
    # Check accuracy of the conversion:
    e = rand(src(op))
    @assert C*e≈op*e
    C
end

"A dense matrix operator"
const DenseMatrixOperator{T} = MatrixOperator{T,Array{T,2}}

DenseMatrixOperator(args...) = DenseMatrixOperator{deduce_eltype(args...)}(args...)

DenseMatrixOperator{T}(A::AbstractMatrix) where {T} = DenseMatrixOperator{T}(A, Asrc(A), Adest(A))
DenseMatrixOperator{T}(A::AbstractMatrix, src::Dictionary) where {T} = DenseMatrixOperator{T}(A, src, src)

ArrayOperator(A::AbstractArray) =
    ArrayOperator(A, Asrc(A), Adest(A))
ArrayOperator(A::AbstractArray, src::Dictionary) =
    ArrayOperator(A, src, src)
ArrayOperator(A::UniformScaling, src::Dictionary) =
    ArrayOperator(A, src, src)

ArrayOperator(args...) = ArrayOperator{deduce_eltype(args...)}(args...)

ArrayOperator{T}(A::Matrix{T}, src::Dictionary, dest::Dictionary) where {T} =
    DenseMatrixOperator{T}(A, src, dest)
ArrayOperator{T}(A::Diagonal, src::Dictionary, dest::Dictionary) where {T} =
    DiagonalOperator{T}(A, src, dest)
ArrayOperator{T}(A::ExtensionIndexMatrix, src::Dictionary, dest::Dictionary) where {T} =
    IndexExtension{T}(A, src, dest)
ArrayOperator{T}(A::RestrictionIndexMatrix, src::Dictionary, dest::Dictionary) where {T} =
    IndexRestriction{T}(A, src, dest)
ArrayOperator{T}(A::HorizontalBandedMatrix{T}, src::Dictionary, dest::Dictionary) where {T} =
    HorizontalBandedOperator{T}(A, src, dest)
ArrayOperator{T}(A::VerticalBandedMatrix{T}, src::Dictionary, dest::Dictionary) where {T} =
    VerticalBandedOperator{T}(A, src, dest)
ArrayOperator{T}(A::Circulant, src::Dictionary, dest::Dictionary) where {T} =
    CirculantOperator{T}(A, src, dest)
ArrayOperator{T}(A::Zeros, src::Dictionary, dest::Dictionary) where {T} =
    ZeroOperator{T}(src, dest)

ArrayOperator{T}(A::AbstractMatrix, src::Dictionary, dest::Dictionary) where {T} =
    MatrixOperator{T}(A, src, dest)

export SparseMatrixOperator
"A sparse matrix operator"
const SparseMatrixOperator{T} = MatrixOperator{T,SparseMatrixCSC{T,Int}}
function SparseMatrixOperator(op; options...)
    A = sparse(op; options...)
    MatrixOperator(A, src(op), dest(op))
end

sparse(A::SparseMatrixOperator) = A.A


text:
module MultifidelityTrajectories
# this directory
const topdirectory = normpath(joinpath(@__DIR__, ".."))

# import AircraftSystems
include(joinpath(topdirectory, "..", "AircraftSystems", "src", "AircraftSystems.jl"))
const AS = AircraftSystems

# import Altro
const Dates = AS.Dates
# import Dates
# const RD = Altro.RobotDynamics
const ST = AS.StaticArrays
# const TO = Altro.TrajectoryOptimization
# using LinearAlgebra, StaticArrays
const LA = AS.LinearAlgebra
const LS = AS.LaTeXStrings
const PP = AS.PyPlot
import SNOW

# set constants
const TODAY = replace(string(Dates.today()),"-" => "")

# include files
include("CRC3_geometry.jl")
include("CRC3_system.jl")
include("coordinate_frames.jl")
include("CRC3_dynamics.jl")
include("CRC3_aerodynamics.jl")
include("utilities.jl")
include("integrate_dynamics.jl")
include("trim.jl")
include("optimization.jl")
include("plots.jl")

end # module

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
function fib(n)
    a = 0
    b = 1
    while n > 0
        b = b + a
        a = b - a
        n = n - 1
    end
    a
end

for n = 0:10
    println(fib(n))
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
include("solver_preamble.jl")
import ProxSDP
factory = with_optimizer(ProxSDP.Optimizer, log_verbose=false)
config = MOI.Test.Config(atol=2e-2, rtol=2e-2)
@testset "Linear" begin
    Tests.linear_test(factory, config, [
        "dsos_options_pricing",
        # Expression: termination_status(model) == MOI.OPTIMAL
        #  Evaluated: MathOptInterface.INFEASIBLE_OR_UNBOUNDED == MathOptInterface.OPTIMAL
        "dsos_concave_then_convex_cubic",
        # Expression: JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        #  Evaluated: MathOptInterface.INFEASIBLE_POINT == MathOptInterface.FEASIBLE_POINT
        "dsos_quartic_comparison",
        "dsos_horn"
    ])
end
@testset "SOC" begin
    Tests.soc_test(factory, config, [
        "sdsos_options_pricing",
        # Expression: termination_status(model) == MOI.OPTIMAL
        #  Evaluated: MathOptInterface.INFEASIBLE_OR_UNBOUNDED == MathOptInterface.OPTIMAL
        "sdsos_concave_then_convex_cubic",
        # ConstraintDual not implemented for VoV-in-SOC
        "sdsos_univariate_quadratic",
        "sdsos_horn", "sdsos_bivariate_quadratic"
    ])
end
@testset "SDP" begin
    Tests.sd_test(factory, config, [
        "sos_options_pricing",
        # ProxSDP does not return duals for nonlinear Conic constraints. Only linear constraints (equalities and inequalities) can be queried.
        "quartic_ideal", "quartic_ideal_4", "quartic_ideal_4_rem",
        # with γ=3.9 it should be infeasible: Test Failed at /home/blegat/.julia/dev/SumOfSquares/test/Tests/maxcut.jl:37
        # Expression: JuMP.dual_status(model) == MOI.INFEASIBILITY_CERTIFICATE
        #  Evaluated: MathOptInterface.INFEASIBLE_POINT == MathOptInterface.INFEASIBILITY_CERTIFICATE
        "maxcut",
        # ITERATION_LIMIT
        "choi", "motzkin",
        # ArgumentError: ModelLike of type ProxSDP.Optimizer does not support accessing the attribute MathOptInterface.ConstraintDual(1)
        "sos_univariate_quadratic",
        "sos_horn",
        "sos_bivariate_quadratic",
        "quadratic_infeasible_lyapunov_switched_system",
        "quartic_infeasible_lyapunov_switched_system",
        "sosdemo5_infeasible"
    ])
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Documenter, NetCDF

makedocs(
    modules = [NetCDF],
    format = Documenter.HTML(),
    sitename = "NetCDF.jl",
    authors = "Fabian Gans and contributors",
    pages = [
        "Home" => "index.md",
        "Manual" => ["quickstart.md", "highlevel.md", "intermediate.md", "strings.md"],
    ],
)

deploydocs(repo = "github.com/JuliaGeo/NetCDF.jl.git")


text:
using LogDensity
using Test

@testset "LogDensity.jl" begin
    r = [ parse(Float64, fl) for fl in readlines("boundestdata.csv") ]
    ( logf, β ) = logdensity(r, collect(0.0:0.1:0.2), 0.5; anal=false )
    @test β ≈ [ -0.409345755364715,-0.321324897846474,-0.401367070582909 ]
    @test logf ≈ [ -0.739190491046256,-0.78965235304968,-0.818590378954698 ]
    β2 = logdensityderivatives(r, collect(0.0:0.1:0.2), 0.5; anal=false )
    @test β ≈ β2

    ( logf, β ) = logdensity(r, collect(0.0:0.1:0.2), 0.5; S=2, anal=false )
    @test β ≈ [0.19052091185219694 -2.434548533489367; -0.3401926625357224 0.09748925928955005; -0.2893470695141677 -0.7973089783812821]
    @test logf ≈ [-0.7904810573794027, -0.7890805128340765, -0.8132211046368243]
    β2 = logdensityderivatives(r, collect(0.0:0.1:0.2), 0.5; S=2, anal=false )
    @test β ≈ β2

    ( logf, β ) = logdensity(r, collect(0.0:0.1:0.2), 0.5; S=3, mz=LogDensity.quartic,anal=false )
    @test β ≈ [-1.572088272730215 15.39013993699334 -71.87012956581209; -0.18776905560267698 -2.786004255503017 14.675398492375955; -0.28141402629073176 -0.10470410790118856 -4.789506265957136]
    @test logf ≈ [-0.7185403492935303, -0.7818410094136379, -0.8190838536404783]
    β2 = logdensityderivatives(r, collect(0.0:0.1:0.2), 0.5; S=3, anal=false )
    @test β ≈ β2

end


text:
using ProgressMeter
using TimerOutputs
using ParticleInCell

include("run_function.jl")

t, energy = run( 1 ) # trigger building
@show nstep = 1000
@time t, energy = run( nstep, npm = 500 )
show(to)

open("results_jl.dat", "w") do f

    for i in 1:nstep
        println(f, t[i], " ", energy[i])
    end

end

println()


text:
__precompile__()
module Remove_eq

export detect_eq_kurtosis, stalta, remove_eq

using SeisIO, DSP, StatsBase, ORCA, SeisIO, Printf, PlotlyJS

"""
    detect_eq_kurtosis(data::SeisChannel,tw::Float64=60.0, threshold::Float64=3.0, overlap::Float64=30)

find earthquake by kurtosis threshold

# Input:
    - `data::SeisChannel`    : SeisData from SeisIO
    - `tw::Float64`  : time window to evaluate if earthquake is contained.
    - `threshold::Float64` : kurtosis threshold: if kurtosis > threshold, the time window contains earthquake
    - `overlap::Float64`           : overlap of time window to control the margin of earthquake removal. (large overlap assigns large margin before and after earthquake.)

    kurtosis evaluation following Baillard et al.(2013)
"""
function detect_eq_kurtosis(data::SeisChannel; tw::Float64=60.0, kurtosis_threshold::Float64=3.0, overlap::Float64=30)

    #convert window lengths from seconds to samples
    twsize = trunc(Int, tw * data.fs)
    overlapsize = trunc(Int, overlap * data.fs)

    #calculate how much to move beginning of window each iteration
    slide = twsize-overlapsize

    #kurtosis of timeseries
    ku1 = data.misc["kurtosis"]#[:]
    #reset long window counter and triggers for current channel
    i = 0

    #loop through current channel by sliding
    while i < length(ku1) - twsize

        #check if last window and change long window length if so
        if length(ku1) - i < twsize
            twsize = length(ku1)-i
        end

        #define chunk of data based on long window length and calculate long-term average
        twtrace = @views ku1[i+1:i+twsize]

        if !isnothing(findfirst(x -> x > kurtosis_threshold, twtrace))
            #this time window includes earthquake
            for tt= i+1:i+twsize
                data.misc["eqtimewindow"][tt] = false
            end
        end

        #advance long window
        i = i + slide

    end

    return data

end

"""
    detect_eq_stalta(data::SeisChannel,tw::Float64=60.0, threshold::Float64=3.0, overlap::Float64=30)

find earthquake and tremors by STA/LTA

# Input:
    - `data::SeisChannel`  : SeisChannel from SeisIO
    - `longwin::Float64`   : long time window
    - `shortwin::Float64`  : short time window
    - `threshold::Float64` : STA/LTA threshold: if STA/LTA > threshold, the time window contains earthquake
    - `overlap::Float64`   : overlap of time window to control the margin of earthquake removal. (large overlap assigns large margin before and after earthquake.)
    - `stalta_absoluteclip::Float64`   : clip if maximum absolute value exceeds this number (for the purpose of removing incoherent noise)

    original code written by Seth Olinger. For our purpose, overlap is applied for short time window
"""
function detect_eq_stalta(data::SeisChannel,longWinLength::Float64, shortWinLength::Float64, threshold::Float64, overlap::Float64, stalta_absoluteclip::Float64;
                            datagap_eps::Float64=1e-8)

    #convert window lengths from seconds to samples
    longWin = trunc(Int,longWinLength * data.fs)
    shortWin = trunc(Int,shortWinLength * data.fs)
    overlapWin = trunc(Int,overlap * data.fs)
    trace = @view data.x[:]

    #calculate how much to move beginning of window each iteration
    slide = longWin

    #save weight
    eqweight = deepcopy(data.misc["eqtimewindow"])
    #print("before manipulation"); println(count(eqweight))

    # manipulate weight to avoid underestimation of lta due to data gap
    # threshold zero signal by 0.001% of lta within the timeseries, which encompasses small computational error such as 1e-21
    ltaall = StatsBase.mean(abs.(trace))
    for i = 1:length(trace)
        if isless(abs(trace[i]), ltaall*datagap_eps)
            # this is regarded as zero signal (data gap)
            eqweight[i] = false
        end
    end

    #print("after manipulation"); println(count(eqweight))

    #define long window counter and empty trigger vector
    i = 1

    #loop through current channel by sliding
    while i < length(trace)

        #check if last window and change long window length if so
        if length(trace) - i < longWin
            longWin = length(trace)-i
        end

        #define chunk of data based on long window length and calculate long-term average
        longTrace = trace[i:i+longWin]

        lta = StatsBase.mean(abs.(longTrace), weights(eqweight[i:i+longWin]))

        #reset short window counter
        n = 0
        breakflag = false
        #loop through long window in short windows
        while n <= longWin - shortWin

            #define chunk of data based on short window length and calculate short-term average
            shortTrace = @view trace[i+n:i+n+shortWin]
            shortTraceWeight = @view eqweight[i+n:i+n+shortWin]
            #sta = mean(abs.(shortTrace))
            sta = StatsBase.mean(abs.(shortTrace), weights(shortTraceWeight))
            #sta = StatsBase.mean(abs.(shortTrace))
            stamax = maximum(abs.(shortTrace))
            #calculate sta/lta ration
            staLta = sta/lta

            #record detection time if sta/lta ratio exceeds threshold
            if staLta > threshold || stamax > stalta_absoluteclip
                #triggers[i+n] = 1
                #this time window includes earthquake
                for tt= i+n:i+n+shortWin
                    data.misc["eqtimewindow"][tt] = false
                end
            end

            if breakflag
                break;
            end

            #advance short window
            n = n + shortWin - overlapWin

            #adjust for the last edge of long window
            if n >= longWin - shortWin
                n = longWin - shortWin
                breakflag = true
            end

        end

        #advance long window
        i = i + slide

    end

    return data

end

"""
    remove_eq(data::SeisChannel)

remove earthquake by kurtosis and STA/LTA threshold

# Input:
    - `data::SeisData`    : SeisData from SeisIO

"""
function remove_eq(data::SeisChannel, data_origin::SeisChannel, plot_kurtosis_α::Float64, max_taper_dur::Float64,
    plot_boxheight::Float64, plot_span::Int64, plot_fmt::String, fodir::String, tstamp::String, tvec::Array{Float64,1}, IsSaveFig::Bool)

    eqidlist = data.misc["eqtimewindow"][:]
    nx = length(data.x)

    i = 1

    t1 = []
    t2 = []
    y1 = []
    y2 = []

    tt1 = 0
    tt2 = 0
    tt3 = 0
    tt4 = 0
    tt5 = 0

    while i <= nx
        if !eqidlist[i]
            push!(t1, tvec[i])

            t1id = i

            #find next id
            tt1 += @elapsed nexttrueid = findfirst(x -> x == true, eqidlist[t1id:end])

            if isnothing(nexttrueid)
                # all data is removed
                t2id = nx
                iinc = nx
            else
                t2id = (t1id - 1) + nexttrueid - 1 # first false id to end false id
                iinc = t2id - i + 1 # first true as end of false + 1
            end

            push!(t2, tvec[t2id])

            max_wintaper_duration = Int(data.fs*max_taper_dur)

            # compute α given maximum taper length
            eq_length = t2id-t1id+1
            taper_length = 2*max_wintaper_duration
            tukey_length = eq_length + taper_length
            invert_tukey_α = taper_length/tukey_length

            # define inverted tukey window
            tt2 += @elapsed invtukeywin = -tukey(Int(tukey_length), invert_tukey_α) .+ 1

            # slice tukey window if it exceeds array bounds
            tt3 += @elapsed if t1id < max_wintaper_duration + 1
                left_overflow = (max_wintaper_duration-t1id)+1
                invtukeywin = @views invtukeywin[left_overflow+1:end]
                # leftmost t
                left = 1
            else
                # full taper length
                left = t1id - max_wintaper_duration
            end

            tt4 += @elapsed if (nx - t2id + 1) < max_wintaper_duration + 1
                right_overflow = (max_wintaper_duration-(nx - t2id + 1))+1
                invtukeywin = @views invtukeywin[1:end-right_overflow]
                # rightmost t
                right = nx
            else
                # full taper length
                right = t2id + max_wintaper_duration
            end

            # apply tukey window
            tt5 += @elapsed data.x[left:right] .*= invtukeywin

            #boxsize
            push!(y1, -plot_boxheight)
            push!(y2, plot_boxheight)

        else
            iinc = 1
        end

        i += iinc

    end

    #println([tt1, tt2, tt3, tt4, tt5])

    if IsSaveFig

        normalized_amp = 0.5 * maximum(data_origin.x[1:plot_span:end])

        trace1 = scatter(;x=tvec[1:plot_span:end], y=data_origin.x[1:plot_span:end]./ normalized_amp,
         mode="lines", name="raw data", line_color="black")

        trace2 = scatter(;x=tvec[1:plot_span:end], y=data.x[1:plot_span:end]./ normalized_amp,
         mode="lines", name="after remove", line_color="blue")

        #to plot kurtosis computed points
        kurtid = findall(x -> !iszero(data.misc["kurtosis"][x]), 1:length(data.misc["kurtosis"]))

        if isempty(filter(!isnan, data.misc["kurtosis"][kurtid]))
            kurtnormalize = 1.0
        else
            kurtnormalize = maximum(filter(!isnan, data.misc["kurtosis"][kurtid]))
        end

        trace3 = scatter(;x=tvec[kurtid], y=data.misc["kurtosis"][kurtid] ./ kurtnormalize .* plot_kurtosis_α,
         line_color="red", mode="lines", name="kurtosis")

        if !isempty(t1)
         shapes = PlotlyJS.rect(t1, t2, y1, y2; fillcolor="#ff99cc", opacity=0.3, line_width=0)
         layout = Layout(shapes=shapes, width=1200, height=600,
             xaxis=attr(title="Time [hour]"),
             yaxis=attr(title="Normalized velocity"),
             font =attr(size=13),
             showlegend=true,
             title = @sprintf("%s %s", data.id, tstamp))

             p = plot([trace1; trace2; trace3],layout)

        else
         layout = Layout(width=1200, height=600,
             xaxis=attr(title="Time [hour]"),
             yaxis=attr(title="Normalized velocity"),
             font =attr(size=12),
             showlegend=true,
             title = @sprintf("%s %s", data.id, tstamp))

         p = plot([trace1; trace2; trace3],layout)

        end

        figdir = joinpath(fodir, "fig")
        mkpath(figdir)
        figname = @sprintf("%s/%s_%s.%s", figdir, data.id, tstamp, plot_fmt)
        savefig(p, figname)
        #display(p)
        #println("press return for next plot...")
        #readline()
    end

    return data

end

end


text:
using Postprocessing
using Test
using LinearAlgebra, Random

@testset "total variation denoising" begin
  println(devnull, reference(total_variation_denoising))

  for T in (Float32, Float64)
    x = range(zero(T), 2*one(T), length=1000)
    y_true = sinpi.(x)
    y_noisy = y_true .+ 0.1 .* randn(T, length(x))
    y_denoised = total_variation_denoising(y_noisy, 1)
    @test norm(y_denoised - y_true) < norm(y_noisy - y_true)
    @test_throws ArgumentError total_variation_denoising(y_noisy, -1)
    @test_throws DimensionMismatch total_variation_denoising!(similar(y_noisy, T, length(y_noisy)-1), y_noisy, 1)
  end
end

@testset "group sparse total variation denoising" begin
  println(devnull, reference(group_sparse_total_variation_denoising))

  for T in (Float32, Float64)
    x = range(zero(T), 2*one(T), length=1000)
    y_true = sinpi.(x)
    y_noisy = y_true .+ 0.1 .* randn(T, length(x))
    y_denoised = group_sparse_total_variation_denoising(y_noisy, 1, group_size=3)
    @test norm(y_denoised - y_true) < norm(y_noisy - y_true)
    @test total_variation(y_denoised) < total_variation(y_noisy)
    @test_throws ArgumentError group_sparse_total_variation_denoising(y_noisy, -1, group_size=3)
    @test_throws DimensionMismatch group_sparse_total_variation_denoising!(similar(y_noisy, T, length(y_noisy)-1), y_noisy, 1)
  end
end


text:
using JuMP, EAGO, Gurobi

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

# ----- Variables ----- #
x_Idx = Any[1, 2]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[1], 0.0)
JuMP.set_lower_bound(x[2], 0.0)
JuMP.set_upper_bound(x[1], 3.0)
JuMP.set_upper_bound(x[2], 4.0)

# ----- Constraints ----- #
@NLconstraint(m, e2, 8* (x[1])^3-2* (x[1])^4-8* (x[1])^2+x[2] <= 2.0)
@NLconstraint(m, e3, 32* (x[1])^3-4* (x[1])^4-88* (x[1])^2+96*x[1]+x[2] <= 36.0)

# ----- Objective ----- #
@objective(m, Min, -x[1] -x[2])

JuMP.optimize!(m)
run_time = backend(m).optimizer.model.optimizer._run_time
objval = objective_value(m)
println("run time: $run_time")


text:
@testset "Test getreducedmargin d = 1" begin
    # d = 1

    idx = reshape([0],(1,1))

    reducedmargin = getreducedmargin(idx)

    @test reducedmargin == reshape([1],(1,1))

    idx = reshape([0 ,1 , 2 ,3],(4,1))

    reducedmargin = getreducedmargin(idx)

    @test reducedmargin == reshape([4],(1,1))

    # Reorder variable
    for i=1:1000
        p = randperm(size(idx,1))

        idxshuffle = idx[p,:]
        reducedmargin = getreducedmargin(idx)
        @test reducedmargin == reshape([4],(1,1))
    end
end

@testset "Test getreducedmargin d = 2" begin
    # d = 2

    idx = [0 0]
    reducedmargin = getreducedmargin(idx)
    @test reducedmargin == [0 1; 1 0]

    idx = [0 0; 0 1; 0 2; 0 3; 1 0; 1 1; 2 0; 3 0]
    reducedmargin = getreducedmargin(idx)
    @test reducedmargin == [0 4; 1 2; 2 1; 4 0]

    # Reorder variable
    for i=1:1000
        p = randperm(size(idx,1))

        idxshuffle = idx[p,:]
        reducedmargin = getreducedmargin(idx)
        @test reducedmargin == [0 4; 1 2; 2 1; 4 0]
    end
end

@testset "Test getreducedmargin d = 3" begin
    # d = 3

    idx = [0 0 0]

    reducedmargin = getreducedmargin(idx)

    @test reducedmargin == [0 0 1; 0 1 0; 1 0 0]

    idx = [0 0 0; 0 0 1; 0 1 0; 0 1 1; 1 1 0; 1 0 1]

    reducedmargin = getreducedmargin(idx)

    @test reducedmargin ==
     [0     0     2;
      0     2     0;
      1     0     0;
      1     1     1]

    # Reorder variable
    for i=1:1000
        p = randperm(size(idx,1))

        idxshuffle = idx[p,:]
        reducedmargin = getreducedmargin(idx)
        @test reducedmargin ==
         [0     0     2;
          0     2     0;
          1     0     0;
          1     1     1]
    end

    # Difficult test
    idx = [0 0 0; 0 0 1; 0 1 0; 1 0 0; 0 0 2; 0 0 3; 1 1 0; 1 2 0; 1 1 0]

    reducedmargin = getreducedmargin(idx)

    @test reducedmargin == [ 0 0 4; 0 1 1; 0 2 0; 1 0 1; 2 0 0]

end

@testset "Test getreducedmargin d = 4" begin
    # d = 4

    idx = [0 0 0 0]

    reducedmargin = getreducedmargin(idx)

    @test reducedmargin == [0 0 0 1; 0 0 1 0;0 1 0 0; 1 0 0 0]
end

@testset "Test update reduced margin d = 2" begin

    lowerset = [0 0; 0 1; 0 2; 0 3; 1 0; 1 1; 2 0; 3 0]
    reduced_margin = getreducedmargin(lowerset)
    reduced_margin0 = deepcopy(reduced_margin)
    lowerset0 = deepcopy(lowerset)

    for idx = 1:size(reduced_margin,1)
        # @show idx
        new_lowerset, new_reduced_margin = updatereducedmargin(lowerset, reduced_margin, idx)
        # @show reduced_margin[idx,:]
        # @show new_reduced_margin
        # @show new_lowerset
        @test lowerset == lowerset0
        @test reduced_margin == reduced_margin0

        @test new_lowerset[1:end-1,:] == lowerset0

        @test new_lowerset[end,:] == reduced_margin[idx,:]

        @test all([any(x in eachslice(getreducedmargin(new_lowerset);dims = 1)) for x in eachslice(new_reduced_margin; dims = 1)])
    end
end

@testset "Test update reduced margin d = 3" begin

    lowerset = [0 0 0; 0 0 1; 0 1 0; 0 2 0; 1 0 0; 0 0 2; 0 0 3; 1 1 0; 1 2 0; 1 1 0]
    reduced_margin = getreducedmargin(lowerset)
    reduced_margin0 = deepcopy(reduced_margin)
    lowerset0 = deepcopy(lowerset)

    for idx = 1:size(reduced_margin,1)
        # @show idx
        new_lowerset, new_reduced_margin = updatereducedmargin(lowerset, reduced_margin, idx)
        # @show reduced_margin[idx,:]
        # @show new_reduced_margin
        # @show new_lowerset
        @test lowerset == lowerset0
        @test reduced_margin == reduced_margin0

        @test new_lowerset[1:end-1,:] == lowerset0

        @test new_lowerset[end,:] == reduced_margin[idx,:]

        @test all([any(x in eachslice(getreducedmargin(new_lowerset);dims = 1)) for x in eachslice(new_reduced_margin; dims = 1)])
    end
end


text:
using Test
import DefaultApplication

TRAVIS = get(ENV, "TRAVIS", "false") == "true"

function isrunning(program; details = false)
    processes = read(`ps axco command`, String)
    if occursin(program, processes)
        true
    else
        if details
            @info("process $(program) not found in processes")
            for p in sort(split(processes, "\n"))
                @info p
            end
        end
        false
    end
end

if Sys.islinux()
    if TRAVIS
        # Test that a text file is opened by jmacs (setup in .travis.yml).
        @test !isrunning("jmacs") # check that it is not running accidentally
        testfile = "/tmp/test.txt"
        write(testfile, "test text")
        ## uncomment lines below for debug information
        # @info("environment",
        #       XDGMIMETYPE = chomp(read(`xdg-mime query filetype $(testfile)`, String)),
        #       XDGMIMEDEFAULT = chomp(read(`xdg-mime query default text/plain`, String)))
        @info "opening $(testfile)"
        DefaultApplication.open(testfile)
        sleep(1)
        @info "jmacs should now be running"
        @test isrunning("jmacs"; details = true)
    else
        @warn "Tests are only ran in Travis VM."
    end
end


text:
using Distributions
using Plots
using QuadGK

Σ₁₁, Σ₁₂, Σ₁₃ = 1.0, 0.5, 0.1
Σ₂₁, Σ₂₂, Σ₂₃ = 0.5, 2.0, 0.2
Σ₃₁, Σ₃₂, Σ₃₃ = 0.1, 0.2, 0.5
μ₁=0
μ₂=-1.0
μ₃=2.0

μ = [
    μ₁
    μ₂
    μ₃
]

Σ = [
    Σ₁₁ Σ₁₂ Σ₁₃
    Σ₂₁ Σ₂₂ Σ₂₃
    Σ₃₁ Σ₃₂ Σ₃₃
]

μ_marginal = [
    μ₁
    μ₂
]

Σ_marginal = [
    Σ₁₁ Σ₁₂
    Σ₂₁ Σ₂₂
]

d = MvNormal(μ, Σ)
data = rand(d,250)

"""
p(x,y) = \\int p(x,y,z) dz
"""
function marginal(d, x, y)
    integral, err = quadgk(
        z->pdf(d,[x,y,z]),
        μ[2] - maximum(abs.(Σ)),
        μ[2] + maximum(abs.(Σ)),
        rtol = 1e-8,
    )
    return integral
end

marginal(y,z)=marginal(d,y,z)
contour(collect(-5:0.1:5), collect(-5:0.1:5), marginal)
scatter!(data[1,:], data[2,:], label="numeric")

# the pdf function p(x,y) should be 𝒩(μ_marginal, Σ_marginal).
d = MvNormal(μ_marginal, Σ_marginal)
data = rand(d,250)
scatter!(data[1,:],data[2,:],label="theoretial")

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# This file was generated, do not modify it. # hide
@load XGBoostRegressor

xgb_mdl = XGBoostRegressor(num_round=10, max_depth=10)
xgb = machine(xgb_mdl, X, y)
fit!(xgb, rows=train)

ypred = predict(xgb, rows=test)
round(rms(ypred, y[test]), sigdigits=3)


text: abstract type ModelSelectionResult end


text:
#!/usr/bin/env julia

export pe25

function pe25(n::Integer=1000)
    fib(n::Integer) = ([1 1;1 0]^n)[1,2]
    i=big(1)
    while ndigits(fib(i))<n
        i+=1
    end
    i
end

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe25()|>print
end

max_stars_repo_path: 25. 1000-digit Fibonacci number/25.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 25. 1000-digit Fibonacci number/25.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 25. 1000-digit Fibonacci number/25.jl

max_forks_repo_name: miRoox/ProjectEuler