# Pojedyncze wyporzyczenie i oddanie roweru
function single_sim(wynik_sim=wynik_sim_init,
                    prawdopodobienstwo_stacji= prawdopodobienstwo_stacji)
    stacja_odebrania = rand(Categorical(prawdopodobienstwo_stacji))
    stacja_oddania = rand(Categorical(prawdopodobienstwo_stacji))
    while (stacja_odebrania == stacja_oddania)
        stacja_odebrania = rand(Categorical(prawdopodobienstwo_stacji))
        stacja_oddania = rand(Categorical(prawdopodobienstwo_stacji)) 
    end
    if wynik_sim[1][stacja_odebrania] >= 1
        wynik_sim[1][stacja_odebrania] -= 1
        wynik_sim[1][stacja_oddania] += 1
    else
        wynik_sim[2][stacja_odebrania] += 1
    end
    return wynik_sim
end

# x razy single_sim
function run_day(liczba_dzienna,
                    wynik_sim = wynik_sim_init,
                    prawdopodobienstwo_stacji=prawdopodobienstwo_stacji)
    wynik_sim_init = deepcopy(wynik_sim)
    for i in range(1,liczba_dzienna)
        wynik_sim = single_sim(wynik_sim_init, prawdopodobienstwo_stacji)
    end
    return wynik_sim
end

# x razy single_sim, ale zapamietuje wyniki z konca symulacji
function run_day_cum(liczba_dzienna,
                    wynik_sim = wynik_sim_init,
                    prawdopodobienstwo_stacji=prawdopodobienstwo_stacji)
    for i in range(1,liczba_dzienna)
        wynik_sim = single_sim(wynik_sim, prawdopodobienstwo_stacji)
    end
    return wynik_sim
end

# x razy single_sim, ale kazdy "dzien" liczony od nowa
function run_sims(liczba_sim = 30,
                    liczba_dzienna = 40,
                    wynik_sim = wynik_sim_init,
                    prawdopodobienstwo_stacji = prawdopodobienstwo_stacji)
    wynik_sim_init = deepcopy(wynik_sim)
    wyniki_okresu = DataFrame(dzien = Int[], nr_stacji = Int[], liczba_na_koniec_dnia = Int[], braki_na_koniec_dnia = Int[])
    for i in range(1,liczba_sim)
        wynik_dnia = run_day(liczba_dzienna, wynik_sim, prawdopodobienstwo_stacji)
        temp_df = DataFrame(dzien = [i,i,i,i,i], nr_stacji = repeat(1:5,1),liczba_na_koniec_dnia = wynik_dnia[1], braki_na_koniec_dnia = wynik_dnia[2])
        wyniki_okresu = vcat(wyniki_okresu, temp_df)
        wynik_sim = wynik_sim_init
    end
    wynik_sim = wynik_sim_init
    return wyniki_okresu
end

# x razy single_sim, ale dni są "połączone" ze sobą
function run_sims_cum(liczba_sim = 30,
                    liczba_dzienna = 40,
                    wynik_sim = wynik_sim_init,
                    prawdopodobienstwo_stacji = prawdopodobienstwo_stacji)
    wyniki_okresu = DataFrame(dzien = Int[], nr_stacji = Int[], liczba_na_koniec_dnia = Int[], braki_na_koniec_dnia = Int[])
    for i in range(1,liczba_sim)
        wynik_dnia = run_day_cum(liczba_dzienna, wynik_sim_init, prawdopodobienstwo_stacji)
        temp_df = DataFrame(dzien = [i,i,i,i,i], nr_stacji = repeat(1:5,1),liczba_na_koniec_dnia = wynik_dnia[1], braki_na_koniec_dnia = wynik_dnia[2])
        wyniki_okresu = vcat(wyniki_okresu, temp_df)
    end
    wynik_sim = wynik_sim_init
    return wyniki_okresu
end

# Wyniki skumulowanych wyników
liczba_rowerow_na_stacji = [5,5,5,5,5]
brak_roweru = [0,0,0,0,0]
prawdopodobienstwo_stacji = [0.3, 0.3, 0.2, 0.1, 0.1]
wynik_sim_init = [deepcopy(liczba_rowerow_na_stacji), deepcopy(brak_roweru)]
wyniki_okresu_cum = run_sims_cum()
# Wyniki nieskumulowanych wyników
liczba_rowerow_na_stacji = [5,5,5,5,5]
brak_roweru = [0,0,0,0,0]
prawdopodobienstwo_stacji = [0.3, 0.3, 0.2, 0.1, 0.1]
wynik_sim_init = [deepcopy(liczba_rowerow_na_stacji), deepcopy(brak_roweru)]
wyniki_okresu = run_sims()
# wykres 1
wyfiltrowane = filter(row -> row.braki_na_koniec_dnia > 0, wyniki_okresu)
bar(wyfiltrowane[:,1], wyfiltrowane[:,4], alpha = 0.7, group = wyfiltrowane[:,2])
title!("Braki rowerów na stacjach w poszczególnych dniach")
xlabel!("Dzień")
ylabel!("Niezrealizowane zamówienie")
# suma brakow z NIESKUMULOWANYCH symulacji
braki_30_dni = combine(groupby(wyfiltrowane, :nr_stacji), :braki_na_koniec_dnia => sum)
# wykres 2
bar(braki_30_dni[:,1], braki_30_dni[:,2])
title!("Skumulowana liczba braków realizacji zamówienia")
xaxis!("Numer stacji")
print("Zrealizowane zamówienia : ", round(sum(braki_30_dni[:,2])/(40*30), digits=3), "%") # ilosc dni * srednia dzienna

There might be some flows in this Julia simulation of bike stations. Detect them.