text:
import FLOWUnsteady
uns = FLOWUnsteady
vlm = uns.vlm
import GeometricTools
gt = GeometricTools
using PyPlot
using LinearAlgebra: norm, dot, cross
using LaTeXStrings                      #Need for setting up plots
using Dates

FLOWUnsteady.vlm.VLMSolver._mute_warning(true)
#Mute the warnings

# velocities = [];       # Total Velocity over wing

today = replace(string(Dates.today()),"-" => "")
savefolder = "data"

disp_plot = true

function main(;disp_plot = true)

    #Define aspects of wing. Single, symmetric wing.

    #Forward Swept Wing
    span = 1.0;              #wing span
    aspectratio = 10.0;      #wing aspect ratio
    taperratio = 0.5;        #wing taper ratio
    wingtwist = 4.0;         #wing twist
    wingsweep = -5.0;        #wing sweep in degrees
    wingdihedral = 6.0;      #wing dihedral in degrees
    savepath = "PlotCollective";

    #Glider
    # span = 1.0;
    # aspectratio = 15.0;
    # taperratio = 0.9;
    # wingtwist = 1.0;
    # wingsweep = 1.0;
    # wingdihedral = 1.0;

    #This bulids the wing using the variables we gave it.
    mainwing = vlm.simpleWing(span,aspectratio,taperratio,wingtwist,wingsweep,wingdihedral);

    #Now we create a wing system
    system = vlm.WingSystem();       #creates an empty system

    vlm.addwing(system, "mainwing", mainwing);  #Adds wing to wystem

    #We will now save the system to a .vtk file. MAKE SURE YOU ARE IN A NON-IMPORANT LOCATION TO SAVE THIS.
    #I added an empty folder for this.

    #Prep file
    Vinf(X,t) = [1,0,0];        #non-dimensional function defining free stream velocity
    vlm.setVinf(system, Vinf)   #set freestream velocity for the system
    "Freestream Velocity set"
    #"$Vinf"

    run_name = "tutorial"           #define identifier at beginning of file names
    save_path = "./simplewing/"     #define directory where files will be saved

    run(`rm -rf $save_path`)        #clear out directory where files will be saved
    run(`mkdir $save_path`)         #re-create directory fresh

    vlm.save(system, run_name; path=save_path)  #save geometry in a .vtk file format

    #now we will show the wing using paraview

    # run(`paraview --data="$(save_path)/$(run_name)_mainwing_vlm.vtk"`)

    #We now have the wing set up and running in ParaView.

    #We will now add a rotor
    rotor_file = "apc10x7.csv"          # Rotor geometry
    data_path = uns.def_data_path       # Path to rotor database

    #this creates a rotor, and generates a heck-a-ton of output, which we suppress.
    R, B = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]] #get the radius for later

    Vinf(X,t) = [1,0,0];
    # Vinf = [1,0,0];

    rotor = uns.generate_rotor(rotor_file; pitch=0.0,
                                            n=10, CW=true, ReD=1.5e6,
                                            verbose=true, xfoil=false,
                                            data_path=data_path,
                                            plot_disc=false);

    #Create a rotor object bsed on parameters above.
    rotors = vlm.Rotor[rotor];
    #We now set the location of the vehicle to be in line with the rotor.

    #First we set the vehicle at the origin.
    vehicleorigin = [0.0; 0.0; 0.0];
    vehicleaxis = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0];

    #Now we set the rotor origin in line with the vehicle
    rotororigin = [-0.1; 0.0; 0.0];
    RPM_Val = 1000

    #And now line the rest of the rotor up, part by part.
    for rotor in rotors
        vlm.setcoordsystem(rotor, rotororigin, vehicleaxis; user=true);
        vlm.setRPM(rotor, RPM_Val)
    end

    #Now we will store the rotor in a tuple
    rotor_systems = (rotors,);
    # vlm.setRPM(rotors, RPM);       #Maybe this.

    #finally, we add the rotor to the wing system. Again, part by part.
    for rotor in rotors
        vlm.addwing(system, run_name, rotor);
    end

    #Now we run it to see the visualization
    run(`rm -rf $save_path`)
    run(`mkdir $save_path`)
    #
    vlm.setVinf(system, Vinf)

    # vlm.setRPM(system, RPM);

    vlm.save(system, run_name; path=save_path)
    # run(`paraview --data="$(save_path)/tutorial_mainwing_vlm.vtk;tutorial_tutorial_Blade1_vlm.vtk;tutorial_tutorial_Blade2_vlm.vtk;tutorial_tutorial_Blade1_loft.vtk;tutorial_tutorial_Blade2_loft.vtk;"`)

    #create a vlm system to add our main wing to.
    vlm_system = vlm.WingSystem();
    vlm.addwing(vlm_system, "mainwing", mainwing);

    #Create a wake system. Add VLM system and rotor to this wake system.
    wake_system = vlm.WingSystem();
    # vlm.addwing(wake_system, "SolveVlm", vlm_system)

    # for rotor in rotors
    #     vlm.addwing(wake_system, run_name, rotor);
    # end

    #We need to define tilting objects (as they are part of the code) but ours will be empty

    tilting_systems = ();   #Empty tuple

    #The vehicle, rotor, and freestream are all set up now.

    #We will now define a maneuver.

    #Start with Vehicle velocity function
    Vvehicle(t) = [-1.0,0.0,0.0];

    #Define Vehicle angle
    anglevehicle(t) = zeros(3); #Flat angle

    #define angle of tilting system. We have no tilting system, so make it empty
    angle = ();

    #Now, lastly, we define the rotation rate function for the rotors
    RPM_fun(t) = 10.0;   #1.0 is constant rotation
    RPM = (RPM_fun,);

    #using each individual element, we can create a maneuver now.
    maneuver = uns.KinematicManeuver(angle, RPM, Vvehicle, anglevehicle);

    #We can plot the maneuver now, which is good.
    uns.plot_maneuver(maneuver);

    #This is cool, but lets put it in paraview.
    vlm.setVinf(system, Vinf);
    vehicle = uns.QVLMVehicle( system;
                              tilting_systems = tilting_systems,
                              rotor_systems = rotor_systems,
                              vlm_system = vlm_system,
                              wake_system = wake_system,
                            );
    #Now we need some reference information: reference velocity, total time, and initial conditions
    Forces_wing = [];      # Total force over wing.
    forces_per_span = [];  # Force per unit span on wing.
    cL_wing = [];          # Coefficient of lift of wing.
    L_wing = [];           # Lift over wing.
    cD_wing = [];          # Coefficint of drag over wing
    D_wing = [];           # Effective drag over wing.
    y2b_vec = [];
    CdCD_vec = [];
    Gamma_vec = [];
    V_vec = [];
    W_vec = [];
    dt = 0;
    Time = 0;
    nsteps = 300;

    Flag = false;           #testing for the forces output

    Vref = 10.0         #define a reference velocity for the vehicle
    ttot = 1.0          #define a total simulation time, in seconds
    nsteps = 300        #define the number of steps the simulation will take
    prev_wing = deepcopy(mainwing)
    rhoinf = 1.23
    b = span
    qinf = 0.5*Vref^2 * rhoinf
    ar = aspectratio
    wake_coupled = false

    function monitor(sim, PFIELD, T, DT; figname="monitor_plot", nsteps_plot=1, disp_plot = disp_plot, figsize_factor = 1)

            if PFIELD.nt!=0 && PFIELD.nt%nsteps_plot==0 && disp_plot
                # figure(figname)

                # Force at each VLM element
                Ftot = uns.calc_aerodynamicforce(mainwing, prev_wing, PFIELD, Vinf, DT,
                                                                rhoinf; t=PFIELD.t)
                # Add Forces to forces vector
                Fx = [v[1] for v in Ftot]
                Fy = [v[2] for v in Ftot]
                Fz = [v[3] for v in Ftot]
                L, D, S = uns.decompose(Ftot, [0,0,1], [-1,0,0]);
                vlm._addsolution(mainwing, "L", L)
                vlm._addsolution(mainwing, "D", D)
                vlm._addsolution(mainwing, "S", S)
                resultant_Force = [sum(Fx),sum(Fy),sum(Fz)]
                push!(Forces_wing, resultant_Force);

                # #Add Lift to lift vector
                push!(L_wing, L);

                # Force per unit span at each VLM element
                ftot = uns.calc_aerodynamicforce(mainwing, prev_wing, PFIELD, Vinf, DT,
                                            rhoinf; t=PFIELD.t, per_unit_span=true);
                fx = [v[1] for v in ftot]
                fy = [v[2] for v in ftot]
                fz = [v[3] for v in ftot]
                l, d, s = uns.decompose(ftot, [0,0,1], [-1,0,0]);
                resultant_force = fx + fy + fz;
                push!(forces_per_span, resultant_force);

                # Lift of the wing
                Lwing = norm(sum(L))
                CLwing = Lwing/(qinf*b^2/ar)
                ClCL = norm.(l) / (Lwing/b)

                # Add cl to vector
                push!(cL_wing, ClCL);

                # Drag of the wing
                Dwing = norm(sum(D))
                CDwing = Dwing/(qinf*b^2/ar)
                CdCD = [sign(dot(this_d, [1,0,0])) for this_d in d].*norm.(d) / (Dwing/b) # Preserves the sign of drag

                # Add Dwing and CDwing to vectors
                push!(D_wing, Dwing)
                push!(cD_wing, CDwing);
                push!(CdCD_vec, CdCD);

                vlm._addsolution(mainwing, "Cl/CL", ClCL)
                vlm._addsolution(mainwing, "Cd/CD", CdCD)

                #May need to pass this.

                y2b = 2*mainwing._ym/b
                # y2b = sum(y2b)             #sums each part of the the forces over the wing.
                push!(y2b_vec, y2b)
                push!(Gamma_vec, mainwing.sol["Gamma"])

                #Need Time vecotr, V vector, and W vector for this plot
                # figure(figname*"_3")
                # subplot(121)
                # plot([sim.t], [sim.vehicle.V[1]], ".r", label=L"V_x", alpha=0.5)
                # plot([sim.t], [sim.vehicle.V[2]], ".g", label=L"V_y", alpha=0.5)
                # plot([sim.t], [sim.vehicle.V[3]], ".b", label=L"V_z", alpha=0.5)
                # if PFIELD.nt==1; legend(loc="best", frameon=false); end;
                # subplot(122)
                # plot([sim.t], [sim.vehicle.W[1]], ".r", label=L"\Omega_x", alpha=0.5)
                # plot([sim.t], [sim.vehicle.W[2]], ".g", label=L"\Omega_y", alpha=0.5)
                # plot([sim.t], [sim.vehicle.W[3]], ".b", label=L"\Omega_z", alpha=0.5)
                # if PFIELD.nt==1; legend(loc="best", frameon=false); end;

                #Add Each componnet on V to the vecotr we return
                V = [sim.vehicle.V[1], sim.vehicle.V[2], sim.vehicle.V[3]];
                push!(V_vec,V);

                W = [sim.vehicle.W[1], sim.vehicle.W[2], sim.vehicle.W[3]];
                push!(W_vec, W);

                #
                # #ADDED CODE HERE.
                # #The forces aren't coming through, and I don't know why.
                #
                # #This Plot has a problem
                # figure(figname*"_4")
                # subplot(311)
                # println(sim.t, sum(Fx))
                # scatter(sim.t, sum(Fx), label="Lift forces on wing", alpha = 0.5, color = clr)
                # subplot(312)
                # scatter(sim.t, sum(Fy), label="Drag forces on wing", alpha = 0.5, color = clr)
                # subplot(313)
                # scatter(sim.t, sum(Fz), label="Side forces on wing", alpha = 0.5, color = clr)
                #
                #
                # #All works here except for the last ylabel. It doesn't show up.
                # figure(figname*"_5")
                # subplot(311)
                # plot(y2b, fx, label = "Lift per unit span", alpha = 0.5, color = clr)
                # subplot(312)
                # plot(y2b, fy, label = "Drag per unit span", alpha = 0.5, color = clr)
                # subplot(313)
                # plot(y2b, fz, label = "Side per unit span", alpha = 0.5, color = clr)
            end

            prev_wing = deepcopy(mainwing)

            return false

    end
    #initial conditions
    tinit = 0.0                                  #initial time
    Vinit = Vref*maneuver.Vvehicle(tinit/ttot)   #initial linear velocity
    Winit = zeros(3)                             #initial angular velocity
    RPMreference = RPM_Val

    #Everything is now defined! Let's create the simulation

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMreference, ttot; Vinit=Vinit, Winit=Winit, t=tinit);

    #Visulaize the simulation in paraview
    # files = uns.visualize_kinematics(   simulation, nsteps, save_path;
    #                                     run_name=run_name,
    #                                     prompt=false,
    #                                     paraview=false,
    #                                 )

    # run(`paraview --data="$save_path/$files"`)

    #The simulation also has some outputs we can recieve. The following code does this.

    # nullfunc(args...) = false
    pfield = uns.run_simulation(simulation, nsteps;
                                        surf_sigma=R/10,
                                        Vinf=Vinf,
                                        save_path=save_path,
                                        run_name=run_name,
                                        prompt=false,
                                        verbose=true,
                                        extra_runtime_function=monitor
                                )

    #run(`paraview --data="$save_path/$(files);tutorial_pfield...vtk"`)

    dt = 1/nsteps;
    time_vec = dt:dt:1;

    return cD_wing,cL_wing,D_wing,L_wing,time_vec,Forces_wing,forces_per_span, y2b_vec, CdCD_vec, mainwing, wake_coupled, Gamma_vec, V_vec, W_vec
end

function plot_Data(cD_wing, cL_wing, D_wing, L_wing, time_vec, Forces_wing,
    forces_per_span, y2b_vec, CdCD_vec, mainwing, wake_coupled, Gamma_vec, V_vec, W_vec, figname = "monitor_plot",
    figsize_factor = 1)
    #Goal, put all plotting here. clear from main.
    #make sure all variables you need are sent here as well.
    #Reset axis on Force plot to start from 0 and run to 1.1 * maximum(vector) on y

    #figure out how to implement colored points

    #Put the save paths in here as well.

    figure(figname, figsize=[7*2, 5*2]*figsize_factor)
    cla()
    subplot(221)
    xlim([0,1])
    xlabel(L"$\frac{2y}{b}$")
    ylabel(L"$\frac{Cl}{CL}$")
    title("Spanwise lift distribution")

    subplot(222)
    xlim([0,1])
    xlabel(L"$\frac{2y}{b}$")
    ylabel(L"$\frac{Cd}{CD}$")
    title("Spanwise drag distribution")

    subplot(223)
    xlabel("Simulation time (s)")
    ylabel(L"Lift Coefficient $C_L$")
    ylim([0,10]);

    subplot(224)
    xlabel("Simulation time (s)")
    ylabel(L"Drag Coefficient $C_D$")
    ylim([0,0.1]);

    figure(figname*"_2", figsize=[7*2, 5*1]*figsize_factor)
    cla()
    subplot(121)
    xlabel(L"$\frac{2y}{b}$")
    ylabel(L"Circulation $\Gamma$")
    subplot(122)
    xlabel(L"$\frac{2y}{b}$")
    ylabel(L"Effective velocity $V_\infty$")

    figure(figname*"_3", figsize=[7*2, 5*1]*figsize_factor)
    cla()
    subplot(121)
    xlabel("Simulation time")
    ylabel("Velocity")
    subplot(122)
    xlabel("Simulation time")
    ylabel("Angular velocity")

    figure(figname*"_4", figsize=[7*2, 5*1]*figsize_factor)
    cla()
    subplot(311)
    ylabel("Lift Force")
    subplot(312)
    ylabel("Drag Force")
    subplot(313)
    xlabel("Simulation time")
    ylabel("Side Force")
    ylim([0,10])

    figure(figname*"_5", figsize = [7*2, 5*1]*figsize_factor)
    cla()
    subplot(311)
    ylabel("Lift per unit span")
    subplot(312)
    ylabel("Drag per unit span")
    subplot(313)
    xlabel("Simulation time")
    ylabel("Side per unit span")

    #PLOTS ARE NOW MADE, LETS ADD THE POINTS

    #Plot data on figure
    size = length(y2b_vec[1]);

    figure(figname)
    subplot(221)
    for i = 1:size
        aux = (i/size);
        clr = (1-aux, 0, aux)
        plot(y2b_vec[i], cL_wing[i], "-", label="FLOWVLM", alpha=0.5, color = clr)
    end

    subplot(222)
    for i = 1:size
        aux = (i/size);
        clr = (1-aux, 0, aux)
        plot(y2b_vec[i], CdCD_vec[i], "-", label="FLOWVLM", alpha=0.5, color = clr)
    end

    subplot(223)
    size = length(time_vec)
    for i = 1:size
        aux = (i/size);
        clr = (aux, 0, 1-aux)
        plot(time_vec[i], norm(sum(cL_wing[i,:])), "o", label="FLOWVLM", alpha=0.5, color=clr)
    end

    subplot(224)
    for i = 1:size
        aux = (i/size);
        clr = (aux, 0, 1-aux);
        plot(time_vec[i], norm(sum(cD_wing[i,:])), "o", label="FLOWVLM", alpha=0.5, color = clr)
    end

    figure(figname*"_2")
    subplot(121)

    for i = 1:size
        cratio = i/size;
        y = Gamma_vec[i]
        clr = [cratio,0,1-cratio]
        plot(y2b_vec[i], y, "-", label="FLOWVLM", alpha=0.5, color = clr)
    end

    #Figure this one out. We may only need the wake_coupled part, which may need to be passed.
    if wake_coupled
        subplot(122)
        for i = 1:size
            plot(y2b_vec, norm.(mainwing.sol["Vkin"])/magVinf, "-", label="FLOWVLM", alpha=0.5)
            #add color = [clr[1], 1, clr[3]]
            if VehicleType==uns.VLMVehicle
                plot(y2b_vec, norm.(mainwing.sol["Vvpm"]), "-", label="FLOWVLM", alpha=0.5)
                #add color = clr
            end
            plot(y2b_vec, [norm(Vinf(vlm.getControlPoint(mainwing, i), T)) for i in 1:vlm.get_m(mainwing)],
                                                        "-k", label="FLOWVLM", alpha=0.5)
        end
    end
    #

    # println("V_vec[1] = ", length(V_vec));

    size = length(V_vec)
    # println("V_vec: ", size)
    # println("W_vec: ", length(W_vec))
    # println("t: ", length(time_vec))

    figure(figname*"_3")
    subplot(121)
    # for i = 1:size
    #     # println(time_vec[i], "  ", V_vec[i][1]);
    #     # println(V_vec[i])
    #     # plot(time_vec[i], V_vec[i], [".r",".g",".b"], label = [L"V_x", L"V_y",L"V_z",], alpha=0.5)
    #     plot(time_vec[i], V_vec[i][1], ".r", label=L"V_x", alpha=0.5)
    #     plot(time_vec[i], V_vec[i][2], ".g", label=L"V_y", alpha=0.5)
    #     plot(time_vec[i], V_vec[i][3], ".b", label=L"V_z", alpha=0.5)
    #     # if PFIELD.nt==1;
    #     legend(loc="best", frameon=false);
    #     # end;
    #     subplot(122)
    #     # println(W_vec[i])
    #     # plot(time_vec[i], W_vec[i], [".r",".g",".b"], label=[L"\Omega_x",L"\Omega_y",L"\Omega_z"], alpha = 0.5)
    #     plot(time_vec[i], W_vec[i][1], ".r", label=L"\Omega_x", alpha=0.5)
    #     plot(time_vec[i], W_vec[i][2], ".g", label=L"\Omega_y", alpha=0.5)
    #     plot(time_vec[i], W_vec[i][3], ".b", label=L"\Omega_z", alpha=0.5)
    #     # if PFIELD.nt==1;
    #     legend(loc="best", frameon=false);
    #     # end;
    # end
    # V_x = V_vec[:][1]
    # V_y = V_vec[2]
    # V_z = V_vec[3]
    # println(V_x)
    for i = 1:size
        # println(V_vec[i]);
        # println(V_vec[1][1]);
        plot(time_vec[i], (V_vec[i])[1], ".r", label=L"V_x", alpha=0.5);
        plot(time_vec[i], (V_vec[i])[2], ".g", label=L"V_y", alpha=0.5);
        plot(time_vec[i], (V_vec[i])[3], ".b", label=L"V_z", alpha=0.5);
        legend(loc="best", frameon=false);
        subplot(122)
        plot(time_vec[i], (W_vec[i])[1], ".r", label=L"W_x", alpha=0.5);
        plot(time_vec[i], (W_vec[i])[2], ".g", label=L"W_x", alpha=0.5);
        plot(time_vec[i], (W_vec[i])[3], ".b", label=L"W_x", alpha=0.5);
        legend(loc="best", frameon=false);
    end

    figure(figname*"_4")
    size = length(time_vec)
    for i = 1:size
        cratio = i/size;
        clr = [cratio,0,1-cratio];
        # println(time_vec[i], [Forces_wing[i]][1])
        subplot(311)
        scatter(time_vec[i], (Forces_wing[i])[1], label="Lift forces on wing", alpha = 0.5, color = clr)
        #add color = clr
        subplot(312)
        scatter(time_vec[i], (Forces_wing[i])[2], label="Drag forces on wing", alpha = 0.5, color = clr)
        #add color = clr
        subplot(313)
        scatter(time_vec[i], (Forces_wing[i])[3], label="Side forces on wing", alpha = 0.5, color = clr)
        #add color = clr
    end
    #
    figure(figname*"_5")
    size = length(y2b_vec)
    for i = 1:size
        cratio = i/size;
        clr = [cratio,0,1-cratio];
        subplot(311)
        scatter(y2b_vec[i], (forces_per_span[i])[1], label = "Lift per unit span", alpha = 0.5, color = clr)
        subplot(312)
        plot(y2b_vec[i], (forces_per_span[i])[2], label = "Drag per unit span", alpha = 0.5, color = clr)
        subplot(313)
        plot(y2b_vec[i], (forces_per_span[i])[3], label = "Side per unit span", alpha = 0.5, color = clr)
    end
    #
    # #Set save file directory
    #
    # topdirectory = normpath(joinpath(@__DIR__)) #Make sure the path is set for inside of the current folder.
    # savepath = joinpath(topdirectory,savepath)  #add a folder into the current path to save the figures.
    # if !(isdir(savepath)); mkdir(savepath); end     #Make sure that savepath is there. Makes one if not
    #
    # figname = "monitor_plot"
    #
    # figure(figname)
    # savefig(joinpath(savepath, figname*"_"*today*".pdf"))
    #
    # figure(figname*"_2")
    # savefig(joinpath(savepath,figname*"_2_"*today*".pdf"))
    #
    # figure(figname*"_3")
    # savefig(joinpath(savepath,figname*"_3_"*today*".pdf"))
    #
    # figure(figname*"_4")
    # savefig(joinpath(savepath, figname*"_4_"*today*".pdf"))
    #
    # figure(figname*"_5")
    # savefig(joinpath(savepath,figname*"_5_"*today*".pdf"))
end

cD_wing,cL_wing,D_wing,L_wing,time_vec,Forces_wing,forces_per_span, y2b_vec, CdCD_vec, mainwing, wake_coupled, Gamma_vec, V_vec, W_vec = main()

if disp_plot == true
    plot_Data(cD_wing, cL_wing, D_wing, L_wing, time_vec, Forces_wing, forces_per_span, y2b_vec, CdCD_vec, mainwing, wake_coupled, Gamma_vec, V_vec, W_vec)
end

max_stars_repo_name: brigoste/FlowVLMForces

max_issues_repo_name: brigoste/FlowVLMForces

max_forks_repo_name: brigoste/FlowVLMForces


text:
#Code written by Skipper af Rosenborg and Esben Bager
using StatsBase
using DataFrames
using CSV

#Set training size
#trainingSizeInput = parse(Int64, ARGS[1])
trainingSize = 12

#Set path
#path = "/zhome/9f/d/88706/SpecialeCode/Thesis/ML/Lasso_Test"
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/ML"
#cd("$(homedir())/Documents/GitHub/Thesis/Data")
#path = "$(homedir())/Documents/GitHub/Thesis/Data"
@everywhere include("ParallelModelGeneration.jl")
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Setup data
possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

modelMatrix = zeros(industriesTotal, possibilities)
testModel = [0 1 0 0 0]
for i=2:industriesTotal
    modelMatrix[i, :] = testModel
end

#Set output path
#path = "$(homedir())/Documents/GitHub/Thesis/Data/IndexDataDiff/"

#Load data
XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)

riskAversions = linspace(0, 2.4, 10)
XArrays, YArrays = generateXandYs(industries, modelMatrix)

nRows = size(XArrays[1])[1]
w1N = repeat([0.1], outer = 10) #1/N weights
return1NMatrix = zeros(nRows-trainingSize)
returnSAAMatrix = zeros(nRows-trainingSize)

weightsSAA = zeros(nRows-trainingSize, 10)
forecastRows = zeros(nRows-trainingSize, 10)

startPoint = 241 #194608
endPoint = 1080 #201607

for g = 1:10
    fileName = "Results"
    gammaRisk = riskAversions[g] #riskAversion in MV optimization
    total = (endPoint-trainingSize)
    for t = (startPoint-trainingSize):(endPoint-trainingSize)
        println("time $t / $total, gammaRisk $g / 10 ")
        trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)
        expectedReturns = zeros(industriesTotal)
        #Get expected return as mean of past training data
        for i=1:10
            expectedReturns[i] = mean(trainingYArrays[i][:])
        end

        Sigma =  cov(trainingXArrays[1][:,1:10])

        #A=U^(T)U where U is upper triangular with real positive diagonal entries
        F = lufact(Sigma)
        U = F[:U]  #Cholesky factorization of Sigma

        #getting the actual Y values for each industry
        valY = zeros(10)
        for i = 1:10
            valY[i] = validationY[i][1]
        end
        return1N, returnSAA, wSAA, forecastRow = performMVOptimization(expectedReturns[:], U, gammaRisk, valY, valY)
        weightsSAA[t, 1:10]    = wSAA
        forecastRows[t, 1:10]  = forecastRow
        return1NMatrix[t]      = return1N
        returnSAAMatrix[t]     = returnSAA
    end
    fileName = fileName*"_train"*string(trainingSize)*"_"*string(gammaRisk)
    writedlm(fileName*"saaWeights.csv", weightsSAA,",")
    writedlm(fileName*"SAAReturns.csv", returnSAAMatrix,",")
    writedlm(fileName*"saa1N.csv", return1NMatrix,",")
end

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/ML"
#cd(path)
println("Finished")

ext: jl

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
using Dolo

model = Model("examples/models/rbc.yaml")

F = Dolo.Euler(model, ignore_constraints=false);

sol = improved_time_iteration(model, details=false)

# P, μ = Dolo.ergodic_distribution(model, sol.dr, F.grid_exo, F.grid_endo, F.dprocess);

x0 = Dolo.MSM([sol.dr(i,F.s0) for i in 1:length(sol.dr.grid_exo.nodes)])

@time Π, dΠ = Dolo.new_transition_dev(model, F.dprocess, x0, F.grid_exo, F.grid_endo);

using StaticArrays

m = @SMatrix [1.0 2.0 ; 3.0 4.0]
x = @SVector [0.4, 0.5]

1 ./ x .* m

m


text:
module AudioIO

# export the basic API
export play, stop, get_audio_devices

# default stream used when none is given
_stream = nothing

################## Types ####################

typealias AudioSample Float32
# A frame of audio, possibly multi-channel
typealias AudioBuf Array{AudioSample}

# used as a type parameter for AudioNodes. Subtypes handle the actual DSP for
# each node
abstract AudioRenderer

# A stream of audio (for instance that writes to hardware). All AudioStream
# subtypes should have a root and info field
abstract AudioStream
samplerate(str::AudioStream) = str.info.sample_rate
bufsize(str::AudioStream) = str.info.buf_size

# An audio interface is usually a physical sound card, but could
# be anything you'd want to connect a stream to
abstract AudioInterface

# Info about the hardware device
type DeviceInfo
    sample_rate::Float32
    buf_size::Integer
end

type AudioNode{T<:AudioRenderer}
    active::Bool
    end_cond::Condition
    renderer::T

    AudioNode(renderer::AudioRenderer) = new(true, Condition(), renderer)
    AudioNode(args...) = AudioNode{T}(T(args...))
end

function render(node::AudioNode, input::AudioBuf, info::DeviceInfo)
    # TODO: not sure if the compiler will infer that render() always returns an
    # AudioBuf. Might need to help it
    if node.active
        result = render(node.renderer, input, info)
        if length(result) < info.buf_size
            node.active = false
            notify(node.end_cond)
        end
        return result
    else
        return AudioSample[]
    end
end

# Get binary dependencies loaded from BinDeps
include( "../deps/deps.jl")
include("nodes.jl")
include("portaudio.jl")
include("sndfile.jl")
include("operators.jl")

############ Exported Functions #############

# Play an AudioNode by adding it as an input to the root mixer node
function play(node::AudioNode, stream::AudioStream)
    push!(stream.root, node)
    return node
end

# If the stream is not given, use the default global PortAudio stream
function play(node::AudioNode)
    global _stream
    if _stream == nothing
        _stream = PortAudioStream()
    end
    play(node, _stream)
end

# Stream not given, but sampling rate of the array is specified.
function play(node::AudioNode, sample_rate::Integer)
    stream = PortAudioStream(sample_rate)
    play(node, stream)
end

function stop(node::AudioNode)
    node.active = false
    notify(node.end_cond)
end

function Base.wait(node::AudioNode)
    if node.active
        wait(node.end_cond)
    end
end

function get_audio_devices()
    return get_portaudio_devices()
end

end # module AudioIO


text:
x, y = rand(300), rand(300)

ctest = ConvergentCrossMappingTest(timeseries_lengths = [50, 60, 75, 80], n_reps = 10)

# One set of cross mappings per time series length
@test causality(x, y, ctest) isa Vector{Vector{T}} where T

# Four time series lengths, so there should be four vectors of cross map values
@test causality(x, y, ctest) |> length == 4

# There should be ten cross map values per time series length
@test causality(x, y, ctest)[1] |> length == 10


text:
module SupportChecking

include("../src/Jaynes.jl")
using .Jaynes

model1 = @jaynes function foo()
    5 > 10 ? y ~ Normal(0.0, 1.0) : y ~ Bernoulli(0.5)
end

model2 = @jaynes function bar(x::Int)
    y ~ Normal(0.0, 1.0)
    y ~ Normal(0.5, 1.0)
    x > 10 ? y ~ Normal(0.0, 1.0) : y ~ Normal(0.3, 3.0)
end (DefaultPipeline)

end # module


text:
SameXYZLength = Union{Cube,Sphere}
SameXYLength  = Union{Cylinder,SphericalCap,SquarePyramid}
NoSameLength  = Union{Ellipsoid,EllipticCylinder,Parallelepiped,
                      RectangularPyramid}

halflengths(s::SameXYZLength   ) = Vec(s[1], s[1], s[1])
halflengths(s::SameXYLength    ) = Vec(s[1], s[1], s[2])
halflengths(s::NoSameLength    ) = Vec(s[1], s[2], s[3])
halflengths(s::HollowCylinder  ) = Vec(s[1], s[1], s[3])
halflengths(s::TSP             ) = Vec(s[1], s[1], s[2] * s[3])
halflengths(s::Ring            ) = Vec(s[1] + s[2], s[1] + s[2], s[3])
halflengths(s::Torus           ) = Vec(s[1] + s[2], s[1] + s[2], s[2])
halflengths(s::TriangularToroid) = Vec(s[1] + s[2], s[1] + s[2], s[3])

"""    in(p::Point, s::Shape)

Returns `true` if the point `p` lies inside or in the boundary of `s`,
or `false` otherwise.
"""
in(p::Point, s::Sphere            ) = p.x^2 + p.y^2 + p.z^2 ≤ s.r^2

in(p::Point, s::Ellipsoid         ) = (p.x / s.a)^2 + (p.y / s.b)^2 +
                                      (p.z / s.c)^2 ≤ 1

in(p::Point, s::Cylinder          ) = (p.x^2 + p.y^2 ≤ s.r^2) &&
                                      (abs(p.z) ≤ s.c)

in(p::Point, s::EllipticCylinder  ) = (p.x / s.a)^2 + (p.y / s.b)^2 ≤ 1 &&
                                      (abs(p.z) ≤ s.c)

in(p::Point, s::Cube              ) = (abs(p.x) ≤ s.a) && (abs(p.y) ≤ s.a) &&
                                      (abs(p.z) ≤ s.a)

in(p::Point, s::Parallelepiped    ) = (abs(p.x) ≤ s.a) && (abs(p.y) ≤ s.b) &&
                                      (abs(p.z) ≤ s.c)

in(p::Point, s::SquarePyramid     ) = (p.z + s.m * abs(p.x) ≤ s.b) &&
                                      (p.z + s.m * abs(p.y) ≤ s.b) &&
                                      (abs(p.z) ≤ s.b)

in(p::Point, s::RectangularPyramid) = (p.z + s.ma * abs(p.x) ≤ s.c) &&
                                      (p.z + s.mb * abs(p.y) ≤ s.c) &&
                                      (abs(p.z) ≤ s.c)

function in(p::Point, s::TSP)
    z = s.r * s.b
    b = 2 * s.b - z
    return (p.z + s.m * abs(p.x) ≤ b) &&
           (p.z + s.m * abs(p.y) ≤ b) &&
           (abs(p.z) ≤ z)
end


text:
module GRASP
using DocStringExtensions

export angularmomentum, parity # from module Symmetries
export nelectrons, nexcitations, maxelectrons

using AtomicLevels: AtomicLevels, Orbital, RelativisticOrbital, IntermediateTerm, degeneracy,
    κ2ℓ, κ2j
using HalfIntegers: HalfInteger, HalfInt

#
# Includes and imports related to compiled Fortran code.
# ------------------------------------------------------------------------------------------

let grasp_path_file = joinpath(@__DIR__, "../deps/grasp-path.jl")
    isfile(grasp_path_file) || error("deps/grasp-path.jl does not exist. Run `Pkg.build(\"GRASP\")`.")
    include(grasp_path_file)
end

const libgrasp_so = joinpath(@__DIR__, "..", "deps", "libgrasp.so")
isfile(libgrasp_so) || error("$(libgrasp_so) does not exist. Run `Pkg.build(\"GRASP\")`.")

#
# Additional includes and imports.
# ------------------------------------------------------------------------------------------

include("Symmetries.jl")
import .Symmetries: Parity, AngularMomentum, AngularSymmetry
import .Symmetries: parity, angularmomentum, absdiff

const SPECTROSCOPIC_NAMES = "s p d f g h i k l m n o q r t u v" |> split
function specname(l::Integer)
    @assert l >= 0
    SPECTROSCOPIC_NAMES[l+1]
end

"""
    nelectrons(obj)

Calculate the number of electrons in `obj`.
"""
function nelectrons end

"""
    nexcitations(from, to)

Calculate the number of excitations needed to go from `from` to `to`.
"""
function nexcitations end

"""
    maxelectrons(obj)

Return the maximum number of electrons.
"""
function maxelectrons end

include("rcsfs.jl")
include("Configurations.jl")
include("rmix.jl")
include("rwfn.jl")
include("rcsfparser.jl")

# Fixed Humanize.digitsep --- fixed for negative numbers
# from: https://github.com/IainNZ/Humanize.jl
function digitsep(value::Integer, sep = ",", k = 3)
    isnegative = value < 0
    value = string(abs(value))
    n = length(value)
    starts = reverse(collect(n:-k:1))
    groups = [value[max(x-k+1, 1):x] for x in starts]
    return (isnegative ? "-" : "") * join(groups, sep)
end

"""
    $(SIGNATURES)

Converts the energy value `x` from Hartrees to Kaysers (aka 1/cm),
using the conversion ``1~\\textrm{Ha} = 219474.63137~\\textrm{cm}^{-1}``.

# Examples

```jldoctest
julia> GRASP.hartree2kayser(1.0)
219474.63137
```
"""
hartree2kayser(x) = 219474.63137*x

"""
    $(SIGNATURES)

Like [`hartree2kayser`](@ref), but produces a humanized string with thousands separators.

# Examples

```jldoctest
julia> GRASP.h2k_humanize(1.0)
"219,475"
```
"""
h2k_humanize(x) = digitsep(round(Int, hartree2kayser(x)))

end # module


text:
module PlayerSeeker

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Player Seeker" => Ahorn.EntityPlacement(
        Maple.PlayerSeeker
    )
)

sprite = "decals/5-temple/statue_e.png"

function Ahorn.selection(entity::Maple.PlayerSeeker)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.PlayerSeeker, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end


text:
function execute_sts_weekly_api_call(user_model::AVKUserModel, stock_symbol::String; data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing)::(Union{T, Nothing} where T<:Any)

    # same error checks as daily api call
    # is user_model valid?
    check_result = _check_user_model(user_model)
    if (check_result != nothing && typeof(check_result.value) == AVKError)
        return check_result
    end

    # do we have the API key?
    check_result = _check_missing_api_key(user_model)
    if (check_result != nothing && typeof(check_result.value) == AVKError)
        return check_result
    end

    # do we have a stock_symbol -
    check_result = _check_missing_symbol(stock_symbol)
    if (check_result != nothing && typeof(check_result.value) == AVKError)
        return check_result
    end

    # get key's and id -
    api_key = user_model.alphavantage_api_key

    #use alpha_vantage_api to access data change function to monthly in this case
    #change this line for different time periods
    url = "$(alphavantage_api_url_string)?function=TIME_SERIES_WEEKLY&symbol=$(stock_symbol)&apikey=$(api_key)&datatype=$(string(data_type))&outputsize=$(string(outputsize))"
    api_call_result = _http_get_call_with_url(url)
    if (typeof(api_call_result.value) == AVKError)
        return api_call_result
    end#keep the same as daily call

    #check that result value is a string
    #this is also new but doesn't change between time periods
    if (typeof(api_call_result.value) isa String == false)
        println("Call result is not valid type")#error from not correct data output
    else
        api_call_raw_data = api_call_result.value
    end

    #call to logger
    if logger != nothing
        debug_log_api_call(logger, user_model, url)
    end
    if (data_type == :json)
        #process json
        #this string line should match the .json of the type we are calling
        data_series_key = "Weekly Time Series"
        return _process_raw_json_api_data_sts(api_call_raw_data, data_series_key)

    elseif(data_type == :csv)
        return _process_raw_csv_api_data(api_call_raw_data)
    else
        # formulate the error message -
        error_message = "$(data_type) is not supported. Supported values are {:json,:csv}"

        return Result{AVKError}(AVKError(error_message))
    end
end

max_stars_repo_name: teamalphavantage/AlphaVantageAPI

max_issues_repo_name: teamalphavantage/AlphaVantageAPI

max_forks_repo_name: teamalphavantage/AlphaVantageAPI


text:
module BenchSeqSimpleSum

import Enzyme
import Zygote
using BenchmarkTools
using FLoops

function sum_loop(xs)
    acc = 0.0
    for x in xs
        acc += x
    end
    return acc
end

function sum_floop(xs)
    @floop begin
        acc = 0.0
        for x in xs
            acc += x
        end
    end
    return acc
end

make_input() = [range(0, 1; length = 1000);]

zygote(sum, xs) = Zygote.gradient(sum, xs)
function enzyme(sum, xs)
    x̄s = similar(xs, float(eltype(xs)))
    Enzyme.autodiff(sum, Enzyme.Duplicated(xs, x̄s)), x̄s
end

function setup()
    suite = BenchmarkGroup()
    for ad in [zygote, enzyme]
        s1 = suite["ad=:$ad"] = BenchmarkGroup()
        for (impl, sum) in [:loop => sum_loop, :floop => sum_floop]
            s1["impl=:$impl"] = @benchmarkable($ad($sum, xs), setup = (xs = make_input()))
        end
    end
    return suite
end

clear() = nothing

end  # module


text:
# BigBed Writer
# =============

mutable struct WriterState
    # section info
    chromid::UInt32
    chromstart::UInt32
    chromend::UInt32

    # last record info
    chromstart_prev::UInt32
    chromend_prev::UInt32

    # section state
    started::Bool
    buffer::IOBuffer
    compressed::Vector{UInt8}
    intervals::Vector{Interval{Nothing}}

    # record buffer
    recordbuffer::IOBuffer

    # global state and stats
    blocks::Vector{BBI.Block}
    nfields::Int
    count::UInt64
    cov::UInt32
    min::Float32
    max::Float32
    sum::Float32
    ssq::Float32

    function WriterState()
        return new(
            # section info
            0, 0, 0,
            # last record info
            0, 0,
            # section state
            false, IOBuffer(), UInt8[], Interval{Nothing}[],
            # record buffer
            IOBuffer(),
            # global info
            BBI.Block[], 0, 0,
            0, +Inf32, -Inf32, 0.0f0, 0.0f0)
    end
end

struct Writer <: BioCore.IO.AbstractWriter
    # output stream
    stream::IO

    # number of zoom levels
    zoomlevels::Int

    # maximum size of uncompressed buffer
    max_buffer_size::UInt64

    # file offsets
    summary_offset::UInt64
    chrom_tree_offset::UInt64
    data_offset::UInt64

    # chrom name => (ID, length)
    chroms::Dict{String,Tuple{UInt32,UInt32}}

    # chrom ID => name
    chromnames::Dict{UInt32,String}

    # mutable state
    state::WriterState

    # zoom buffer
    zoombuffer::BBI.ZoomBuffer
end

const ZOOM_SCALE_FACTOR = 4

"""
    BigBed.Writer(output::IO, chromlist; binsize=64)

Create a data writer of the bigBed file format.

Arguments
---------

* `output`: data sink
* `chromlist`: chromosome list with length
* `binsize=64`: size of a zoom with the highest resolution

Examples
--------

```julia
output = open("data.bb", "w")
writer = BigBed.Writer(output, [("chr1", 12345), ("chr2", 9100)])
write(writer, ("chr1", 101, 150, "gene 1"))
write(writer, ("chr2", 211, 250, "gene 2"))
close(writer)
```
"""
function Writer(output::IO, chromlist::Union{AbstractVector,AbstractDict};
                binsize::Integer=64)
    # write dummy header (filled later)
    write_zeros(output, BBI.HEADER_SIZE)

    # write dummy zoom headers (filled later)
    chromlist_with_id = BBI.add_chrom_ids(chromlist)
    maxlen = Base.maximum(x[3] for x in chromlist_with_id)
    zoomlevels = BBI.determine_zoomlevels(maxlen, binsize, ZOOM_SCALE_FACTOR)
    write_zeros(output, BBI.ZOOM_HEADER_SIZE * zoomlevels)

    # write dummy total summary (filled later)
    summary_offset = position(output)
    write_zeros(output, BBI.SUMMARY_SIZE)

    # write chromosome B+-tree
    chrom_tree_offset = position(output)
    BBI.write_btree(output, chromlist_with_id)

    # write dummy data count (filled later)
    data_offset = position(output)
    write_zeros(output, sizeof(UInt64))

    # initialize zoom buffer (use temporary file?)
    max_block_size = 16 * 2^10
    chromlens = Dict(id => len for (name, id, len) in chromlist_with_id)
    zoombuffer = BBI.ZoomBuffer(chromlens, binsize, max_block_size)

    return Writer(
        output,
        zoomlevels,
        max_block_size,
        summary_offset,
        chrom_tree_offset,
        data_offset,
        Dict(name => (id, len) for (name, id, len) in chromlist_with_id),
        Dict(id => name for (name, id, len) in chromlist_with_id),
        WriterState(),
        zoombuffer)
end

function Base.write(writer::Writer, record::Tuple{String,Integer,Integer,Vararg})
    chromname, chromstart, chromend = record
    optionals = Base.tail(Base.tail(Base.tail(record)))
    chromid = writer.chroms[chromname][1]
    return write_impl(writer, chromid, UInt32(chromstart - 1), UInt32(chromend), optionals)
end

function Base.write(writer::Writer, interval::GenomicFeatures.Interval{Nothing})
    chromid = writer.chroms[seqname(interval)][1]
    return write_impl(writer, chromid, UInt32(leftposition(interval) - 1), UInt32(rightposition(interval)), ())
end

function Base.write(writer::Writer, interval::GenomicFeatures.Interval{T}) where T<:Tuple
    chromid = writer.chroms[seqname(interval)][1]
    return write_impl(writer, chromid, UInt32(leftposition(interval) - 1), UInt32(rightposition(interval)), GenomicFeatures.metadata(interval)...)
end

function Base.close(writer::Writer)
    state = writer.state
    if state.started
        finish_section!(writer)
    end

    # write data index
    stream = writer.stream
    data_index_offset = position(stream)
    BBI.write_rtree(writer.stream, state.blocks)

    # fill header
    header = BBI.Header(
        BBI.BED_MAGIC,
        3,  # version
        writer.zoomlevels,
        writer.chrom_tree_offset,
        writer.data_offset,
        data_index_offset,
        state.nfields,
        state.nfields,  # predefined field count (should be same as above?)
        0,  # autoSQL offset (0 for bigWig?)
        writer.summary_offset,
        writer.max_buffer_size,
        0   # reserved
    )
    seekstart(stream)
    write(stream, header)

    # fill summary
    seek(stream, writer.summary_offset)
    write(stream, BBI.Summary(state.cov, state.min, state.max, state.sum, state.ssq))

    # fill data count
    seek(stream, writer.data_offset)
    write(stream, state.count)

    # write zoom
    seekend(stream)
    zoomheaders = BBI.write_zoom(stream, writer.zoombuffer, writer.zoomlevels, ZOOM_SCALE_FACTOR)

    # fill zoom headers
    seek(stream, BBI.HEADER_SIZE)
    for zheader in zoomheaders
        write(stream, zheader)
    end

    close(stream)
    return
end

function BioCore.IO.stream(writer::Writer)
    return writer.stream
end

function write_impl(writer::Writer, chromid::UInt32, chromstart::UInt32, chromend::UInt32, optionals::Tuple)
    # infer the number of fields from the first record
    state = writer.state
    if state.nfields == 0
        state.nfields = 3 + length(optionals)
    end

    # check new record
    check_interval(state, chromid, chromstart, chromend, optionals)

    # write record to record buffer in order to estimate the section buffer size
    truncate(state.recordbuffer, 0)
    write(state.recordbuffer, chromid, chromstart, chromend)
    write_optionals(state.recordbuffer, optionals)
    write(state.recordbuffer, 0x00)

    if state.started && (chromid != state.chromid || position(state.buffer) + position(state.recordbuffer) > writer.max_buffer_size)
        finish_section!(writer)
    end
    if !state.started
        start_section!(writer, chromid, chromstart, chromend)
    end

    # write data to the buffer
    seekstart(state.recordbuffer)
    n = write(state.buffer, state.recordbuffer)

    # update state
    state.chromend = max(state.chromend, chromend)
    state.chromstart_prev = chromstart
    state.chromend_prev = chromend

    push!(state.intervals, Interval(writer.chromnames[chromid], chromstart + 1, chromend))

    return n
end

function check_interval(state::WriterState, chromid::UInt32, chromstart::UInt32, chromend::UInt32, optionals::Tuple)
    # NOTE: BigBed allows empty intervals.
    if chromid < state.chromid || (chromid == state.chromid && chromstart < state.chromstart_prev)
        throw(ArgumentError("disordered intervals"))
    end
    if state.nfields != 3 + length(optionals)
        throw(ArgumentError("inconsistent field counts"))
    end
end

# Write optional fields to stream.
function write_optionals(stream::IO, optionals::Tuple)
    n = length(optionals)
    if n ≥ 1  # name
        print(stream, optionals[1])
    end
    if n ≥ 2  # score
        print(stream, '\t', optionals[2])
    end
    if n ≥ 3  # strand
        print(stream, '\t', optionals[3])
    end
    if n ≥ 4  # thickstart
        print(stream, '\t', optionals[4] - 1)
    end
    if n ≥ 5  # thickend
        print(stream, '\t', optionals[5])
    end
    if n ≥ 6  # itemrgb
        print(stream, '\t')
        print_rgb(stream, optionals[6])
    end
    if n ≥ 7  # blockcount
        print(stream, '\t', optionals[7])
    end
    if n ≥ 8  # blocksizes
        print(stream, '\t')
        for x in optionals[8]
            print(stream, x, ',')
        end
    end
    if n ≥ 9  # blockstarts
        print(stream, '\t')
        for x in optionals[9]
            print(stream, x - 1, ',')
        end
    end
end

function print_rgb(stream::IO, rgb::ColorTypes.RGB)
    print_rgb(stream, convert(ColorTypes.RGB{N0f8}, rgb))
end

function print_rgb(stream::IO, rgb::ColorTypes.RGB{N0f8})
    print(stream, reinterpret(UInt8, rgb.r), ',', reinterpret(UInt8, rgb.g), ',', reinterpret(UInt8, rgb.b))
end

function start_section!(writer::Writer, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    state = writer.state
    @assert !state.started

    # initialize section info
    state.chromid = chromid
    state.chromstart = chromstart
    state.chromend = chromend

    # initialize last record info
    state.chromstart_prev = 0
    state.chromend_prev = 0

    # initialize data buffer
    truncate(state.buffer, 0)
    resize!(state.compressed, div(writer.max_buffer_size * 11, 10))

    # set flag
    state.started = true
    empty!(state.intervals)
    return
end

function finish_section!(writer::Writer)
    state = writer.state

    # write compressed data
    datasize = BBI.compress!(state.compressed, take!(state.buffer))
    offset = position(writer.stream)
    unsafe_write(writer.stream, pointer(state.compressed), datasize)

    # record block
    push!(
        state.blocks,
        BBI.Block((state.chromid, state.chromstart), (state.chromid, state.chromend), offset, datasize))

    # update global stats
    count, cov, min, max, sum, ssq = compute_section_summary(state.intervals)
    state.count += count
    state.cov += cov
    state.min = Base.min(state.min, min)
    state.max = Base.max(state.max, max)
    state.sum += sum
    state.ssq += ssq

    state.started = false
    return
end

function compute_section_summary(intervals::Vector{Interval{Nothing}})
    cov = 0
    min = Inf
    max = -Inf
    sum = 0.0
    ssq = 0.0
    for range in GenomicFeatures.coverage(intervals)
        depth = GenomicFeatures.metadata(range)
        size = range.last - range.first + 1
        cov += size
        min = Base.min(min, depth)
        max = Base.max(max, depth)
        sum += depth   * size
        ssq += depth^2 * size
    end
    return length(intervals), cov, min, max, sum, ssq
end

function write_zeros(stream::IO, n::Integer)
    for _ in 1:n
        write(stream, 0x00)
    end
    return n
end


text:
using OmegaCore
using OmegaTest
using OmegaCore.Traits

using Test
using Distributions
using Random

function testsimple1()
  ω = defΩ()()
  ω = OmegaCore.tagrng(ω, Random.GLOBAL_RNG)
  ω2 = defΩ()()
  ω2 = OmegaCore.tagrng(ω2, Random.GLOBAL_RNG)
  x = 1 ~ StdNormal{Float64}()
  map(x, [ω, ω2])
end

function testsimple2()
  x = 1 ~ StdNormal{Float64}()
  rng = Random.GLOBAL_RNG
  OmegaCore.randsample(rng, x, 1; alg = RejectionSample)
end

function testsimple3a()
  x = 1 ~ StdNormal{Float64}()
  ω = defω()
  map(x, typeof(ω)[ω, ω])
end

function testsimple3()
  x = 1 ~ StdNormal{Float64}()
  ω = defω()
  x.([ω, ω])
end

@testset "infer types" begin
  @test isinferred(testsimple1)
  @test isinferred(testsimple2)
  @test isinferred(testsimple3)
  @test isinferred(testsimple3a)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
"""
src contains dst
"""
struct Contains <: Relationship
    src
    dst
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using StatsBase
using Plots
using DataFrames
using Underscores
using CSV
using HTTP, JSON
using BSON
using Dates
## ============================================================================
#                    apple
#  ============================================================================
regions = ["Hungary", "United Kingdom", "France", "Sweden", "Denmark", "Germany", "Italy", "United States"]
old2new = Dict(
    "Hungary" => "Hungary",
    "United Kingdom" => "UK",
    "France" => "France",
    "Sweden" => "Sweden",
    "Denmark" => "Denmark",
    "Germany" => "Germany",
    "Italy" => "Italy",
    "United States" => "USA"
)

function replacemissing(xs)
    for i in eachindex(xs)
        if ismissing(xs[i])
            xs[i] = @_ vcat( xs[i-2:i-1], xs[i+1:i+2] ) |>
                skipmissing |>
                mean
        end
    end
    convert(Vector{Float64}, xs)
end

function custom_mean(xs)
    n = 0
    s = 0.
    for x in xs
        if !ismissing(x)
            n += 1
            s += x
        end
    end
    if n == 0
        @warn "all missing"
        return missing
    else
        m = s / n
        @assert isfinite(m)
        return m
    end
end

apple = @_ projectdir("data/apple/applemobilitytrends.csv") |>
    CSV.read(__, DataFrame) |>
    filter(:geo_type => ==("country/region") , __) |>
    filter(:region => ∈(regions) , __) |>
    # transform(__, :region => ByRow(x -> old2new[x]) => :region) |>
    select(__, Not(["geo_type", "alternative_name", "sub-region", "country"])) |>
    DataFrames.stack(__, Not([:region, :transportation_type]), variable_name="date", value_name="mobility") |>
    unstack(__, :transportation_type, :mobility) |>
    select(__, :date => ByRow(Date) => :date, Not(:date), AsTable([:driving,:walking,:transit]) => ByRow(x->custom_mean(x)) => :apple) |>
    rename(__, [:driving,:walking,:transit] .=> [:apple_driving,:apple_walking,:apple_transit])

## ============================================================================
#                    google
#  ============================================================================
const mobility_names = [
    "retail_and_recreation_percent_change_from_baseline",
    "grocery_and_pharmacy_percent_change_from_baseline",
    "workplaces_percent_change_from_baseline",
    "residential_percent_change_from_baseline"
]

const new_mobility_names = [
    "google_retail",
    "google_grocery",
    "google_workplaces",
    "google_residential"
]

google = @_ "/home/and/data/google-mobility/Global_Mobility_Report.csv" |>
    CSV.read(__, DataFrame) |>
    filter(:country_region => x -> x ∈ regions, __) |>
    filter( :sub_region_1 => ismissing, __) |>
    select(__, :date,
        :country_region => :region, #ByRow(x -> old2new[x]) =>
        mobility_names .=> new_mobility_names) |>
    select(__, All(), AsTable([:google_retail,:google_grocery,:google_workplaces]) => ByRow(mean) => :google)

## ============================================================================
#                    combine google and apple
#  ============================================================================
mobil = @_ outerjoin(apple, google; on=[:date,:region]) |>
    sort(__, :date) |>
    filter( Date("2020-02-15") <= _.date <= Date("2021-06-30"), __)
CSV.write(projectdir("data/mobility/", "mobility_international.csv"), mobil)

## regional data sets

function rollingmean(xs, window)
    n  = length(xs)-window+1
    ys = Vector{Float64}(undef, n)
    for i in 1:n
        x = xs[i:i+window-1]
        @assert sum( ismissing.(x) ) < 4 "nr of missing: $(sum( ismissing.(x) )) at i=$i"
        ys[i] = mean(skipmissing(x))
        @assert isfinite(ys[i]) "infinite at i=$i"
    end
    return ys
end

numericcols(df::DataFrame) =
    filter(x->nonmissingtype( eltype( df[!,x] ) ) <: Number, names(df))

function postprocessing(df, window)
    @assert isodd(window)
    Δ = window ÷ 2
    smoothed = DataFrame( :date => df.date[1+Δ:end-Δ] )
    cols = filter(x->nonmissingtype( eltype( df[!,x] ) ) <: Number, names(df)) #[eltype(col) <: Number for col in eachcol(mobil)]
    for col in cols
        @info col
        smoothed[!,col] = rollingmean(df[:,col], window)
        smoothed[!,col] ./= smoothed[1,col] == 0 ? 1 : abs(smoothed[1,col])
        smoothed[!,col] .-= 1
        smoothed[!,col] ./= StatsBase.std(smoothed[!,col])
        if !all( isfinite.(smoothed[!,col]) )
            @error "is not finite at $col:"
            println("$(findall(x-> !isfinite(x), smoothed[!,col]))")
        end
    end
    smoothed
end

function assertconsistency(df)
    @assert isequal( df.date |> length,
                     df.date[1]:Day(1):df.date[end] |> length)
    @assert all( [!any( @. isnan(col) | ismissing(col) | ! isfinite(col) ) for col in eachcol(df) if eltype(col) <: Number] )
    return df
end

function drop_empy_columns(df)
    drop = []
    for n in names(df)
        col = df[!,n]
        all(ismissing.(col)) && push!(drop, n)
    end
    return select(df, Not(drop))
end

regionalmobil = let
    regionalmobil = Dict()
    for region in regions
        @info region
        regional = @_ mobil |> filter(:region => ==(region), __) |>
            sort(__, :date) |>
            drop_empy_columns |>
            postprocessing(__, 7) |>
            select(__, All(), AsTable([:google,:apple]) => ByRow(mean) => :applegoogle) |>
            disallowmissing |>
            assertconsistency

        CSV.write(projectdir("data/mobility/", "mobility_country=$(old2new[region]).csv"), regional)
        regionalmobil[old2new[region]] = regional
    end
    BSON.bson(projectdir("data/mobility/", "mobility_international.bson"), regionalmobil)
    regionalmobil
end

## plot
using Plots
plotlyjs()
let
    ps = []
    nplots = length(regionalmobil)
    for (region, mob) in pairs(regionalmobil)
        mob = regionalmobil[region]
        p = plot(title="$region", xticks=:native)
        plot!(mob.date, mob.apple, lab="apple")
        plot!(mob.date, mob.google, lab="google")
        plot!(mob.date, mob.applegoogle, lab="apple-google")
        push!(ps,p)
    end
    plot(ps..., layout=(nplots,1), sharex=true, link=:x, size=(1000, nplots*250))
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
#=
https://github.com/JuliaPy/SymPy.jl/blob/master/examples/tutorial.md
plot(ex1, ex2, a, b) will plot a parametric plot of the two expressions over the interval [a,b]
=#
using Plots
using SymPy

gr()

@vars x

p1 = plot(
    cos(x), sin(x), 0, 2pi,
    aspect_ratio=1,legend=false
)

p2 = plot(
    sin(2x), cos(3x), 0, 4pi,
    aspect_ratio=1, legend=false
)

plot(p1, p2, layout=(1, 2))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Test
using SearchLight
using Item

@testset "Item unit tests" begin
    @test 1 == 1
end

max_stars_repo_name: jaantollander/genie-webapp-db

max_issues_repo_name: jaantollander/genie-webapp-db

max_forks_repo_name: jaantollander/genie-webapp-db


text:
# https://en.wikipedia.org/wiki/Kernel_(statistics)

## these functions inspired by KernelFunctions.jl
# --- begin
abstract type AbstractKernel end

"""
    (k::AbstractKernel)(x, xg)

Evaluate the kernel for the pair of values `(x,xg)`.
"""
(k::AbstractKernel)(x, xg) = kern(k, x, xg)
kern(k::AbstractKernel, x, xg) = scaling(k) * unscaledkern(k,x,xg)
bandwidth(k::AbstractKernel) = k.h
scaling(k::AbstractKernel) = k.constant / k.h

function Base.show(io::IO, k::T) where {T<:AbstractKernel}
    return print(io, "$T (h = ", k.h, ")")
end
# --- end

# Gaussian
struct GaussianKernel{T<:Real,S<:Real} <: AbstractKernel
    h::T
    constant::S
end
GaussianKernel(h::T; constant::S=1/sqrt(2*pi)) where {T,S} = GaussianKernel{T,S}(h,constant)
unscaledkern(k::GaussianKernel, x, xg) =  exp(-((x-xg)/k.h)^2 / 2)

"""
    NormalKernel(h)

Alias of [`GaussianKernel`](@ref).
"""
const NormalKernel = GaussianKernel

# Uniform
struct UniformKernel{T<:Real,S<:Real} <: AbstractKernel
    h::T
    constant::S
end
UniformKernel(h::T; constant::S=1/2) where {T,S} = UniformKernel{T,S}(h,constant)
unscaledkern(k::UniformKernel, x::T, xg) where {T} = abs((x-xg)/k.h) <= 1 ? one(T) : zero(T)

# Epanechnikov
struct EpanechnikovKernel{T<:Real,S<:Real} <: AbstractKernel
    h::T
    constant::S
end
EpanechnikovKernel(h::T; constant::S=3/4) where {T,S} = EpanechnikovKernel{T,S}(h,constant)
unscaledkern(k::EpanechnikovKernel, x::T, xg) where {T} = abs((x-xg)/k.h) <= 1 ? 1 - ((x-xg)/k.h)^2 : zero(T)


text:
module AbInitioSoftwareBase

abstract type AbInitioSoftware end

include("fileops.jl")
include("Inputs.jl")
include("Commands.jl")

end


text:
export center_extract, center_set!, get_indices_around_center, center_pos
export rr, rr_3D
export fftpos
export normabs2
export get_index_for_pos

"""
    get_index_for_pos(pos_arr, pos)

For a arr containing `range` like data,
return the `index` which is closest to the `pos` value

# Example
```julia-repl
julia> pos = collect(fftpos(1, 5))
5-element Vector{Float64}:
 -0.5
 -0.25
  0.0
  0.25
  0.5

julia> get_index_for_pos(pos, 0.124)
3

julia> get_index_for_pos(pos, 0.1251)
4

julia> pos = collect(fftpos(1, 6))
6-element Vector{Float64}:
 -0.5
 -0.3333333333333333
 -0.16666666666666666
  2.2204460492503052e-17
  0.16666666666666669
  0.33333333333333337

julia> get_index_for_pos(pos, 0)
4

julia> get_index_for_pos(pos, 0.1)
5
```
"""
function get_index_for_pos(pos_arr, pos)
    left = pos_arr[begin]
    right = pos_arr[end]
    Δ = (right - left) / (length(pos_arr) - 1)

    ind = 1 + round(Int, (pos - left) / Δ)
    return ind
end

"""
    normabs2(arr)

Takes abs2 and normalizes maximum to 1

 # Examples
```jldoctest
julia> normabs2([10.0, 12.1])
2-element Array{Float64,1}:
 0.6830134553650707
 1.0
```
"""
function normabs2(arr)
    arr2 = abs2.(arr)
    return arr2 ./ maximum(arr2)
end

"""
    fftpos(L, N)

Construct a range from -L/2 to L/2.
However, we ensure that everything is centered around the center
in a way that a FFT interpretes it correctly.
For odd sequences it is indeed in the real center.
For even sequences the center is at `N/2 + 1`.

 # Examples
```jldoctest
julia> collect(fftpos(1, 4))
4-element Array{Float64,1}:
 -0.5
 -0.25
  0.0
  0.25

julia> collect(fftpos(1, 5))
5-element Array{Float64,1}:
 -0.5
 -0.25
  0.0
  0.25
  0.5
```
"""
function fftpos(l, N)
    if iszero(l) || isone(N)
        return zero(l):zero(l)
    end
    if N % 2 == 0
        dx = l / N
        return range(-l/2, l/2-dx, length=N)
    else
        return range(-l/2, l/2, length=N)
    end
end

function fftpos_old(l, N)
    if iszero(l) && isone(N)
        return l:l
    end
    if N % 2 == 0
        dx = l / N
        return -l/2:dx:l/2-dx
    else
        dx = l / (N-1)
        return -l/2:dx:l/2
    end
end

"""
    get_indices_around_center(i_in, i_out)
A function which provides two output indices `i1` and `i2`
where `i2 - i1 = i_out`
The indices are chosen in a way that the set `i1:i2`
cuts the interval `1:i_in` in a way that the center frequency
stays at the center position.
Works for both odd and even indices
"""
function get_indices_around_center(i_in, i_out)
    if (mod(i_in, 2) == 0 && mod(i_out, 2) == 0
     || mod(i_in, 2) == 1 && mod(i_out, 2) == 1)
        x = (i_in - i_out) ÷ 2
        return 1 + x, i_in - x
    elseif mod(i_in, 2) == 1 && mod(i_out, 2) == 0
        x = (i_in - 1 - i_out) ÷ 2
        return 1 + x, i_in - x - 1
    elseif mod(i_in, 2) == 0 && mod(i_out, 2) == 1
        x = (i_in - (i_out - 1)) ÷ 2
        return 1 + x, i_in - (x - 1)
    end
end

"""
    center_extract(arr, new_size_array)
Extracts a center of an array.
`new_size_array` must be list of sizes indicating the output
size of each dimension. Centered means that a center frequency
stays at the center position. Works for even and uneven.
If `length(new_size_array) < length(ndims(arr))` the remaining dimensions
are untouched and copied.
# Examples
```jldoctest
julia> center_extract([1 2; 3 4], [1])
1×2 Array{Int64,2}:
 3  4
julia> center_extract([1 2; 3 4], [1, 1])
1×1 Array{Int64,2}:
 4
julia> center_extract([1 2 3; 3 4 5; 6 7 8], [2 2])
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
function center_extract(arr::AbstractArray, new_size_array)
    new_size_array = collect(new_size_array)

    # we construct two lists
    # the reason is, that we don't change higher dimensions which are not
    # specified in new_size_array
    out_indices1 = [get_indices_around_center(size(arr)[x], new_size_array[x])
                    for x = 1:length(new_size_array)]

    out_indices1 = [x[1]:x[2] for x = out_indices1]

    # out_indices2 contains just ranges covering the full size of each dimension
    out_indices2 = [1:size(arr)[i] for i = (1 + length(new_size_array)):ndims(arr)]
    return arr[out_indices1..., out_indices2...]
end

"""
    center_set!(arr_large, arr_small)
Puts the `arr_small` central into `arr_large`.
The convention, where the center is, is the same as the definition
as for FFT based centered.
Function works both for even and uneven arrays.
# Examples
```jldoctest
julia> center_set!([1, 1, 1, 1, 1, 1], [5, 5, 5])
6-element Array{Int64,1}:
 1
 1
 5
 5
 5
 1
```
"""
function center_set!(arr_large, arr_small)
    out_is = []
    for i = 1:ndims(arr_large)
        a, b = get_indices_around_center(size(arr_large)[i], size(arr_small)[i])
        push!(out_is, a:b)
    end

    #rest = ones(Int, ndims(arr_large) - 3)
    arr_large[out_is...] = arr_small

    return arr_large
end

"""
    center_pos(x)
Calculate the position of the center frequency.
Size of the array is `x`
# Examples
```jldoctest
julia> center_pos(3)
2
julia> center_pos(4)
3
```
"""
function center_pos(x::Integer)
    # integer division
    return div(x, 2) + 1
end

function rr_3D(s)
    rarr = zeros((s...))
    for k = 1:s[3]
        for j = 1:s[2]
            for i = 1:s[1]
                rarr[i, j, k] = sqrt( (i-center_pos(s[1]))^2 + (j-center_pos(s[2]))^2 + (k-center_pos(s[3]))^2)
            end
        end
    end
    return rarr
end

"""
    rr(s; L=nothing)

Generate a array with values being the distance to the center pixel.
`s` specifies the output size of the 2D array.
`L` specifies the length. Can be a single number indicating
the same length for both dimensions or a tuple specifying a length for each
dimension separately.

# Examples
```julia-repl
julia> PhysicalOptics.rr((4,4))
4×4 Matrix{Float64}:
 2.82843  2.23607  2.0  2.23607
 2.23607  1.41421  1.0  1.41421
 2.0      1.0      0.0  1.0
 2.23607  1.41421  1.0  1.41421

julia> PhysicalOptics.rr((4,4), L=1)
4×4 Matrix{Float64}:
 0.707107  0.559017  0.5   0.559017
 0.559017  0.353553  0.25  0.353553
 0.5       0.25      0.0   0.25
 0.559017  0.353553  0.25  0.353553

julia> PhysicalOptics.rr((4,4), L=(5, 1))
4×4 Matrix{Float64}:
 2.54951  1.34629  0.5   1.34629
 2.51247  1.27475  0.25  1.27475
 2.5      1.25     0.0   1.25
 2.51247  1.27475  0.25  1.27475

julia> PhysicalOptics.rr((5,5), L=(1, 1))
5×5 Matrix{Float64}:
 0.707107  0.559017  0.5   0.559017  0.707107
 0.559017  0.353553  0.25  0.353553  0.559017
 0.5       0.25      0.0   0.25      0.5
 0.559017  0.353553  0.25  0.353553  0.559017
 0.707107  0.559017  0.5   0.559017  0.707107
```
"""
function rr_old(s; L=nothing)
    # in case of odd/even array we need to fix the size
    s2 = s[2] %2 == 0 ? s[2] : s[2] -1
    s1 = s[1] %2 == 0 ? s[1] : s[1] -1
    x = fftpos(s2, s[2])'
    y = fftpos(s1, s[1])

    if isnothing(L) == false
        if typeof(L) <: Number
            x *= L ./ 2 ./ abs(x[1])
            y *= L ./ 2 ./  abs(y[1])
        else
            x *= L[1] ./ 2 ./ abs(x[1])
            y *= L[2] ./ 2 ./ abs(y[1])
        end
    end

    r = sqrt.(x.^2 .+ y .^ 2)
    return r
end

function rr(s, arr_type=Array; L=nothing)
    x, y = _get_rr_x_y(L, s, arr_type)
    r = sqrt.(x.^2 .+ y .^ 2)
    return r
end

function _get_rr_x_y(L::Number, s, arr_type)
    x = to_gpu_or_cpu(arr_type, fftpos(L, s[2])')
    y = to_gpu_or_cpu(arr_type, fftpos(L, s[1]))
    return x, y
end

function _get_rr_x_y(L, s, arr_type)
    x = to_gpu_or_cpu(arr_type, fftpos(L[2], s[2])')
    y = to_gpu_or_cpu(arr_type, fftpos(L[1], s[1]))
    return x, y
end

function _get_rr_x_y(L::Nothing, s, arr_type)
    s2 = s[2] %2 == 0 ? s[2] : s[2] -1
    s1 = s[1] %2 == 0 ? s[1] : s[1] -1
    x = to_gpu_or_cpu(arr_type, fftpos(s2, s[2])')
    y = to_gpu_or_cpu(arr_type, fftpos(s1, s[1]))
    return x, y
end

"""
    jinc(x)

Computes the jinc function which is \$\\text{jinc} = \\frac{J_1(x)}{x}\$
where \$J_1\$ being the first Bessel function.

In future, probably switch to: https://github.com/JuliaMath/SpecialFunctions.jl/issues/264

# Examples
```julia-repl
julia> jinc(0.0)
0.5

julia> jinc(0.0im)
0.5 + 0.0im

julia> jinc(0.0001im)
0.5000000006250005 - 3.0616170016954073e-17im

julia> jinc(1.0f0)
0.44005057f0

julia> jinc.([0.0 0.5 1.0])
1×3 Array{Float64,2}:
 0.5  0.484537  0.440051
```
"""
function jinc(x::T) where T
    if iszero(x)
		return convert(T, 0.5)
	else
		return besselj1(x) / x
	end
end

"""
    apply_rot_symmetry(arr, xs, xs_out, ys_out)

Takes a 1D array `arr` with rotational symmetry about the center point.
`arr` should have values evaluates at the positions `xs`.

`xs_out` and `ys_out` are the positions for the output array.
We interpolate `arr` and evaluate it at distances `sqrt(xs_out[j]^2+ys_out[i]^2)`.

"""
function apply_rot_symmetry(arr, xs, xs_out, ys_out)
    itp = Interpolations.scale(
            interpolate(arr,
                        BSpline(Cubic(Line(Interpolations.OnCell())))),
            xs)

    out = zeros(eltype(arr), (length(ys_out), length(xs_out)))

    for j in axes(out, 2)
        for i in axes(out, 1)
            r = sqrt(xs_out[j]^2+ys_out[i]^2)
            out[i, j] = itp(r)
        end
    end

    return out

end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
"""
    renderResult(img::Array{Float32}, predictions::Vector{YOLO.PredictLabel}, settings::YOLO.Settings; save_file::String = "")

Display image with predicted bounding boxes overlaid. Optionally save render as a file (should end in .jpg, .png etc.)
"""
function renderResult(img::Array{Float32}, predictions::Vector{YOLO.PredictLabel}, settings::YOLO.Settings; save_file::String = "")
    cols = distinguishable_colors(settings.num_classes+1, [RGB(1,1,1)])[2:end]
    img_w = settings.image_shape[1]
    img_h = settings.image_shape[2]
    img_rgb = imgslice_to_rgbimg(img)
    scene = Makie.Scene(resolution = size(img_rgb') .* 2 )
    Makie.image!(scene, img_rgb, scale_plot=false, show_axis = false, limits = Makie.FRect(0, 0, img_h, img_w))
    Makie.rotate!(scene, -0.5pi)

    for p in predictions
        alpha = clamp(p.conf,0.5,1.0)
        col = cols[p.class]
        pcol = RGBA(red(col),blue(col),green(col),alpha)
        x = p.bbox.x*img_w
        y = p.bbox.y*img_h
        w = p.bbox.w*img_w
        h = p.bbox.h*img_h
        points = [
            Point(y,x+w) => Point(y+h,x+w);
            Point(y,x) => Point(y+h,x);
            Point(y+h,x+w) => Point(y+h,x);
            Point(y,x+w) => Point(y,x);
            ]
        Makie.linesegments!(points, linewidth = 8, color=pcol)
        name = get(settings.numsdic,p.class,"")
        conf_rounded = round(p.conf, digits=2)
        Makie.text!(scene, "$name\n$(conf_rounded)",
             position = (y+3, x+3),
             align = (:left,  :top),
             color = pcol,
             textsize = 8,
             font = "Dejavu Sans"
         )
    end
    save_file != "" && Makie.save(save_file, scene)
    return scene
end

renderResult(img::KnetArray{Float32}, predictions::Vector{YOLO.PredictLabel}, settings::YOLO.Settings; save_file::String = "") =
renderResult(Array{Float32}(img), predictions, settings, save_file = save_file)


text:
function getbasketprice(basket)
    basket = JSON.parse(string("""{
        "date":"2021-06-20",
        "vegetables": [
            {
                "vegetable": "broccoli",
                "quantity": "2",
                "unit": "unit"
            },{
                "vegetable": "tomato",
                "quantity": "1",
                "unit": "kg"
            }
        ]
    }"""))

    #transform basket to a table
    basket_df = reduce(vcat, DataFrame.(basket["vegetables"]))

    #left join with a table of price
    price_table = CSV.read("data/crop plan - jm_yields_and_sales.csv", DataFrame)
    res = @pipe leftjoin(basket_df, price_table, on = :vegetable, makeunique = true) |>
        @transform(_, item_price = parse.(Int64, :quantity) .* :desired_selling_price)

    return sum(res[!, :item_price])
end


text:
module GIFTI

using LightXML
using GeometryTypes: HomogenousMesh, Point3f0, Face, OffsetInteger
using CodecZlib: transcode, GzipDecompressor
using Base64: base64decode

load(io::IO) = parse_gifti_mesh(parse_string(readstring(io)))
load(fname::String) = parse_gifti_mesh(parse_file(fname))

function parse_nifti_datatype(Tstr::String)
    if Tstr == "NIFTI_TYPE_INT32"
        Int32
    elseif Tstr == "NIFTI_TYPE_FLOAT32"
        Float32
    elseif Tstr == "NIFTI_TYPE_INT64"
        Int64
    elseif Tstr == "NIFTI_TYPE_FLOAT64"
        Float64
    else
        error("Unhandled data type: $(attribute(xml, "DataType"))")
    end
end

parse_nifti_datatype(xml::XMLElement) = parse_nifti_datatype(attribute(xml, "DataType"))

function parse_nifti_array_size(xml)
    N = parse(Int, attribute(xml, "Dimensionality"))
    ntuple(i -> parse(Int, attribute(xml, "Dim$(i-1)")), N)
end

function parse_nifti_array_data(xml::XMLElement)
    data_blobs = get_elements_by_tagname(xml, "Data")
    if length(data_blobs) == 0
        error("Could not find any `Data` elements")
    elseif length(data_blobs) > 1
        warn("Found multiple `Data` elements, using only the first one.")
    end
    data = content(first(data_blobs))
    encoding = attribute(xml, "Encoding")
    if encoding == "GZipBase64Binary"
        gzbytes = base64decode(data)
        return transcode(GzipDecompressor, gzbytes)
    elseif encoding == "Base64Binary"
    	return base64decode(data)
    else
        error("Unhandled encoding: $encoding")
    end
end

function parse_nifti_data_array(xml)
    T = parse_nifti_datatype(xml)
    dims = parse_nifti_array_size(xml)
    data = parse_nifti_array_data(xml)
    order = attribute(xml, "ArrayIndexingOrder")
    if order == "RowMajorOrder"
        array = reshape(reinterpret(T, data), reverse(dims))
    else
        @assert order = "ColumnMajorOrder"
        array = reshape(reinterpret(T, data), dims)
    end
    endian = attribute(xml, "Endian")
    if endian == "LittleEndian"
        array .= ltoh.(array)
    else
        @assert endian == "BigEndian"
        array .= btoh.(array)
    end
    array
end

function parse_gifti_mesh(xml::XMLElement)
    arrays = get_elements_by_tagname(xml, "DataArray");
    pointset_arrays = filter(get_elements_by_tagname(xml, "DataArray")) do arr
        attribute(arr, "Intent") == "NIFTI_INTENT_POINTSET"
    end
    if length(pointset_arrays) == 0
        error("Could not find any arrays with intent NIFTI_INTENT_POINTSET")
    elseif length(pointset_arrays) > 1
        warn("Found multiple pointset arrays, using only the first.")
    end
    pointset = first(pointset_arrays)

    triangle_arrays = filter(get_elements_by_tagname(xml, "DataArray")) do arr
        attribute(arr, "Intent") == "NIFTI_INTENT_TRIANGLE"
    end
    if length(triangle_arrays) == 0
        error("Could not find any arrays with intent NIFTI_INTENT_TRIANGLE")
    elseif length(triangle_arrays) > 1
        warn("Found multiple triangle arrays, using only the first.")
    end
    triangles = first(triangle_arrays);

    vert_array = parse_nifti_data_array(pointset) ./ 100
    face_array = parse_nifti_data_array(triangles)
    vertices = reshape(reinterpret(Point3f0, vert_array), (size(vert_array, 2),))
    faces = reshape(reinterpret(Face{3, OffsetInteger{-1, Int32}}, face_array), (size(face_array, 2),));
    mesh = HomogenousMesh(vertices, faces)
end

parse_gifti_mesh(doc::XMLDocument) = parse_gifti_mesh(root(doc))

end


text:
using Luxor

include("environments/collect_points.jl")
include("tools/play.jl")

function kernel_test_initialize(environments, env_seed)

    tx = threadIdx().x
    bx = blockIdx().x

    input = @cuDynamicSharedMem(Float32, environments.number_inputs)
    offset = sizeof(input)

    sync_threads()

    initialize(tx, bx, input, environments, offset, env_seed)

    return
end

function main()

    config_environment = OrderedDict()
    config_environment["maze_columns"] = 5
    config_environment["maze_rows"] = 5
    config_environment["maze_cell_size"] = 80
    config_environment["agent_radius"] = 12
    config_environment["point_radius"] = 8
    config_environment["agent_movement_range"] = 10.0
    config_environment["use_sensors"] = true
    config_environment["reward_per_collected_positive_point"] = 500.0
    config_environment["reward_per_collected_negative_point"] = -700.0
    config_environment["number_time_steps"] = 1000

    env_seed = rand(1:1000)
    number_individuals = 100

    environments = CollectPoints(config_environment, number_individuals)

    shared_memory = get_memory_requirements(environments) + sizeof(Float32) * environments.number_inputs

    CUDA.@cuda threads = 10 blocks = number_individuals shmem = shared_memory kernel_test_initialize(environments, env_seed)

    CUDA.synchronize()

    println(environments.agents_positions[:,1])

    width = environments.maze_cell_size * environments.maze_columns
    height = environments.maze_cell_size * environments.maze_rows
    number_iterations = 100

    @play width height number_iterations begin
        render(environments)
    end

end

main()


text:
using Regularismus
using Base.Test

type countertype
  pass::UInt
  fail::UInt
end
counter = countertype(0, 0)

function verbose_handler(r::Test.Success)
    counter.pass += 1
    println("PASS ", r.expr)
end
function verbose_handler(r::Test.Failure)
    counter.fail += 1
    println("FAIL ", r.expr, " : ", r.resultexpr)
end
verbose_handler(r::Test.Error) = rethrow(r)

silent_handler(r::Test.Success) = counter.pass += 1
silent_handler(r::Test.Failure) = counter.fail += 1
silent_handler(r::Test.Error) = rethrow(r)

Test.with_handler(verbose_handler) do
    @test bots() == bots()
    @test wx"" == wx""
    @test alt(wx"", bots(), tops()) == alt(wx"", bots(), tops())
    @test cat(bots(), tops()) == cat(bots(), tops())
    @test int(star(bots()), tops()) == int(star(bots()), tops())
    @test plus(bots()) == plus(bots())
    @test botr() == botr()
    @test rx(wx"", wx"") == rx(wx"", wx"")
    @test alt(botr(), topr()) == alt(botr(), topr())
    @test cat(botr(), topr()) == cat(botr(), topr())
    @test star(botr()) == star(botr())
    @test plus(botr()) == plus(botr())
end

if counter.fail == 0
    info("Tests PASS: ", counter.pass)
    info("All tests pass")
elseif counter.pass == 0
    info("No tests run")
else
    info("Tests PASS: ", counter.pass)
    info("Tests FAIL: ", counter.fail)
    info("Some tests pass")
end


text:
using Compat

tests = [
    "common",
    "graph",
    "digraph"
    ]

for t in tests
    tp = joinpath(dirname(@__FILE__), "$(t).jl")
    println("Running $(tp)...")
    include(tp)
end


text:
"""
    execute_sts_search_api_call(user_model::PSUserModel, keyword::String;
        data_type::Symbol = :json, logger::Union{Nothing,AbstractLogger} = nothing) -> PSResult

This documentation will be awesome. the best description ever.
"""
function execute_sts_search_api_call(user_model::PSUserModel, keyword::String;
    data_type::Symbol = :json, logger::Union{Nothing,AbstractLogger} = nothing)::PSResult

    try

        # get key's and id -
        api_key = user_model.alphavantage_api_key

        # call to alpha_vantage_api to get data
        url = "$(alphavantage_api_url_string)?function=SYMBOL_SEARCH&keywords=$(keyword)&apikey=$(api_key)&datatype=$(string(data_type))"
        api_call_raw_data = http_get_call_with_url(url) |> checkresult

        # call to logger
        if (isnothing(logger) == false)
            log_api_call(logger, user_model, url)
        end

        # parse json if data called as a .json
        if (data_type == :json)

            # process json
            data_series_key = "bestMatches"
            return process_raw_json_data_sts_search_data(api_call_raw_data, data_series_key)

        elseif (data_type == :csv)

            # return process .csv
            return process_raw_csv_api_data(api_call_raw_data)
        else

            # tell user they requested an unsupported type of data
            error_message = "$(data_type) isn't supported by AlphaVantage. Supported values are {:json, :csv}"
            return PSResult{PSError}(PSError(error_message))
        end

    catch error
        return PSResult(error)
    end
end


text:
using IRTools, MacroTools, InteractiveUtils, Test
using IRTools: @dynamo, IR, meta, isexpr, xcall, self, insertafter!, recurse!,
  argument!, return!, func, var

@dynamo roundtrip(a...) = IR(a...)

@dynamo function passthrough(a...)
  ir = IR(a...)
  ir == nothing && return
  recurse!(ir)
  return ir
end

add(a, b) = a+b
@test roundtrip(add, 2, 3) == 5
@test passthrough(add, 2, 3) == 5

@test @code_ir(passthrough, add(2, 3)) isa IR

IRTools.refresh(roundtrip)

add(a, b) = a*b
@test roundtrip(add, 2, 3) == 6

relu(x) = x > 0 ? x : 0
@test roundtrip(relu, 1) == 1
@test passthrough(relu, 1) == 1

foo(x) = y = x > 0 ? x + 1 : x - 1
@test roundtrip(foo, 1) == 2
@test roundtrip(foo, -1) == -2
@test passthrough(foo, 1) == 2

@test passthrough(() -> [1, 2, 3]) == [1, 2, 3]

@dynamo function mullify(a...)
  ir = IR(a...)
  ir == nothing && return
  ir = MacroTools.prewalk(ir) do x
    x isa GlobalRef && x.name == :(*) && return GlobalRef(Base, :+)
    return x
  end
  for (x, st) in ir
    isexpr(st.expr, :call) || continue
    ir[x] = xcall(Main, :mullify, st.expr.args...)
  end
  return ir
end

@test mullify(prod, [5, 10]) == 15

@dynamo err(a...) = 0//0

@test_throws IRTools.Inner.CompileError err(+, 2, 3)

mutable struct Context
  calls::Int
end

@dynamo function (cx::Context)(a...)
  ir = IR(a...)
  ir == nothing && return
  recurse!(ir)
  calls = pushfirst!(ir, xcall(:getfield, self, QuoteNode(:calls)))
  calls = insertafter!(ir, calls, xcall(:+, calls, 1))
  insertafter!(ir, calls, xcall(:setfield!, self, QuoteNode(:calls), calls))
  return ir
end

@code_lowered Context(0)(add, 2, 3)

cx = Context(0)
@test cx(add, 2, 3.0) == 6
@test cx.calls > 5

ir = IR()
x = argument!(ir)
y = push!(ir, xcall(:*, x, x))
return!(ir, y)

@test IRTools.evalir(ir, 5) == 25

ir = IR()
x = argument!(ir)
y = argument!(ir)
z = push!(ir, xcall(:*, x, y))
return!(ir, z)

@test IRTools.evalir(ir, 5, 3) == 15

function pow(x, n)
  r = 1
  while n > 0
    n -= 1
    r *= x
  end
  return r
end

ir = @code_ir pow(2, 3)

ir[var(8)] = xcall(:+, var(5), var(2))

mul = func(ir)

@test mul(nothing, 10, 3) == 31


text:
f = x-> 0;
u_exact = x-> x;

a = 0;
b = 1;
n = 9;

problem = SparseFiniteDifferenceBVPProblem(0, 1, n, f, NeumannBC(1), NeumannBC(1));
assemble_system!(problem);
u = solve_bvp(problem);
norm(@. u - problem.x - u[1] )<1e-13


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMTruss.jl/blob/master/LICENSE

using FEMBase
using FEMBase.Test
using FEMTruss

@testset "FEMTrussTemperature.jl" begin
    E = 2e11
    A1 = 20e-4
    A2 = 15e-4
    X = Dict(1 => [3.0, 3.0],
             2 => [3.0, 0.0],
             3 => [0.0, 0.0],
             4 => [0.0, 3.0])
    el1 = Element(Seg2, [4,1])
    el2 = Element(Seg2, [3,1])
    el3 = Element(Seg2, [2,1])
    bel1 = Element(Poi1, [3])
    bel2 = Element(Poi1, [2])
    bel3 = Element(Poi1, [4])
    #bel4 = Element(Poi1, [1])
    elements = [el1, el2, el3, bel1, bel2, bel3]
    update!(elements, "geometry", X)

    update!([bel1, bel2, bel3], "fixed displacement 1", 0.0)
    update!([bel1, bel2, bel3], "fixed displacement 2", 0.0)
    #update!(bel4, "nodal force 2", 10)
    update!(elements, "youngs modulus", E)
    update!([el1,el3], "cross section area", A2)
    update!([el2], "cross section area", A1)

    lambda = 1/ 75000
    dt = 40
    update!([el2], "thermal expansion coefficient", lambda)
    update!([el2], "temperature difference", dt)

    problem = Problem(Truss, "test problem", 2)
    add_elements!(problem, elements)

    step = Analysis(Static)
    add_problems!(step, [problem])
    ls, normu, normla = run!(step)
    # Check deflections in top
    u_oracle = [0.00077645, 0.00077645]
    r_oracle = [0.0, -77645.0, +77645.0, +77645.0, -77645.0, 0.0]
    @test isapprox(ls.u[collect(1:2)], u_oracle,rtol=1e-5)
    @test isapprox(ls.la[collect(3:8)], r_oracle*-1.0, rtol=1e-5)

# Find element forces, these
    f1_oracle = [-77645.0, 77645.0]
    f2_oracle = [1.09807e5, -1.09807e5]
    K1,T1 = get_Kg_and_Tg(el1, 2,2, 0.0)
    f1 = get_F_local(el1,2,2, 0.0)
    f_loc1 = T1*K1*ls.u[[3,4, 1,2]]  +f1
    K2,T2 = get_Kg_and_Tg(el2, 2,2, 0.0)
    f2 = get_F_local(el2,2,2, 0.0)
    f_loc2 = T2*K2*ls.u[[5,6,1,2]] + f2
    @test isapprox(f_loc1, f1_oracle,rtol=1e-5)
    @test isapprox(f_loc2, f2_oracle,rtol=1e-5)

    # Lets verify the system
    #Stiffness Matrix assemble by hand
end


text:
function genInitialState(c::Constants, d::Data)
  J = d.J
  K = c.K
  L = c.L
  I = d.I
  N = d.N

  vec_y = vcat(vec.(d.y)...)
  y_neg = filter(y_inj -> !isnan(y_inj) && y_inj < 0, vec_y)

  y_imputed = begin
    local out = [zeros(Float64, N[i], J) for i in 1:I]
    grid_size = 30
    y_grid = collect(range(-7, 0, length=grid_size))
    miss_probs = [[prob_miss(yg, c.beta[:, i]) for yg in y_grid] for i in 1:I]
    for i in 1:I
      for n in 1:N[i]
        for j in 1:J
          if isnan(d.y[i][n, j])
            # out[i][n, j] = rand(Uniform(y_lower, y_upper))
            out[i][n, j] = wsample(y_grid, miss_probs[i])
          else
            out[i][n, j] = d.y[i][n, j]
          end
          @assert !isnan(out[i][n, j])
        end
      end
    end

    out
  end

  alpha = mean(c.alpha_prior)
  v = [mean(Beta(alpha, 1)) for k in 1:K]
  b = cumprod(v)
  Z = [ Bool(rand(Bernoulli(b[k]))) for j in 1:J, k in 1:K ]
  delta = Dict(false => rand(c.delta_prior[0], L[0]),
               true  => rand(c.delta_prior[1], L[1]))
  sig2 = [rand(c.sig2_prior) for i in 1:I]
  W = Matrix{Float64}(hcat([rand(c.W_prior) for i in 1:I]...)')
  lam = [ Int8.(rand(Categorical(W[i, :]), N[i])) for i in 1:I ]
  eta = begin
    function gen(z)
      arrMatTo3dArr([ rand(c.eta_prior[z]) for i in 1:I, j in 1:J ])
    end
    Dict([Bool(z) => gen(z) for z in 0:1])
  end
  gam = [zeros(Int8, N[i], J) for i in 1:I]
  for i in 1:I
    for j in 1:J
      for n in 1:N[i]
        z_lin = Z[j, lam[i][n]]
        gam[i][n, j] = rand(Categorical(eta[z_lin][i, j, :]))
      end
    end
  end

  eps = mean.(c.eps_prior)

  # Create State object
  state = State{Float64}()

  # Initialize State
  state.Z = Z
  state.delta = delta
  state.alpha = alpha
  state.v = v
  state.W = W
  state.sig2 = sig2
  state.eps = eps
  state.eta = eta
  state.lam = lam
  state.gam = gam
  state.y_imputed = y_imputed

  return state
end

include("SmartInit.jl")

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using Catlab
using Catlab.Graphs.BasicGraphs
using Colors
using FreeTypeAbstraction
using LaTeXStrings
using Luxor
using MathTeXEngine

import Catlab.CategoricalAlgebra.FreeSchema
import Catlab.CategoricalAlgebra.@acset
import Catlab.CategoricalAlgebra.@acset_type
import Catlab.CategoricalAlgebra.nparts
import Catlab.CategoricalAlgebra.subpart
import Catlab.CategoricalAlgebra.incident

########################
# FUNCTION DEFINITIONS #
########################

"""
    latex_text_size(lstr::LaTeXString)
Returns the width and height of a latex string.
(Credit: Created by Davi Sales Barreira)
"""
function latex_text_size(lstr::LaTeXString)
    sentence = generate_tex_elements(lstr)
    els = filter(x -> x[1] isa TeXChar, sentence)
    chars = [x[1].char for x in els]
    fonts = [x[1].font for x in els]
    pos_x = [x[2][1] for x in els]
    pos_y = [x[2][2] for x in els]
    scales = [x[3] for x in els]
    extents = [FreeTypeAbstraction.get_extent(f, c) for (f, c) in zip(fonts, chars)]

    textw = []
    texth = []
    for i = 1:length(extents)
        textw = push!(
            textw,
            FreeTypeAbstraction.height_insensitive_boundingbox(
                extents[i],
                fonts[i],
            ).widths[1] * scales[i] + pos_x[i],
        )
        texth = push!(
            texth,
            FreeTypeAbstraction.height_insensitive_boundingbox(
                extents[i],
                fonts[i],
            ).widths[2] * scales[i] - pos_y[i],
        )
    end
    return (maximum(textw), maximum(texth))
end

"""
    Luxor.text(lstr::LaTeXString, valign=:baseline, halign=:left; kwargs...)
Draws LaTeX string using `MathTexEngine.jl`. Hence, uses ModernCMU as font family.
Note that `valign` is not perfect.
This function assumes that the axis are in the standard Luxor directions,
i.e. ↓→.
(Credit: Created by Davi Sales Barreira)
"""
function Luxor.text(
    lstr::LaTeXString,
    pt::Point;
    valign = :baseline,
    halign = :left,
    angle = 0,
)

    # Function from MathTexEngine
    sentence = generate_tex_elements(lstr)

    # Get current font size.
    textsize = get_fontsize()

    textw, texth = latex_text_size(lstr)

    if halign === :left
        translate_x = 0
    elseif halign === :right
        translate_x = -textw * textsize
    elseif halign === :center
        translate_x = -textw / 2 * textsize
    end

    if valign === :baseline
        translate_y = 0
    elseif valign === :top
        translate_y = -textsize
    elseif valign === :bottom
        translate_y = texth * textsize
    elseif valign === :middle
        translate_y = textsize / 2
    end

    # Writes text using ModernCMU font.
    for text in sentence
        if text[1] isa TeXChar
            @layer begin
                translate(translate_x, translate_y)
                fontface(text[1].font.family_name)
                fontsize(textsize * text[3])
                Luxor.text(
                    string(text[1].char),
                    pt + Point(text[2]...) * textsize * (1, -1),
                    angle = angle,
                )
            end
        elseif text[1] isa HLine
            @layer begin
                translate(translate_x, translate_y)
                pointstart = pt + Point(text[2]...) * textsize * (1, -1)
                pointend = pointstart + Point(text[1].width, 0) * textsize
                setline(1 * text[1].thickness * textsize)
                line(pointstart, pointend, :stroke)
            end
        end
    end
end

"""
    make_drawing(width, height, img_path, bkg_color, origin_p)

Creates a Luxor `Drawing` object with given width and height dimensions.
Saves the image to whatever path of choosing.
Allows one to set the background color of an image.
Finally, one can set Luxor's reference origin point
"""
function make_drawing(width, height, img_path, bkg_color, origin_p)
    d = Drawing(width, height, img_path)
    background(bkg_color)
    origin(origin_p)
    return d
end

"""
Converts a coordinate pair in the traditional Cartesian coordinate system to a `Point` in Luxor's coordinate system.
"""
cart2luxor(p) = Point(p.x, -p.y)
cart2luxor(p::Array) = Point(p[1], -p[2])
cart2luxor(x, y) = Point(x, -y)

"""
Converts a Luxor `Point` to a coordinate pair in the traditional Cartesian coordinate system.
"""
luxor2cart(p::Point) = [p.x, -p.y]

"""
Finds points along the parametric representation of a given line.
"""
parametric_point(t, x, y) = x + y * t

"""
Draws text at a given location along a morphism
"""
function morphtext(
    label1::Any,
    label2::Any;
    offset1::Real = 20,
    offset2::Real = 20,
    angle = 0,
    size = nothing,
)
    fontsize(size)
    text(label1, cart2luxor(0, offset1), valign = :middle, halign = :center, angle = angle)
    text(label2, cart2luxor(0, -offset2), valign = :middle, halign = :center, angle = angle)
end

function morphtext(label::Any; offset::Real = 20, angle = 0, size = nothing)
    fontsize(size)
    text(label, cart2luxor(0, offset), valign = :middle, halign = :center, angle = angle)
end

"""
Draws a single morphism between two objects.
"""
function morphism(
    initial_point::Array,
    terminal_point::Array;
    label1::Any = "",
    label2::Any = "",
    offset1::Real = 10,
    offset2::Real = 10,
    size::Real = 16,
    padding::Real = 0.15,
)

     v⃗ = terminal_point .- initial_point
       θ = atan(v⃗[2] / v⃗[1]) + (v⃗[1] < 0 ? π : 0)

    padded_initial_point =
        between(cart2luxor(initial_point), cart2luxor(terminal_point), padding)
    padded_terminal_point =
        between(cart2luxor(initial_point), cart2luxor(terminal_point), 1 - padding)

    arrow(
        padded_initial_point,
        padded_terminal_point,
        arrowheadlength = 15,
        decoration = 0.5,
        decorate = () -> morphtext(
            label1,
            label2,
            offset1 = offset1,
            offset2 = offset2,
            angle = θ,
            size = size,
        ),
        linewidth = 3,
    )

end

"""
Draws two morphism between two objects.
Morphisms can face different directions and be styled differently.
"""
function morphism(
    initial_point_1::Array,
    terminal_point_1::Array,
    initial_point_2::Array,
    terminal_point_2::Array;
    label1::Any = "",
    label2::Any = "",
    style1::String = "solid",
    style2::String = "solid",
    offset1::Real = 10,
    offset2::Real = 10,
    size::Real = 16,
    padding::Real = 0.15,
)

    Set([initial_point_1, terminal_point_1]) != Set([initial_point_2, terminal_point_2]) &&
        @error "Only two different points allowed"

    padded_is = []
    padded_ts = []
    angles = []
    for (i, t) in [[initial_point_1, terminal_point_1], [initial_point_2, terminal_point_2]]
         v⃗ = t .- i
           θ = atan(v⃗[2] / v⃗[1]) + (v⃗[1] < 0 ? π : 0)

        padded_i = between(cart2luxor(i), cart2luxor(t), padding)
        padded_t = between(cart2luxor(i), cart2luxor(t), 1 - padding)

         v⃗ₚ = luxor2cart(padded_i) .- luxor2cart(padded_t)
        x = 1
        y = 1
          if (-1 * v⃗ₚ[1] * x) / v⃗ₚ[2] |> isinf
              x = (v⃗ₚ[2] * y) / (-1 * v⃗ₚ[1])
        else
              y = (-1 * v⃗ₚ[1] * x) / v⃗ₚ[2]
        end
        i_new = parametric_point(
            length(padded_is) == 0 ? 10 : -10,
            luxor2cart(padded_i),
            [x; y],
        )
        t_new = parametric_point(
            length(padded_ts) == 0 ? 10 : -10,
            luxor2cart(padded_t),
            [x; y],
        )
        push!(padded_is, i_new)
        push!(padded_ts, t_new)
        push!(angles, θ)
    end

    setdash(style1)

    arrow(
        cart2luxor(padded_is[1]),
        cart2luxor(padded_ts[1]),
        arrowheadlength = 15,
        decoration = 0.5,
        decorate = () ->
            morphtext(label1, offset = offset1 + 20, angle = angles[1], size = size),
        linewidth = 3,
    )

    setdash(style2)

    arrow(
        cart2luxor(padded_is[2]),
        cart2luxor(padded_ts[2]),
        arrowheadlength = 15,
        decoration = 0.5,
        decorate = () ->
            morphtext(label2, offset = -offset2 - 10, angle = angles[2], size = size),
        linewidth = 3,
    )

    setdash("solid")

end

"""
Draws a Commutative Diagram of a given graph
"""
function commutative_diagram(C; offset::Real = 30, size::Real = 20, padding::Real = 0.15)

    fontsize(20)
    sethue("white")
    for i in 1:nparts(C, :V)
	p = C[i, :loc]
	text(C[i, :v_labels], cart2luxor(p), valign = :middle, halign = :center)
    end

    for i in 1:nparts(C, :E)
        initial_point = C[i, :src] |> s -> C[s, :loc]
	terminal_point = C[i, :tgt] |> s -> C[s, :loc]
	 v⃗ = terminal_point .- initial_point
	   θ = atan(v⃗[2] / v⃗[1]) + (v⃗[1] < 0 ? π : 0)

	padded_initial_point =
	    between(cart2luxor(initial_point), cart2luxor(terminal_point), padding)
	padded_terminal_point =
	    between(cart2luxor(initial_point), cart2luxor(terminal_point), 1 - padding)

	arrow(
	    padded_initial_point,
	    padded_terminal_point,
	    arrowheadlength = 15,
	    decoration = 0.5,
	    decorate = () -> morphtext(C[i, :e_labels], offset = offset, angle = θ, size = size),
	    linewidth = 3,
	)
    end

end

#############
# CONSTANTS #
#############

width = 500
height = 500
path = "catlab_integration.png"
color = RGBA(0.0, 0.0, 0.0)
op = Point(width / 2, height / 2)

p1 = [-125, 125] # Top left quadrant
p2 = [125, 125] # Top right quadrant
p3 = [-125, -125] # Bottom left quadrant
p4 = [125, -125] # Bottom right quadrant

########
# MAIN #
########

my_draw = make_drawing(width, height, path, color, op)

@present TheoryGraph(FreeSchema) begin
    V::Ob
    E::Ob
    src::Hom(E, V)
    tgt::Hom(E, V)
    T::AttrType
    loc::Attr(V, T)
    v_labels::Attr(V, T)
    e_labels::Attr(E, T)
end

@acset_type Graph(TheoryGraph, index = [:src, :tgt])

g = @acset Graph{Any} begin
    V = 4
    E = 5
    src = [1, 1, 2, 2, 3]
    tgt = [2, 3, 3, 4, 4]
    loc = [p1, p2, p3, p4]
    v_labels = [L"A", L"A", L"B", L"B"]
    e_labels = [L"id_{A}", L"f", L"f", L"f", L"id_{B}"]
end

commutative_diagram(g; padding = .1)

finish()


text:
## lines for reproducibility

using Pkg
Pkg.activate(".")


text:
using Test, SparseArrays, SHA
using LDPCStorage

qc_file_path = "$(pkgdir(LDPCStorage))/test/files/test_Hqc.cscmat"
h_file_path = "$(pkgdir(LDPCStorage))/test/files/test_H.cscmat"

"""
Returns a 256 bit hash of an sparse matrix.
This function should only be used for unit tests!
"""
function hash_sparse_matrix(H::SparseMatrixCSC)
    ctx = SHA2_256_CTX()

    update!(ctx, H.colptr .% UInt8)
    update!(ctx, H.rowval .% UInt8)
    update!(ctx, H.nzval .% UInt8)

    return digest!(ctx)
end

@testset "hash_sparse_matrix" begin
    hash = hash_sparse_matrix(sparse([1 2 3 0; 5 0 7 8]))
    @test hash == UInt8[0x82, 0x80, 0xcf, 0xe6, 0x13, 0x4b, 0x54, 0x50,
        0x2a, 0xc4, 0x21, 0x33, 0x6e, 0xd5, 0xe4, 0x2d,
        0x61, 0xa9, 0x30, 0x45, 0x77, 0xa7, 0x45, 0x09,
        0x29, 0xca, 0x0b, 0x5a, 0x5b, 0x9f, 0x0b, 0x0f]
end

@testset "load_cscmat" begin
    Hqc = load_cscmat(qc_file_path; print_file_header=false)
    lifting_factor = 32
    @test hash_sparse_matrix(Hqc) == UInt8[0xf7, 0x41, 0xee, 0x9c, 0xf8, 0x2c, 0x0b, 0xb8, 0x02, 0x2f, 0x4f, 0xa6, 0x40, 0xa5, 0xb0,
        0xcd, 0xf0, 0x7c, 0x11, 0xfc, 0x6d, 0x99, 0xb0, 0x13, 0x9e, 0x18, 0x19, 0xb2, 0xc7, 0x7d, 0x25, 0x92]

    H = load_cscmat(h_file_path; print_file_header=false)
    @test hash_sparse_matrix(H) == UInt8[0xf0, 0xe7, 0x55, 0xa8, 0x38, 0x9c, 0xf6, 0x1f, 0x41, 0x5e, 0x6d, 0xd4, 0xaa, 0x38, 0xee,
        0xb6, 0x77, 0x05, 0xad, 0x73, 0x32, 0xc6, 0xf0, 0x0e, 0x75, 0x2a, 0x1d, 0xe3, 0xc9, 0x88, 0x67, 0xc5]
end

@testset "save_to_cscmat and load_cscmat" begin
    Hqc = load_cscmat(qc_file_path; print_file_header=false)
    lifting_factor = 32

    H = load_cscmat(h_file_path; print_file_header=false)

    for allow_omit_entries_if_only_stored_ones in [false, true]
        for try_hex in [false, true]
            for additional_header_lines in ["", "QC matrix with expansion factor 32", "QC matrix\nwith expansion factor\n\n# 32"]

                target_file = tempname() * ".cscmat"

                @show allow_omit_entries_if_only_stored_ones
                @show try_hex

                save_to_cscmat(
                    H, target_file;
                    allow_omit_entries_if_only_stored_ones, try_hex, additional_header_lines)
                H_read = load_cscmat(target_file; print_file_header=false)
                @test H_read == H

                save_to_cscmat(
                    Hqc, target_file;
                    allow_omit_entries_if_only_stored_ones, try_hex, additional_header_lines)
                H_read = load_cscmat(target_file; print_file_header=true)
                @test H_read == Hqc
                println("\n")
            end
        end
    end

end

@testset "Hqc_to_pcm" begin
    Hqc = load_cscmat(qc_file_path; print_file_header=false)
    lifting_factor = 32

    H = load_cscmat(h_file_path; print_file_header=false)

    @test LDPCStorage.Hqc_to_pcm(Hqc, lifting_factor) == H
end

@testset "load_matrix_from_qc_cscmat_file" begin
    H_32 = load_matrix_from_qc_cscmat_file(qc_file_path; expansion_factor=32)
    H_auto = load_matrix_from_qc_cscmat_file(qc_file_path)
    H_true = load_cscmat(h_file_path; print_file_header=false)

    @test H_true == H_32
    @test H_true == H_auto
end