text:
"""
The pycity_scheduling framework

Copyright (C) 2022,
Institute for Automation of Complex Power Systems (ACS),
E.ON Energy Research Center (E.ON ERC),
RWTH Aachen University

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
documentation files (the "Software"), to deal in the Software without restriction, including without limitation the
rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit
persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the
Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec

import pycity_scheduling.util.factory as factory
from pycity_scheduling.classes import *
from pycity_scheduling.algorithms import *

# This is a simple power scheduling example to demonstrate the integration and impact of building cooling and
# heating loads using the central optimization algorithm.

def main(do_plot=False):
    print("\n\n------ Example 18: Scheduling Heating and Cooling Loads ------\n\n")

    # Scheduling will be performed for a full year:
    env = factory.generate_standard_environment(step_size=3600, op_horizon=24*365, mpc_horizon=None,
                                                mpc_step_width=None, initial_date=(2018, 1, 1), initial_time=(0, 0, 0))

    # Use a standardized thermal load profile for space heating and then 'convert' it into a cooling load
    # (inverted load assumption according to pycity_base):
    sh_slp = SpaceHeating(environment=env, method=1, living_area=150, specific_demand=50)

    sc_slp = SpaceCooling(environment=env,
                          method=0,
                          loadcurve=np.ones(env.timer.timesteps_total)*max(sh_slp.get_power()) / 1000.0 - sh_slp.
                          get_power() / 1000.0)

    # Plot the thermal loads:
    plot_time = list(range(env.timer.timesteps_used_horizon))

    gs = gridspec.GridSpec(2, 1)
    ax0 = plt.subplot(gs[0])
    ax0.plot(plot_time, sh_slp.p_th_heat_schedule)
    plt.xlim((0, env.timer.timesteps_used_horizon - 1))
    plt.ylabel("Space Heating [kW]")

    ax1 = plt.subplot(gs[1], sharex=ax0)
    ax1.plot(plot_time, sc_slp.p_th_cool_schedule)
    plt.xlim((0, env.timer.timesteps_used_horizon - 1))
    plt.grid()
    plt.ylabel("Space Cooling [kW]")
    plt.title("Thermal Load Profiles")
    plt.grid()
    if do_plot:
        figManager = plt.get_current_fig_manager()
        if hasattr(figManager, "window"):
            figManagerWindow = figManager.window
            if hasattr(figManagerWindow, "state"):
                figManager.window.state("zoomed")
        plt.show()

    # Now perform the scheduling of two buildings, where the first building is equipped with the heating load and
    # a thermal heating system and the second building is equipped with the cooling load and a thermal cooling system.

    # City district with district operator objective "price":
    cd = CityDistrict(environment=env, objective='max-consumption')

    # Building no. one with building objective price:
    bd1 = Building(environment=env, objective='price')
    cd.addEntity(entity=bd1, position=[0, 0])
    bes = BuildingEnergySystem(environment=env)
    bd1.addEntity(bes)
    ths = ThermalHeatingStorage(environment=env, e_th_max=16.0, soc_init=0.5, loss_factor=0.05)
    bes.addDevice(ths)
    hp = HeatPump(environment=env, p_th_nom=8.0)
    bes.addDevice(hp)
    ap = Apartment(environment=env)
    bd1.addEntity(ap)
    ap.addEntity(sh_slp)

    # Building no. two with building objective co2:
    bd2 = Building(environment=env, objective='co2')
    cd.addEntity(entity=bd2, position=[0, 1])
    bes = BuildingEnergySystem(environment=env)
    bd2.addEntity(bes)
    tcs = ThermalCoolingStorage(environment=env, e_th_max=10.0, soc_init=0.5, loss_factor=0.1)
    bes.addDevice(tcs)
    ch = Chiller(environment=env, p_th_nom=5.0, cop=3.5)
    bes.addDevice(ch)
    ap = Apartment(environment=env)
    bd2.addEntity(ap)
    ap.addEntity(sc_slp)

    # Perform the scheduling:
    opt = CentralOptimization(city_district=cd)
    results = opt.solve()
    cd.copy_schedule("central")

    # Plot the (thermal) schedules of interest:
    plot_time = list(range(env.timer.timesteps_used_horizon))

    gs = gridspec.GridSpec(4, 1)
    ax0 = plt.subplot(gs[0])
    ax0.plot(plot_time, -hp.p_th_heat_schedule)
    plt.xlim((0, env.timer.timesteps_used_horizon - 1))
    plt.ylim([-15, 15])
    plt.grid()
    plt.title("Thermal Schedules")
    plt.ylabel("Heat Pump [kW]")

    ax1 = plt.subplot(gs[1], sharex=ax0)
    ax1.plot(plot_time, ths.e_th_heat_schedule)
    plt.xlim((0, env.timer.timesteps_used_horizon - 1))
    plt.grid()
    plt.ylabel("Heating Storage [kWh]")

    ax2 = plt.subplot(gs[2], sharex=ax0)
    ax2.plot(plot_time, -ch.p_th_cool_schedule)
    plt.xlim((0, env.timer.timesteps_used_horizon - 1))
    plt.grid()
    plt.ylabel("Chiller [kW]")

    ax3 = plt.subplot(gs[3], sharex=ax0)
    ax3.plot(plot_time, tcs.e_th_cool_schedule)
    plt.xlim((0, env.timer.timesteps_used_horizon - 1))
    plt.grid()
    plt.ylabel("Cooling Storage [kWh]")

    plt.xlabel("Time", fontsize=12)

    if do_plot:
        figManager = plt.get_current_fig_manager()
        if hasattr(figManager, "window"):
            figManagerWindow = figManager.window
            if hasattr(figManagerWindow, "state"):
                figManager.window.state("zoomed")
        plt.show()
    return

# Conclusions:
# In analogy to the previous example (compare example_16_scheduling_pv+battery_system.py), also thermal energy storage
# units can be used as a buffer inside a sector-coupled energy system to provide flexibility to the local electrical
# system. In this example, for instance, the thermal heating and cooling storage units are used to (partly) decouple the
# heating respectively cooling system's thermal generation from its electrical consumption. Thus, it becomes evident
# that one can take advantage from the thermal flexibility to achieve the given local and system level scheduling
# objectives, too.

if __name__ == '__main__':
    # Run example:
    main(do_plot=True)


text:
#!/usr/bin/env python
# coding: utf-8

# # Table of Contents
#

1  Oraux CentraleSupélec PSI - Juin 2019

1.1  Remarques préliminaires

1.1.1  Importez les modules :

1.2  Planche 160

1.3  Planche 162

1.4  Planche 166

1.5  Planche 168

1.6  Planche 170

1.7  Planche 172

1.8  Planche 177

2  À voir aussi

2.1  Les oraux   (exercices de maths avec Python)

2.2  Fiches de révisions pour les oraux

2.3  D'autres notebooks ?

# # Oraux CentraleSupélec PSI - Juin 2019
#
# - Ce [notebook Jupyter](https://www.jupyter.org) est une proposition de correction, en [Python 3](https://www.python.org/), d'exercices d'annales de l'épreuve "maths-info" du [concours CentraleSupélec](http://www.concours-centrale-supelec.fr/), filière PSI.
# - Les exercices viennent de l'[Officiel de la Taupe](http://odlt.fr/), [2017](http://www.odlt.fr/Oraux_2017.pdf) (planches 162 à 177, pages 10 et 11).
# - Ce document a été écrit par [Lilian Besson](http://perso.crans.org/besson/), et est disponible en ligne [sur mon site](https://perso.crans.org/besson/publis/notebooks/Oraux_CentraleSupelec_PSI__Juin_2019.html).

# ## Remarques préliminaires
# - Les exercices sans Python ne sont pas traités.
# - Les exercices avec Python utilisent Python 3, [numpy](http://numpy.org), [matplotlib](http://matplotlib.org), [scipy](http://scipy.org) et [sympy](http://sympy.org), et essaient d'être résolus le plus simplement et le plus rapidement possible. L'efficacité (algorithmique, en terme de mémoire et de temps de calcul), n'est *pas* une priorité. La concision et simplicité de la solution proposée est prioritaire.

# In[1]:

import numpy as np
import matplotlib as mpl
mpl.rcParams['figure.figsize'] = (10, 7)
import matplotlib.pyplot as plt
from scipy import integrate
import numpy.random as rd

# In[16]:

import seaborn as sns
sns.set(context="notebook", style="whitegrid", palette="hls", font="sans-serif", font_scale=1.1)

# ### Importez les modules :

# In[ ]:

import numpy as np

# In[ ]:

import matplotlib.pyplot as plt

# In[ ]:

from scipy import integrate

# In[ ]:

import numpy.random as rd

# ----
# ## Planche 160

# - $I_n := \int_0^1 \frac{1}{(1+t)^n \sqrt{1-t}} \mathrm{d}t$ et $I_n := \int_0^1 \frac{1/2}{(1+t)^n \sqrt{1-t}} \mathrm{d}t$ sont définies pour tout $n$ car leur intégrande est continue et bien définie sur $]0,1[$ et intégrable en $1$ parce qu'on sait (par intégrale de Riemann) que $\frac{1}{\sqrt{u}}$ est intégrable en $0^+$ (et changement de variable $u = 1-t$).
# - On les calcule très simplement :

# In[2]:

def I(n):
    def f(t):
        return 1 / ((1+t)**n * np.sqrt(1-t))
    i, err = integrate.quad(f, 0, 1)
    return i

# In[3]:

def J(n):
    def f(t):
        return 1 / ((1+t)**n * np.sqrt(1-t))
    i, err = integrate.quad(f, 0, 0.5)
    return i

# In[4]:

valeurs_n = np.arange(1, 50)
valeurs_In = np.array([I(n) for n in valeurs_n])

plt.figure()
plt.plot(valeurs_n, valeurs_In, 'ro')
plt.title("Valeurs de $I_n$")
plt.show()

# - On conjecture que $I_n$ est décroissante. C'est évident puisque si on note $f_n(t)$ son intégrande, on observe que $f_{n+1}(t) \leq f_n(t)$ pour tout $t$, et donc par monotonie de l'intégrale, $I_{n+1} \leq I_n$.
# - On conjecture que $I_n \to 0$. Cela se montre très facilement avec le théorème de convergence dominée.

# In[5]:

plt.figure()
plt.plot(np.log(valeurs_n), np.log(valeurs_In), 'go')
plt.title(r"Valeurs de $\ln(I_n)$ en fonction de $\ln(n)$")
plt.show()

# - Ce qu'on observe permet de conjecturer que $\alpha=1$ est l'unique entier tel que $n^{\alpha} I_n$ converge vers une limite non nulle.

# In[16]:

valeurs_n = np.arange(1, 500)
valeurs_In = np.array([I(n) for n in valeurs_n])
valeurs_Jn = np.array([J(n) for n in valeurs_n])
alpha = 0.9

plt.figure()
plt.plot(valeurs_n, valeurs_n**alpha * valeurs_In, 'r+', label=r'$n^{\alpha} I_n$')
plt.plot(valeurs_n, valeurs_n**alpha * valeurs_Jn, 'b+', label=r'$n^{\alpha} J_n$')
plt.legend()
plt.title(r"Valeurs de $n^{\alpha} I_n$ et $n^{\alpha} J_n$")
plt.show()

# - On en déduit qu'il en est de même pour $J_n$, on a $n^{\alpha} J_n \to l$ la même limite que $n^{\alpha} I_n$.

# - Pour finir, on montre mathématiquement que $n^{\alpha} (I_n - J_n)$ tend vers $0$.

# In[22]:

plt.figure()
plt.plot(valeurs_n, valeurs_n**alpha * (valeurs_In - valeurs_Jn), 'g+', label=r'$n^{\alpha} (I_n - J_n)$')
plt.legend()
plt.title(r"Valeurs de $n^{\alpha} (I_n - J_n)$")
plt.show()

# - Puis rapidement, on montre que $\forall x \geq 0, \ln(1 + x) \geq \frac{x}{1+x}$. Ca peut se prouver de plein de façons différentes, mais par exemple on écrit $f(x) = (x+1) \log(x+1) - x$ qui est de classe $\mathcal{C}^1$, et on la dérive. $f'(x) = \log(x+1) + 1 - 1 > 0$ donc $f$ est croissante, et $f(0) = 0$ donc $f(x) \geq f(0) = 0$ pour tout $x \geq 0$.

# In[26]:

X = np.linspace(0, 100, 1000)
plt.plot(X, np.log(1 + X), 'ro-', label=r'$\log(1+x)$', markevery=50)
plt.plot(X, X / (1 + X), 'b+-', label=r'$\frac{x}{1+x}$', markevery=50)
plt.legend()
plt.title("Comparaison entre deux fonctions")
plt.show()

# ---
# ## Planche 162

# On commence par définir la fonction, en utilisant `numpy.cos` et pas `math.cos` (les fonctions de `numpy` peuvent travailler sur des tableaux, c'est plus pratique).

# In[27]:

def f(x):
    return x * (1 - x) * (1 + np.cos(5 * np.pi * x))

Xs = np.linspace(0, 1, 2000)
Ys = f(Xs)

# Pas besoin de lire le maximum sur un graphique :

# In[29]:

M = max_de_f = max(Ys)
print("Sur [0, 1], avec 2000 points, le maximum de f(x) est M =", M)

# In[32]:

i_maximisant_f = np.argmax(Ys)
x_maximisant_f = Xs[i_maximisant_f]
print("Sur ces 2000 points, le maximum est atteint en x =", x_maximisant_f)

# On affiche la fonction, comme demandé, avec un titre :

# In[39]:

plt.figure()
plt.plot(Xs, Ys)
plt.vlines(x_maximisant_f, min(Ys), max(Ys), color="b", linestyles="dotted")
plt.hlines(max(Ys), min(Xs), max(Xs), color="b", linestyles="dotted")
plt.title("Fonction $f(x)$ sur $[0,1]$")
plt.show()

# Pour calculer l'intégrale, on utilise `scipy.integrate.quad` :

# In[40]:

def In(x, n):
    def fn(x):
        return f(x) ** n
    return integrate.quad(fn, 0, 1)[0]

def Sn(x):
    return np.sum([In(Xs, n) * x**n for n in range(0, n+1)], axis=0)

# On vérifie avant de se lancer dans l'affichage :

# In[41]:

for n in range(10):
    print("In(x,", n, ") =", In(Xs, n))

# In[45]:

a = 1/M + 0.1
X2s = np.linspace(-a, a, 2000)

plt.figure()
for n in [10, 20, 30, 40, 50]:
    plt.plot(X2s, Sn(X2s), "-+", label="n =" + str(n), markevery=20)
plt.legend()
plt.show()

# $S_n(x)$ semble diverger pour $x\to2^-$ quand $n\to\infty$.
# Le rayon de convergence de la série $\sum In x^n$ **semble** être $2$.

# In[46]:

def un(n):
    return In(Xs, n + 1) / In(Xs, n)

# In[47]:

for n in range(10):
    print("un =", un(n), "pour n =", n)

# Ici, `un` ne peut pas être utilisé comme une fonction "numpy" qui travaille sur un tableau, on stocke donc les valeurs "plus manuellement" :

# In[50]:

def affiche_termes_un(N):
    valeurs_un = [0] * N
    for n in range(N):
        valeurs_un[n] = un(n)

    plt.figure()
    plt.plot(valeurs_un, 'o-')
    plt.title("Suite $u_n$")
    plt.grid(True)
    plt.show()

# In[51]:

affiche_termes_un(30)

# La suite $u_n$ semble être croissante (on peut le prouver), toujours plus petite que $1$ (se prouve facilement aussi, $I_{n+1} < I_n$), et semble converger.
# Peut-être vers $1/2$, il faut aller regarder plus loin ?

# In[52]:

affiche_termes_un(100)

# Pour conclure, on peut prouver que la suite est monotone et bornée, donc elle converge.
# Il est plus dur de calculer sa limite, et cela sort de l'exercice.

# ---
# ## Planche 166

# In[53]:

case_max = 12
univers = list(range(case_max))

# In[54]:

def prochaine_case(case):
    return (case + rd.randint(1, 6+1)) % case_max

# In[55]:

def Yn(duree, depart=0):
    case = depart
    for coup in range(duree):
        case = prochaine_case(case)
    return case

# Avant de s'en servir pour simuler plein de trajectoirs, on peut vérifier :
#
# - en un coup, on avance pas plus de 6 cases :

# In[56]:

[Yn(1) for _ in range(10)]

# - En 100 coups, on commence à ne plus voir de tendance :

# In[62]:

[Yn(100) for _ in range(10)]

# Pour l'histogramme, on triche un peu en utilisant `numpy.bincount`. Mais on peut le faire à la main très simplement !

# In[64]:

observations = [Yn(100) for _ in range(10)]
print(observations)
np.bincount(observations, minlength=case_max)

# In[65]:

def histogramme(duree, repetitions=5000):
    cases = [Yn(duree) for _ in range(repetitions)]
    frequences = np.bincount(cases, minlength=case_max)
    # aussi a la main si besoin
    frequences = [0] * case_max
    for case in cases:
        frequences[case] += 1
    return frequences / np.sum(frequences)

# In[66]:

histogramme(50)

# On va afficher des histogrammes :

# In[67]:

def voir_histogramme(valeurs_n):
    for n in valeurs_n:
        plt.figure()
        plt.bar(np.arange(case_max), histogramme(n))
        plt.title("Histogramme de cases visitées en " + str(n) + " coups")
        plt.show()

# In[68]:

voir_histogramme([1, 2, 3, 50, 100, 200])

# On s'approche d'une distribution uniforme !

# On a tout simplement l'expression suivante :
# $$\forall n \geq 0, \mathbb{P}(Y_{n+1} = k) = \frac{1}{6} \sum_{\delta = 1}^{6} \mathbb{P}(Y_n = k - \delta \mod 12).$$
# Avec $k - 1 \mod 12 = 11$ si $k = 0$ par exemple.

# On a donc la matrice suivante pour exprimer $U_n = (\mathbb{P}(Y_n = k))_{0\leq k \leq 11}$ en fonction de $U_{n-1}$ :
#
# $$ P = \frac{1}{6} \begin{bmatrix}
# 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 \\
# 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1\\
# 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\\
# 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1\\
# 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1\\
# 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\
# 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
# 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\
# 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\
# 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 \\
# 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 \\
# 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 \\
# \end{bmatrix}$$

# On va la définir rapidement en Python, et calculer ses valeurs propres notamment.

# In[18]:

case_max = 12

# In[19]:

P = np.zeros((case_max, case_max))

# In[26]:

for k in range(case_max):
    for i in range(k - 6, k):
        P[k, i] = 1/6

# In[27]:

P

# In[28]:

import numpy.linalg as LA

# In[29]:

spectre, vecteur_propres = LA.eig(P)

# On a besoin d'éliminer les erreurs d'arrondis, mais on voit que $1$ est valeur propre, associée au vecteur $[1,\dots,1,\dots,1]$ avec un $1$ seulement à la 8ème composante.

# In[34]:

np.round(spectre,10)

# In[31]:

np.round(vecteur_propres[0])

# $P$ n'est pas diagonalisable, **à prouver** au tableau si l'examinateur le demande.

# ----
# ## Planche 168

# - Soit $f(x) = \frac{1}{2 - \exp(x)}$, et $a(n) = \frac{f^{(n)}(0)}{n!}$.
#

# In[76]:

def f(x):
    return 1 / (2 - np.exp(x))

# - Soit $g(x) = 2 - \exp(x)$, telle que $g(x) f(x) = 1$. En dérivant $n > 0$ fois cette identité et en utilisant la formule de Leibniz, on trouve :
# $$ (g(x)f(x))^{(n)} = 0 = \sum_{k=0}^n {n \choose k} g^{(k)}(x) f^{(n-k)}(x).$$
# Donc en $x=0$, on utilise que $g^{(k)}(x) = - \exp(x)$, qui donne que $g^{(k)}(0) = 1$ si $k=0$ ou $-1$ sinon, pour trouver que $\sum_{k=0}^n {n \choose k} f^{(k)}(0) = f^{(n)}(0)$. En écrivant ${n \choose k} = \frac{k! (n-k)!}{n!}$ et avec la formule définissant $a(n)$, cela donne directement la somme recherchée : $$ a(n) = \sum_{k=1}^n \frac{a(n-k)}{k!}.$$

# - Pour calculer $a(n)$ avec Python, on utilise cette formule comme une formule récursive, et on triche un peu en utilisant `math.factorial` pour calculer $k!$.Il nous faut aussi $a(0) = f(0) = 1$ :

# In[77]:

from math import factorial

def a_0an(nMax):
    valeurs_a = np.zeros(nMax+1)
    valeurs_a[0] = 1.0
    for n in range(1, nMax+1):
        valeurs_a[n] = sum(valeurs_a[n-k] / factorial(k) for k in range(1, n+1))
    return valeurs_a

# In[78]:

nMax = 10
valeurs_n = np.arange(0, nMax + 1)
valeurs_a = a_0an(nMax)

for n in valeurs_n:
    print("Pour n =", n, "on a a(n) =", valeurs_a[n])

# In[83]:

plt.figure()
plt.plot(valeurs_n, valeurs_a, 'ro-', label=r'$a(n)$', markersize=12)
plt.plot(valeurs_n, 1 / np.log(2)**valeurs_n, 'g+-', label=r'$1/\log(2)^n$')
plt.plot(valeurs_n, 1 / (2 * np.log(2)**valeurs_n), 'bd-', label=r'$1/(2\log(2)^n)$')
plt.title("$a(n)$ et deux autres suites")
plt.legend()
plt.show()

# - On observe que $a(n)$ est comprise entre $\frac{1}{2(\log(2))^n}$ et $\frac{1}{\log(2)^n}$, donc le rayon de convergence de $S(x) = \sum a(n) x^n$ est $\log(2)$.
# - On va calculer les sommes partielles $S_n(x)$ de la série $S(x)$ :

# In[84]:

def Sn(x, n):
    valeurs_a = a_0an(n)
    return sum(valeurs_a[k] * x**k for k in range(0, n + 1))

# On peut vérifie que notre fonction marche :

# In[85]:

x = 0.5
for n in range(0, 6 + 1):
    print("Pour n =", n, "S_n(x) =", Sn(x, n))

# In[86]:

valeurs_x = np.linspace(0, 0.5, 1000)
valeurs_f = f(valeurs_x)

# Je pense que l'énoncé comporte une typo sur l'intervale ! Vu le rayon de convergence, on ne voit rien si on affiche sur $[0,10]$ !

# In[89]:

plt.figure()
for n in range(0, 6 + 1):
    valeurs_Sn = []
    for x in valeurs_x:
        valeurs_Sn.append(Sn(x, n))
    plt.plot(valeurs_x, valeurs_Sn, 'o:', label='$S_' + str(n) + '(x)$', markevery=50)
plt.plot(valeurs_x, valeurs_f, '+-', label='$f(x)$', markevery=50)
plt.title("$f(x)$ et $S_n(x)$ pour $n = 0$ à $n = 6$")
plt.legend()
plt.show()

# ## Planche 170

# In[90]:

def u(n):
    return np.arctan(n+1) - np.arctan(n)

# In[92]:

valeurs_n = np.arange(50)
valeurs_u = u(valeurs_n)

plt.figure()
plt.plot(valeurs_n, valeurs_u, "o-")
plt.xlabel("Entier $n$")
plt.ylabel("Valeur $u_n$")
plt.title("Premières valeurs de $u_n$")

# On peut vérifier le prognostic quand à la somme de la série $\sum u_n$ :

# In[93]:

pi/2

# In[94]:

sum(valeurs_u)

# In[95]:

somme_serie = pi/2
somme_partielle = sum(valeurs_u)
erreur_relative = abs(somme_partielle - somme_serie) / somme_serie
erreur_relative

# Avec seulement $50$ termes, on a moins de $1.5\%$ d'erreur relative, c'est déjà pas mal !

# $(u_n)_n$ semble être décroisante, et tendre vers $0$. On peut prouver ça mathématiquement.

# On sait aussi que $\forall x\neq0, \arctan(x) + \arctan(1/x) = \frac{\pi}{2}$, et que $\arctan(x) \sim x$, donc on obtient que $u_n \sim \frac{1}{n} - \frac{1}{n+1} = \frac{1}{n(n+1)}$.
# On peut le vérifier :

# In[97]:

valeurs_n = np.arange(10, 1000)
valeurs_u = u(valeurs_n)
valeurs_equivalents = 1 / (valeurs_n * (valeurs_n + 1))

plt.figure()
plt.plot(valeurs_n, valeurs_u / valeurs_equivalents, "-")
plt.xlabel("Entier $n$")
plt.title(r"Valeurs de $u_n / \frac{1}{n(n+1)}$")

# - Pour $e = (e_n)_{n\in\mathbb{N}}$ une suite de nombres égaux à $0$ ou $1$ (*i.e.*, $\forall n, e_n \in \{0,1\}$, $S_n(e) = \sum{i=0}^n e_i u_i$ est bornée entre $0$ et $\sum_{i=0}^n u_i$. Et $u_n \sim \frac{1}{n(n+1)}$ qui est le terme général d'une série convergente (par critère de Cauchy, par exemple, avec $\alpha=2$). Donc la série $\sum u_n$ converge et donc par encadrement, $S_n(e)$ converge pour $n\to\infty$, *i.e.*, $S(e)$ converge. Ces justifications donnent aussi que $$0 \leq S(e) \leq \sum_{n\geq0} u_n = \lim_{n\to\infty} \arctan(n) - \arctan(0) = \frac{\pi}{2}.$$

# - Pour $e = (0, 1, 0, 1, \ldots)$, $S(e)$ peut être calculée avec Python. Pour trouver une valeur approchée à $\delta = 10^{-5}$ près, il faut borner le **reste** de la série, $R_n(e) = \sum_{i \geq n + 1} e_i u_i$. Ici, $R_{2n+1}(e) \leq u_{2n+2}$ or $u_i \leq \frac{1}{i(i+1)}$, donc $R_{2n+1}(e) \leq \frac{1}{(2n+1)(2n+2)}$. $\frac{1}{(2n+1)(2n+2)} \leq \delta$ dès que $2n+1 \geq \sqrt{\delta}$, *i.e.*, $n \geq \frac{\sqrt{\delta}+1}{2}$. Calculons ça :

# In[98]:

from math import ceil, sqrt, pi

# In[99]:

def Se(e, delta=1e-5, borne_sur_n_0=10000):
    borne_sur_n_1 = int(ceil(1 + sqrt(delta)/2.0))
    borne_sur_n = max(borne_sur_n_0, borne_sur_n_1)
    somme_partielle = 0
    for n in range(0, borne_sur_n + 1):
        somme_partielle += e(n) * u(n)
    return somme_partielle

# In[100]:

def e010101(n):
    return 1 if n % 2 == 0 else 0

# In[101]:

delta = 1e-5
Se010101 = Se(e010101, delta)
print("Pour delta =", delta, "on a Se010101(delta) ~=", round(Se010101, 5))

# - Pour inverser la fonction, et trouver la suite $e$ telle que $S(e) = x$ pour un $x$ donné, il faut réfléchir un peu plus.

# In[102]:

def inverse_Se(x, n):
    assert 0 < x < pi/2.0, "Erreur : x doit être entre 0 et pi/2 strictement."
    print("Je vous laisse chercher.")
    raise NotImplementedError

# Ca suffit pour la partie Python.

# ----
# ## Planche 172

# In[103]:

from random import random

def pile(proba):
    """ True si pile, False si face (false, face, facile à retenir)."""
    return random() < proba

# - D'abord, on écrit une fonction pour **simuler** l'événement aléatoire :

# In[104]:

def En(n, p):
    lance = pile(p)
    for i in range(n - 1):
        nouveau_lance = pile(p)
        if lance and nouveau_lance:
            return False
        nouveau_lance = lance
    return True

# In[105]:

import numpy as np

# In[106]:

lances = [ En(2, 0.5) for _ in range(100) ]
np.bincount(lances)

# In[107]:

def pn(n, p, nbSimulations=100000):
    return np.mean([ En(n, p) for _ in range(nbSimulations) ])

# - Par exemple, pour seulement $2$ lancés, on a $1 - p_n = p^2$ car $\overline{E_n}$ est l'événement d'obtenir $2$ piles qui est de probabilité $p^2$.

# In[108]:

pn(2, 0.5)

# - Avec $4$ lancés, on a $p_n$ bien plus petit.

# In[109]:

pn(4, 0.5)

# - On vérifie que $p_n(n, p)$ est décroissante en $p$, à $n$ fixé :

# In[110]:

pn(4, 0.1)

# In[111]:

pn(4, 0.9)

# - On vérifie que $p_n(n, p)$ est décroissante en $n$, à $p$ fixé :

# In[112]:

pn(6, 0.2)

# In[113]:

pn(20, 0.2)

# In[114]:

pn(100, 0.2)

# - Notons que la suite semble converger ? Ou alors elle décroit de moins en moins rapidement.

# - Par récurrence et en considérant les possibles valeurs des deux derniers lancés numérotés $n+2$ et $n+1$, on peut montrer que
#   $$\forall n, p_{n+2} = (1-p) p_{n+1} + p(1-p) p_n$$

# - Si $p_n$ converge, on trouve sa limite $l$ comme point fixe de l'équation précédente. $l = (1-p) l + p(1-p) l$ ssi $1 = 1-p + p(1-p)$ ou $l=0$, donc si $p\neq0$, $l=0$. Ainsi l'événement "on obtient deux piles d'affilé sur un nombre infini de lancers$ est bien presque sûr.

# - Je vous laisse terminer pour calculer $T$ et les dernières questions.

# ----
# ## Planche 177
#
# - Le domaine de définition de $f(x) = \sum_{n \geq 1} \frac{x^n}{n^2}$ est $[-1, 1]$ car $\sum \frac{x^n}{n^k}$ converge si $\sum x^n$ converge (par $k$ dérivations successives), qui converge ssi $|x| < 1$. Et en $-1$ et $1$, on utilise $\sum \frac{1}{n^2} = \frac{\pi^2}{6}$.
#
# - Pour calculer $f(x)$ à $10^{-5}$ près, il faut calculer sa somme partielle $S_n(x) := \sum_{i=1}^n \frac{x^i}{i^2}$ en bornant son reste $S_n(x) := \sum_{i \geq n+1} \frac{x^i}{i^2}$ par (au moins) $10^{-5}$. Une inégalité montre rapidement que $R_n(x) \leq |x|^{n+1}\sum_{i\geq n+1} \frac{1}{i^2} $, et donc $R_n(x) \leq \delta$ dès que $|x|^{n+1} \leq \frac{\pi^2}{6} \delta$, puisque $\sum_{i\geq n+1} \frac{1}{i^2} \leq \sum_{i=0}^{+\infty} \frac{1}{i^2} = \frac{\pi^2}{6}$. En inversant pour trouver $n$, cela donne que le reste est contrôlé par $\delta$ dès que $n \leq \log_{|x|}\left( \frac{6}{\pi^2} \delta \right) - 1$ (si $x\neq 0$, et par $n \geq 0$ sinon).

# In[115]:

from math import floor, log, pi

# In[116]:

delta = 1e-5

def f(x):
    if x == 0: return 0
    borne_sur_n = int(floor(log((6/pi**2 * delta), abs(x)) - 1))
    somme_partielle = 0
    for n in range(1, borne_sur_n + 1):
        somme_partielle += x**n / n**2
    return somme_partielle

# In[117]:

for x in [-0.75, -0.5, 0.25, 0, 0.25, 0.5, 0.75]:
    print("Pour x =", x, "\tf(x) =", round(f(x), 5))

# - L'intégrale $g(t) = \int_0^x \frac{\ln(1 - t)}{t} \mathrm{d}t$ est bien défine sur $D = [-1, 1]$ puisque son intégrande existe, est continue et bien intégrable sur tout interval de la forme $]a, 0[$ ou $]0, b[$ pour $-1 < a < 0$ ou $0 < b < 1$. Le seul point qui peut déranger l'intégrabilité est en $0$, mais $\ln(1-t) \sim t$ quand $t\to0$ donc l'intégrande est $\sim 1$ en $0^-$ et $0^+$ et donc est bien intégrable. De plus, comme "intégrale de la borne supérieure" d'une fonction continue, $g$ est dérivable sur l'intérieur de son domaine, *i.e.*, sur $]-1, 1[$.
#
# - Pour la calculer numériquement, on utilise **évidemment** le module `scipy.integrate` et sa fonction `integrale, erreur = quad(f, a, b)`, qui donne une approximation de la valeur d'une intégrale en dimension 1 et une *borne* sur son erreur :

# In[118]:

from scipy import integrate

# In[119]:

def g(x):
    def h(t):
        return log(1 - t) / t
    integrale, erreur = integrate.quad(h, 0, x)
    return integrale

# - On visualise les deux fonctions $f$ et $g$ sur le domaine $D$ :

# In[120]:

import numpy as np
import matplotlib.pyplot as plt

# In[122]:

domaine = np.linspace(-0.99, 0.99, 1000)

valeurs_f = [f(x) for x in domaine]
valeurs_g = [g(x) for x in domaine]

plt.figure()
plt.plot(domaine, valeurs_f, "+-", label="$f(x)$", markevery=50)
plt.plot(domaine, valeurs_g, "+-", label="$g(x)$", markevery=50)
plt.legend()
plt.grid()
plt.title("Représentation de $f(x)$ et $g(x)$")
plt.show()

# - On conjecture que $g(x) = - f(x)$.
#
# La suite des questions est à faire au brouillon et sans Python :
#
# - On trouve que $f'(x) = \sum_{n\geq 1} \frac{n x^{n-1}}{n^2} = \frac{1}{x} \sum_{n\geq 1} \frac{x^n}{n}$ si $x\neq0$. Or on sait que $\log(1 + x) = \sum_{n\geq 1} \frac{x^n}{n}$ et donc cela montre bien que $g(x) = \int_0^x - f'(t) \mathrm{d}t = f(0) - f(x) = f(x)$ comme observé.
#
# - On trouve que $g(1) = - f(1) = - \frac{\pi^2}{6}$.
#
# - Par ailleurs, un changement de variable $u=1-x$ donne $g(1-x) = \int_x^1 \frac{\ln(u)}{1-u} \mathrm{d} u$, et une intégration par partie avec $a(u) = \ln(u)$ et $b'(u) = \frac{1}{1-u}$ donne $g(1-x) = [\ln(u)\ln(1-u)]_x^1 + \int_x^1 \frac{\ln(1-u)}{u} \mathrm{d}u$ et donc on reconnaît que $$g(1-x) = \ln(x)\ln(1-x) + g(1) - g(x).$$
#
# - Je vous laisse la fin comme exercice !

# ----
# # À voir aussi
#
# ## [Les oraux](http://perso.crans.org/besson/infoMP/oraux/solutions/)   *(exercices de maths avec Python)*
#
# Se préparer aux oraux de ["maths avec Python" (maths 2)](http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/#oMat2) du concours Centrale Supélec peut être utile.
#
# Après les écrits et la fin de l'année, pour ceux qui seront admissibles à Centrale-Supélec, ils vous restera les oraux (le concours Centrale-Supélec a un oral d'informatique, et un peu d'algorithmique et de Python peuvent en théorie être demandés à chaque oral de maths et de SI).
#
# Je vous invite à lire [cette page avec attention](http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/#oMat2), et à jeter un œil aux documents mis à disposition :
#
# ## Fiches de révisions *pour les oraux*
#
# 1. [Calcul matriciel](http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/Python-matrices.pdf), avec [numpy](https://docs.scipy.org/doc/numpy/) et [numpy.linalg](http://docs.scipy.org/doc/numpy/reference/routines.linalg.html),
# 2. [Réalisation de tracés](http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/Python-plot.pdf), avec [matplotlib](http://matplotlib.org/users/beginner.html),
# 3. [Analyse numérique](http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/Python-AN.pdf), avec [numpy](https://docs.scipy.org/doc/numpy/) et [scipy](http://docs.scipy.org/doc/scipy/reference/tutorial/index.html). Voir par exemple [scipy.integrate](http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html) avec les fonctions [scipy.integrate.quad](http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) (intégrale numérique) et [scipy.integrate.odeint](http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) (résolution numérique d'une équation différentielle),
# 4. [Polynômes](http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/Python-polynomes.pdf) : avec [numpy.polynomials](https://docs.scipy.org/doc/numpy/reference/routines.polynomials.package.html), [ce tutoriel peut aider](https://docs.scipy.org/doc/numpy/reference/routines.polynomials.classes.html),
# 5. [Probabilités](http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/Python-random.pdf), avec [numpy](https://docs.scipy.org/doc/numpy/) et [random](https://docs.python.org/3/library/random.html).
#
# Pour réviser : voir [ce tutoriel Matplotlib (en anglais)](http://www.labri.fr/perso/nrougier/teaching/matplotlib/), [ce tutoriel Numpy (en anglais)](http://www.labri.fr/perso/nrougier/teaching/numpy/numpy.html).
# Ainsi que tous les [TP](http://perso.crans.org/besson/infoMP/TPs/solutions/), [TD](http://perso.crans.org/besson/infoMP/TDs/solutions/) et [DS](http://perso.crans.org/besson/infoMP/DSs/solutions/) en Python que j'ai donné et corrigé au Lycée Lakanal (Sceaux, 92) en 2015-2016 !

# ----
# ## D'autres notebooks ?
#
# > Ce document est distribué [sous licence libre (MIT)](https://lbesson.mit-license.org/), comme [les autres notebooks](https://GitHub.com/Naereen/notebooks/) que j'ai écrit depuis 2015.

ext: py

reason: import numpy,from scipy


text:
# AUTOGENERATED! DO NOT EDIT! File to edit: 00_qdrant.ipynb (unless otherwise specified).

__all__ = ['collection_name', 'dim', 'prod_client', 'dev_client']

# Cell
from .tools import *
from .paths import *
from .psql import *

import pandas as pd
import numpy as np
from tqdm import tqdm

from qdrant_client import QdrantClient
import qdrant_client
collection_name = "food_images"

# Cell
dim = 768 #+ onehot.n_dim()
QdrantClient.get_collections = lambda self: [c['name'] for c in self.http.collections_api.get_collections().dict()['result']['collections']]
QdrantClient.collection_len = lambda self, name: self.http.collections_api.get_collection(name).dict()['result']['vectors_count']

# Cell
prod_client = QdrantClient(host='localhost', port=6333)
dev_client  = QdrantClient(host='localhost', port=6334)

reason: import numpy


text:
"""Communicate with the basler camera."""

import numpy as np
from pypylon import genicam, pylon

from .units import u

class BaslerCam:
    """Run a generic Basler camera, copied heavily from grab.py code online."""

    def __init__(self, px_size_um: float = 3.45, debug: bool = True) -> None:
        """Initialize the camera.

        :param px_size_um: Pixel size in um
        """
        # pixel size
        self.px_size = px_size_um * u.um

        self.debug = debug
        self.cam = pylon.InstantCamera(
            pylon.TlFactory.GetInstance().CreateFirstDevice()
        )
        self.cam.Open()

        self.exposure_time = 1000  # expsoure time in us

        self.cam.ExposureAuto.SetValue("Off")
        self.cam.ExposureTime.SetValue(self.exposure_time)

        if self.debug:
            print("Using device ", self.cam.GetDeviceInfo().GetModelName())

    def capture_image(self) -> np.array:
        """Grab an image."""
        # self.cam.Exposur
        result = self.cam.GrabOne(1000)
        return result.Array

    def find_exposure_time(self, min_peak: int = 200, max_peak: int = 255) -> None:
        """find exposure time automatically."""
        pass

if __name__ == "__main__":
    app = BaslerCam()

max_stars_repo_name: RIMS-Code/LaserProfiler

max_issues_repo_name: RIMS-Code/LaserProfiler

max_forks_repo_name: RIMS-Code/LaserProfiler


text:
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import normalize
import os
import pickle
import numpy as np

import provider

def get_gmm(points, n_gaussians, NUM_POINT, type='grid', variance=0.05, n_scales=3, D=3):
    """
    Compute weights, means and covariances for a gmm with two possible types 'grid' (2D/3D) and 'learned'

    :param points: num_points_per_model*nummodels X 3 - xyz coordinates
    :param n_gaussians: scalar of number of gaussians /  number of subdivisions for grid type
    :param NUM_POINT: number of points per model
    :param type: 'grid' / 'leared' toggle between gmm methods
    :param variance: gaussian variance for grid type gmm
    :return gmm: gmm: instance of sklearn GaussianMixture (GMM) object Gauassian mixture model
    """
    if type == 'grid':
        #Generate gaussians on a grid - supports only axis equal subdivisions
        if n_gaussians >= 32:
            print('Warning: You have set a very large number of subdivisions.')
        if not(isinstance(n_gaussians, list)):
            if D == 2:
                gmm = get_2d_grid_gmm(subdivisions=[n_gaussians, n_gaussians], variance=variance)
            elif D == 3:
                gmm = get_3d_grid_gmm(subdivisions=[n_gaussians, n_gaussians, n_gaussians], variance=variance)
            else:
                ValueError('Wrong dimension. This supports either D=2 or D=3')

    elif type == 'learn':
        #learn the gmm from given data and save to gmm.p file, if already learned then load it from existing gmm.p file for speed
        if isinstance(n_gaussians, list):
            raise  ValueError('Wrong number of gaussians: non-grid value must be a scalar')
        print("Computing GMM from data - this may take a while...")
        info_str = "g" + str(n_gaussians) + "_N" + str(len(points)) + "_M" + str(len(points) / NUM_POINT)
        gmm_dir = "gmms"
        if not os.path.exists(gmm_dir):
            os.mkdir(gmm_dir)
        filename = gmm_dir + "/gmm_" + info_str + ".p"
        if os.path.isfile(filename):
            gmm = pickle.load(open(filename, "rb"))
        else:
            gmm = get_learned_gmm(points, n_gaussians, covariance_type='diag')
            pickle.dump(gmm, open( filename, "wb"))
    else:
        ValueError('Wrong type of GMM [grid/learn]')

    return gmm

def get_learned_gmm(points, n_gaussians, covariance_type='diag'):
    """
    Learn weights, means and covariances for a gmm based on input data using sklearn EM algorithm

    :param points: num_points_per_model*nummodels X 3 - xyz coordinates
    :param n_gaussians: scalar of number of gaussians /  3 element list of number of subdivisions for grid type
    :param covariance_type: Specify the type of covariance mmatrix : 'diag', 'full','tied', 'spherical' (Note that the Fisher Vector method relies on diagonal covariance matrix)
        See sklearn documentation : http://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html
    :return gmm: gmm: instance of sklearn GaussianMixture (GMM) object Gauassian mixture model
    """
    gmm = GaussianMixture(n_components = n_gaussians, covariance_type=covariance_type)
    gmm.fit(points.astype(np.float64))
    return gmm

def get_3d_grid_gmm(subdivisions=[5,5,5], variance=0.04):
    """
    Compute the weight, mean and covariance of a gmm placed on a 3D grid
    :param subdivisions: 2 element list of number of subdivisions of the 3D space in each axes to form the grid
    :param variance: scalar for spherical gmm.p
    :return gmm: gmm: instance of sklearn GaussianMixture (GMM) object Gauassian mixture model
    """
    # n_gaussians = reduce(lambda x, y: x*y,subdivisions)
    n_gaussians = np.prod(np.array(subdivisions))
    step = [1.0/(subdivisions[0]),  1.0/(subdivisions[1]),  1.0/(subdivisions[2])]

    means = np.mgrid[ step[0]-1: 1.0-step[0]: complex(0, subdivisions[0]),
                      step[1]-1: 1.0-step[1]: complex(0, subdivisions[1]),
                      step[2]-1: 1.0-step[2]: complex(0, subdivisions[2])]
    means = np.reshape(means, [3, -1]).T
    covariances = variance*np.ones_like(means)
    weights = (1.0/n_gaussians)*np.ones(n_gaussians)
    gmm = GaussianMixture(n_components=n_gaussians, covariance_type='diag')
    gmm.weights_ = weights
    gmm.covariances_ = covariances
    gmm.means_ = means
    from sklearn.mixture.gaussian_mixture import _compute_precision_cholesky
    gmm.precisions_cholesky_ = _compute_precision_cholesky(covariances, 'diag')
    return gmm

def get_2d_grid_gmm(subdivisions=[5, 5], variance=0.04):
    """
    Compute the weight, mean and covariance of a 2D gmm placed on a 2D grid

    :param subdivisions: 2 element list of number of subdivisions of the 2D space in each axes to form the grid
    :param variance: scalar for spherical gmm.p
    :return gmm: gmm: instance of sklearn GaussianMixture (GMM) object Gauassian mixture model
    """
    # n_gaussians = reduce(lambda x, y: x*y,subdivisions)
    n_gaussians = np.prod(np.array(subdivisions))
    step = [1.0/(subdivisions[0]),  1.0/(subdivisions[1])]

    means = np.mgrid[step[0]-1: 1.0-step[0]: complex(0, subdivisions[0]),
            step[1]-1: 1.0-step[1]: complex(0, subdivisions[1])]
    means = np.reshape(means, [2,-1]).T
    covariances = variance*np.ones_like(means)
    weights = (1.0/n_gaussians)*np.ones(n_gaussians)
    gmm = GaussianMixture(n_components=n_gaussians, covariance_type='diag')
    gmm.weights_ = weights
    gmm.covariances_ = covariances
    gmm.means_ = means
    from sklearn.mixture.gaussian_mixture import _compute_precision_cholesky
    gmm.precisions_cholesky_ = _compute_precision_cholesky(covariances, 'diag')
    return gmm

def get_fisher_vectors(points,gmm, normalization=True):
    """
    Compute the fisher vector representation of a point cloud or a batch of point clouds

    :param points: n_points x 3 / B x n_points x 3
    :param gmm: sklearn MixtureModel class containing the gmm.p parameters.p
    :return: fisher vector representation for a single point cloud or a batch of point clouds
    """

    if len(points.shape) == 2:
        # single point cloud
        fv = fisher_vector(points, gmm, normalization=normalization)
    else:
        # Batch of  point clouds
        fv = []
        n_models = points.shape[0]
        for i in range(n_models):
            fv.append(fisher_vector(points[i], gmm, normalization=True))
        fv = np.array(fv)
    return fv

def fisher_vector(xx, gmm, normalization=True):
    """
    Computes the Fisher vector on a set of descriptors.
    code from : https://gist.github.cnsom/danoneata/9927923
    Parameters
    ----------
    xx: array_like, shape (N, D) or (D, )
        The set of descriptors

    gmm: instance of sklearn mixture.GMM object
        Gauassian mixture model of the descriptors.

    Returns
    -------
    fv: array_like, shape (K + 128 * D * K, )
        Fisher vector (derivatives with respect to the mixing weights, means
        and variances) of the given descriptors.

    Reference
    ---------
    Sanchez, J., Perronnin, F., Mensink, T., & Verbeek, J. (2013).
    Image classification with the fisher vector: Theory and practice. International journal of computer vision, 105(64), 222-245.
    https://hal.inria.fr/hal-00830491/file/journal.pdf

    """
    xx = np.atleast_2d(xx)
    n_points = xx.shape[0]
    D = gmm.means_.shape[1]
    tiled_weights = np.tile(np.expand_dims(gmm.weights_, axis=-1), [1, D])

    #start = time.time()
    # Compute posterior probabilities.
    Q = gmm.predict_proba(xx)  # NxK
    #mid = time.time()
    #print("Computing the probabilities took ", str(mid-start))
    #Compute Derivatives

    # Compute the sufficient statistics of descriptors.
    s0 = np.sum(Q, 0)[:, np.newaxis] / n_points
    s1 = np.dot(Q.T, xx) / n_points
    s2 = np.dot(Q.T, xx ** 2) / n_points

    d_pi = (s0.squeeze() - n_points * gmm.weights_) / np.sqrt(gmm.weights_)
    d_mu = (s1 - gmm.means_ * s0 ) / np.sqrt(tiled_weights*gmm.covariances_)
    d_sigma = (
        + s2
        - 2 * s1 * gmm.means_
        + s0 * gmm.means_ ** 2
        - s0 * gmm.covariances_
        ) / (np.sqrt(2*tiled_weights)*gmm.covariances_)

    #Power normaliation
    alpha = 0.5
    d_pi = np.sign(d_pi) * np.power(np.absolute(d_pi),alpha)
    d_mu = np.sign(d_mu) * np.power(np.absolute(d_mu), alpha)
    d_sigma = np.sign(d_sigma) * np.power(np.absolute(d_sigma), alpha)

    if normalization == True:
        d_pi = normalize(d_pi[:,np.newaxis], axis=0).ravel()
        d_mu = normalize(d_mu, axis=0)
        d_sigma = normalize(d_sigma, axis=0)
    # Merge derivatives into a vector.

    #print("comnputing the derivatives took ", str(time.time()-mid))

    return np.hstack((d_pi, d_mu.flatten(), d_sigma.flatten()))

def fisher_vector_per_point( xx, gmm):
    """
    see notes for above function - performs operations per point

    :param xx: array_like, shape (N, D) or (D, )- The set of descriptors
    :param gmm: instance of sklearn mixture.GMM object - Gauassian mixture model of the descriptors.
    :return: fv_per_point : fisher vector per point (derivative by w, derivative by mu, derivative by sigma)
    """
    xx = np.atleast_2d(xx)
    n_points = xx.shape[0]
    n_gaussians = gmm.means_.shape[0]
    D = gmm.means_.shape[1]

    sig2 = np.array([gmm.covariances_.T[0, :], gmm.covariances_.T[1, :], gmm.covariances_.T[2,:]]).T
    sig2_tiled = np.tile(np.expand_dims(sig2, axis=0), [n_points, 1, 1])

    # Compute derivativees per point and then sum.
    Q = gmm.predict_proba(xx)  # NxK
    tiled_weights = np.tile(np.expand_dims(gmm.weights_, axis=-1), [1, D])
    sqrt_w = np.sqrt(tiled_weights)

    d_pi = (Q - np.tile(np.expand_dims(gmm.weights_, 0), [n_points, 1])) / np.sqrt(np.tile(np.expand_dims(gmm.weights_, 0), [n_points, 1]))
    x_mu = np.tile( np.expand_dims(xx, axis=2), [1, 1, n_gaussians]) - np.tile(np.expand_dims(gmm.means_.T, axis=0), [n_points, 1, 1])
    x_mu = np.swapaxes(x_mu, 1, 2)
    d_mu = (np.tile(np.expand_dims(Q, -1), D) * x_mu) / (np.sqrt(sig2_tiled) * sqrt_w)

    d_sigma =   np.tile(np.expand_dims(Q, -1), 3)*((np.power(x_mu,2)/sig2_tiled)-1)/(np.sqrt(2)*sqrt_w)

    fv_per_point = (d_pi, d_mu, d_sigma)
    return fv_per_point

def l2_normalize(v, dim=1):
    """
    Normalize a vector along a dimension

    :param v: a vector or matrix to normalize
    :param dim: the dimension along which to normalize
    :return: normalized v along dim
    """
    norm = np.linalg.norm(v, axis=dim)
    if norm.all() == 0:
       return v
    return v / norm

def get_3DmFV(points, w, mu, sigma, normalize=True):
    """
       Compute the 3D modified fisher vectors given the gmm model parameters (w,mu,sigma) and a set of points
       For faster performance (large batches) use the tensorflow version

       :param points: B X N x 3 tensor of XYZ points
       :param w: B X n_gaussians tensor of gaussian weights
       :param mu: B X n_gaussians X 3 tensor of gaussian cetnters
       :param sigma: B X n_gaussians X 3 tensor of stddev of diagonal covariance
       :return: fv: B X 20*n_gaussians tensor of the fisher vector
       """
    n_batches = points.shape[0]
    n_points = points.shape[1]
    n_gaussians = mu.shape[0]
    D = mu.shape[1]

    # Expand dimension for batch compatibility
    batch_sig = np.tile(np.expand_dims(sigma, 0), [n_points, 1, 1])  # n_points X n_gaussians X D
    batch_sig = np.tile(np.expand_dims(batch_sig, 0), [n_batches, 1, 1, 1])  # n_batches X n_points X n_gaussians X D
    batch_mu = np.tile(np.expand_dims(mu, 0), [n_points, 1, 1])  # n_points X n_gaussians X D
    batch_mu = np.tile(np.expand_dims(batch_mu, 0), [n_batches, 1, 1, 1])  # n_batches X n_points X n_gaussians X D
    batch_w = np.tile(np.expand_dims(np.expand_dims(w, 0), 0), [n_batches, n_points,
                                                                1])  # n_batches X n_points X n_guassians X D  - should check what happens when weights change
    batch_points = np.tile(np.expand_dims(points, -2), [1, 1, n_gaussians,
                                                        1])  # n_batchesXn_pointsXn_gaussians_D  # Generating the number of points for each gaussian for separate computation

    # Compute derivatives
    w_per_batch_per_d = np.tile(np.expand_dims(np.expand_dims(w, 0), -1),
                                [n_batches, 1, 3*D])  # n_batches X n_gaussians X 3*D (D for min and D for max)

    # Define multivariate noraml distributions
    # Compute probability per point
    p_per_point = (1.0 / (np.power(2.0 * np.pi, D / 2.0) * np.power(batch_sig[:, :, :, 0], D))) * np.exp(
        -0.5 * np.sum(np.square((batch_points - batch_mu) / batch_sig), axis=3))

    w_p = p_per_point
    Q = w_p  # enforcing the assumption that the sum is 1
    Q_per_d = np.tile(np.expand_dims(Q, -1), [1, 1, 1, D])

    d_pi_all = np.expand_dims((Q - batch_w) / (np.sqrt(batch_w)), -1)
    d_pi = np.concatenate([np.max(d_pi_all, axis=1), np.sum(d_pi_all, axis=1)], axis=2)

    d_mu_all = Q_per_d * (batch_points - batch_mu) / batch_sig
    d_mu = (1 / (np.sqrt(w_per_batch_per_d))) * np.concatenate([np.max(d_mu_all, axis=1), np.min(d_mu_all, axis=1), np.sum(d_mu_all, axis=1)], axis=2)

    d_sig_all = Q_per_d * (np.square((batch_points - batch_mu) / batch_sig) - 1)
    d_sigma = (1 / (np.sqrt(2 * w_per_batch_per_d))) * np.concatenate([np.max(d_sig_all, axis=1), np.min(d_sig_all, axis=1), np.sum(d_sig_all, axis=1)], axis=2)

    # number of points  normaliation
    d_pi = d_pi / n_points
    d_mu = d_mu / n_points
    d_sigma =d_sigma / n_points

    if normalize:
        # Power normaliation
        alpha = 0.5
        d_pi = np.sign(d_pi) * np.power(np.abs(d_pi), alpha)
        d_mu = np.sign(d_mu) * np.power(np.abs(d_mu), alpha)
        d_sigma = np.sign(d_sigma) * np.power(np.abs(d_sigma), alpha)

        # L2 normaliation
        d_pi = np.array([l2_normalize(d_pi[i, :, :], dim=0) for i in range(n_batches)])
        d_mu = np.array([l2_normalize(d_mu[i, :, :], dim=0) for i in range(n_batches)])
        d_sigma = np.array([l2_normalize(d_sigma[i, :, :], dim=0) for i in range(n_batches)])

    fv = np.concatenate([d_pi, d_mu, d_sigma], axis=2)
    fv = np.transpose(fv, axes=[0, 2, 1])

    return fv

if __name__ == "__main__":

    model_idx = 0
    num_points = 1024
    gmm = get_3d_grid_gmm(subdivisions=[5, 5, 5], variance=0.04)
    points = provider.load_single_model(model_idx=model_idx, train_file_idxs=0, num_points=num_points)
    points = np.tile(np.expand_dims(points, 0), [128, 1, 1])

    fv_gpu = get_fisher_vectors(points, gmm, normalization=True)


text:
import torch
import torch.nn as nn
from torch.distributions import Normal
import numpy as np
import torch.nn.functional as F

def hidden_init(layer):
    fan_in = layer.weight.data.size()[0]
    lim = 1. / np.sqrt(fan_in)
    return (-lim, lim)

class Actor(nn.Module):
    """Actor (Policy) Model."""

    def __init__(self, state_size, action_size, seed, device, hidden_size=32, init_w=3e-3, log_std_min=-20, log_std_max=2):
        """Initialize parameters and build model.
        Params
        ======
            state_size (int): Dimension of each state
            action_size (int): Dimension of each action
            seed (int): Random seed
            fc1_units (int): Number of nodes in first hidden layer
            fc2_units (int): Number of nodes in second hidden layer
        """
        super(Actor, self).__init__()
        self.seed = torch.manual_seed(seed)
        self.device = device
        self.log_std_min = log_std_min
        self.log_std_max = log_std_max

        self.fc1 = nn.Linear(state_size, hidden_size)
        self.fc2 = nn.Linear(hidden_size, hidden_size)

        self.mu = nn.Linear(hidden_size, action_size)
        self.log_std_linear = nn.Linear(hidden_size, action_size)

    def reset_parameters(self):
        self.fc1.weight.data.uniform_(*hidden_init(self.fc1))
        self.fc2.weight.data.uniform_(*hidden_init(self.fc2))
        self.mu.weight.data.uniform_(-init_w, init_w)
        self.log_std_linear.weight.data.uniform_(-init_w, init_w)

    def forward(self, state):

        x = F.relu(self.fc1(state))
        x = F.relu(self.fc2(x))
        mu = self.mu(x)

        log_std = self.log_std_linear(x)
        log_std = torch.clamp(log_std, self.log_std_min, self.log_std_max)
        return mu, log_std

    def evaluate(self, state, epsilon=1e-6):
        mu, log_std = self.forward(state)
        std = log_std.exp()
        dist = Normal(mu, std)
        e = dist.rsample().to(self.device)
        action = torch.tanh(e)
        log_prob = (dist.log_prob(e) - torch.log(1 - action.pow(2) + epsilon)).sum(1, keepdim=True)

        return action, log_prob

    def get_action(self, state):
        """
        returns the action based on a squashed gaussian policy. That means the samples are obtained according to:
        a(s,e)= tanh(mu(s)+sigma(s)+e)
        """
        mu, log_std = self.forward(state)
        std = log_std.exp()
        dist = Normal(mu, std)
        e = dist.rsample().to(self.device)
        action = torch.tanh(e)
        return action.detach().cpu()

    def get_det_action(self, state):
        mu, log_std = self.forward(state)
        return torch.tanh(mu).detach().cpu()

class Critic(nn.Module):
    """Critic (Value) Model."""

    def __init__(self, state_size, action_size, seed, device, hidden_size=32):
        """Initialize parameters and build model.
        Params
        ======
            state_size (int): Dimension of each state
            action_size (int): Dimension of each action
            seed (int): Random seed
            hidden_size (int): Number of nodes in the network layers

        """
        super(Critic, self).__init__()
        self.seed = torch.manual_seed(seed)
        self.device = device
        self.fc1 = nn.Linear(state_size+action_size, hidden_size)
        self.fc2 = nn.Linear(hidden_size, hidden_size)
        self.fc3 = nn.Linear(hidden_size, 1)
        self.reset_parameters()

    def reset_parameters(self):
        self.fc1.weight.data.uniform_(*hidden_init(self.fc1))
        self.fc2.weight.data.uniform_(*hidden_init(self.fc2))
        self.fc3.weight.data.uniform_(-3e-3, 3e-3)

    def forward(self, state, action):
        """Build a critic (value) network that maps (state, action) pairs -> Q-values."""
        x = torch.cat((state, action), dim=1)
        x = F.relu(self.fc1(x))
        x = F.relu(self.fc2(x))
        return self.fc3(x)

class DeepActor(nn.Module):
    """Actor (Policy) Model."""

    def __init__(self, state_size, action_size, seed, device, hidden_size=32, init_w=3e-3, log_std_min=-20, log_std_max=2):
        """Initialize parameters and build model.
        Params
        ======
            state_size (int): Dimension of each state
            action_size (int): Dimension of each action
            seed (int): Random seed
            fc1_units (int): Number of nodes in first hidden layer
            fc2_units (int): Number of nodes in second hidden layer
        """
        super(DeepActor, self).__init__()
        self.seed = torch.manual_seed(seed)
        self.device = device
        self.log_std_min = log_std_min
        self.log_std_max = log_std_max

        in_dim = hidden_size+state_size

        self.fc1 = nn.Linear(state_size, hidden_size)
        self.fc2 = nn.Linear(in_dim, hidden_size)
        self.fc3 = nn.Linear(in_dim, hidden_size)
        self.fc4 = nn.Linear(in_dim, hidden_size)

        self.mu = nn.Linear(hidden_size, action_size)
        self.log_std_linear = nn.Linear(hidden_size, action_size)
        #self.reset_parameters() # check if this improves training

    def reset_parameters(self, init_w=3e-3):
        self.fc1.weight.data.uniform_(*hidden_init(self.fc1))
        self.fc2.weight.data.uniform_(*hidden_init(self.fc2))
        self.fc3.weight.data.uniform_(*hidden_init(self.fc3))
        self.fc4.weight.data.uniform_(*hidden_init(self.fc4))
        self.mu.weight.data.uniform_(-init_w, init_w)
        self.log_std_linear.weight.data.uniform_(-init_w, init_w)

    def forward(self, state: torch.tensor)-> (float, float):

        x = F.relu(self.fc1(state))
        x = torch.cat([x, state], dim=1)
        x = F.relu(self.fc2(x))
        x = torch.cat([x, state], dim=1)
        x = F.relu(self.fc3(x))
        x = torch.cat([x, state], dim=1)
        x = F.relu(self.fc4(x))

        mu = self.mu(x)

        log_std = self.log_std_linear(x)
        log_std = torch.clamp(log_std, self.log_std_min, self.log_std_max)
        return mu, log_std

    def evaluate(self, state, epsilon=1e-6):
        mu, log_std = self.forward(state)
        std = log_std.exp()
        dist = Normal(mu, std)
        e = dist.rsample().to(self.device)
        action = torch.tanh(e)
        log_prob = (dist.log_prob(e) - torch.log(1 - action.pow(2) + epsilon)).sum(1, keepdim=True)

        return action, log_prob

    def get_action(self, state):
        """
        returns the action based on a squashed gaussian policy. That means the samples are obtained according to:
        a(s,e)= tanh(mu(s)+sigma(s)+e)
        """
        #state = torch.FloatTensor(state).to(device) #.unsqzeeze(0)
        mu, log_std = self.forward(state)
        std = log_std.exp()
        dist = Normal(mu, std)
        e = dist.rsample().to(self.device)
        action = torch.tanh(e)
        return action.detach().cpu()

    def get_det_action(self, state):
        mu, log_std = self.forward(state)
        return torch.tanh(mu).detach().cpu()

class DeepCritic(nn.Module):
    """Critic (Value) Model."""

    def __init__(self, state_size, action_size, seed, device, hidden_size=32):
        """Initialize parameters and build model.
        Params
        ======
            state_size (int): Dimension of each state
            action_size (int): Dimension of each action
            seed (int): Random seed
            hidden_size (int): Number of nodes in the network layers

        """
        super(DeepCritic, self).__init__()
        self.seed = torch.manual_seed(seed)
        self.device = device
        in_dim = hidden_size+action_size+state_size
        self.fc1 = nn.Linear(state_size+action_size, hidden_size)
        self.fc2 = nn.Linear(in_dim, hidden_size)
        self.fc3 = nn.Linear(in_dim, hidden_size)
        self.fc4 = nn.Linear(in_dim, hidden_size)

        self.fc5 = nn.Linear(hidden_size, 1)
        self.reset_parameters()

    def reset_parameters(self):
        self.fc1.weight.data.uniform_(*hidden_init(self.fc1))
        self.fc2.weight.data.uniform_(*hidden_init(self.fc2))
        self.fc3.weight.data.uniform_(*hidden_init(self.fc2))
        self.fc4.weight.data.uniform_(*hidden_init(self.fc2))
        self.fc5.weight.data.uniform_(-3e-3, 3e-3)

    def forward(self, state, action):
        """Build a critic (value) network that maps (state, action) pairs -> Q-values."""
        xu = torch.cat((state, action), dim=1)
        x = F.relu(self.fc1(xu))
        x = torch.cat([x, xu], dim=1)
        x = F.relu(self.fc2(x))
        x = torch.cat([x, xu], dim=1)
        x = F.relu(self.fc3(x))
        x = torch.cat([x, xu], dim=1)
        x = F.relu(self.fc4(x))

        return self.fc5(x)

class IQN(nn.Module):
    def __init__(self, state_size, action_size, layer_size, seed, N, device="cuda:0"):
        super(IQN, self).__init__()
        self.seed = torch.manual_seed(seed)
        self.input_shape = state_size
        self.action_size = action_size
        self.N = N
        self.n_cos = 64
        self.layer_size = layer_size
        self.pis = torch.FloatTensor([np.pi*i for i in range(1,self.n_cos+1)]).view(1,1,self.n_cos).to(device) # Starting from 0 as in the paper
        self.device = device

        # Network Architecture

        self.head = nn.Linear(self.action_size+self.input_shape, layer_size)
        self.cos_embedding = nn.Linear(self.n_cos, layer_size)
        self.ff_1 = nn.Linear(layer_size, layer_size)
        self.ff_2 = nn.Linear(layer_size, 1)
        #weight_init([self.head_1, self.ff_1])

    def calc_input_layer(self):
        x = torch.zeros(self.input_shape).unsqueeze(0)
        x = self.head(x)
        return x.flatten().shape[0]

    def calc_cos(self, batch_size, n_tau=32):
        """
        Calculating the cosinus values depending on the number of tau samples
        """
        taus = torch.rand(batch_size, n_tau).unsqueeze(-1).to(self.device) #(batch_size, n_tau, 1)  .to(self.device)
        cos = torch.cos(taus*self.pis)

        assert cos.shape == (batch_size,n_tau,self.n_cos), "cos shape is incorrect"
        return cos, taus

    def forward(self, input, action, num_tau=32):
        """
        Quantile Calculation depending on the number of tau

        Return:
        quantiles [ shape of (batch_size, num_tau, action_size)]
        taus [shape of ((batch_size, num_tau, 1))]

        """
        batch_size = input.shape[0]

        x = torch.cat((input, action), dim=1)
        x = torch.relu(self.head(x  ))

        cos, taus = self.calc_cos(batch_size, num_tau) # cos shape (batch, num_tau, layer_size)
        cos = cos.view(batch_size*num_tau, self.n_cos)
        cos_x = torch.relu(self.cos_embedding(cos)).view(batch_size, num_tau, self.layer_size) # (batch, n_tau, layer)

        # x has shape (batch, layer_size) for multiplication –> reshape to (batch, 1, layer)
        x = (x.unsqueeze(1)*cos_x).view(batch_size*num_tau, self.layer_size)  #batch_size*num_tau, self.cos_layer_out
        # Following reshape and transpose is done to bring the action in the same shape as batch*tau:
        # first 32 entries are tau for each action -> thats why each action one needs to be repeated 32 times
        # x = [[tau1   action = [[a1
        #       tau1              a1
        #        ..               ..
        #       tau2              a2
        #       tau2              a2
        #       ..]]              ..]]
        #action = action.repeat(num_tau,1).reshape(num_tau,batch_size*self.action_size).transpose(0,1).reshape(batch_size*num_tau,self.action_size)
        #x = torch.cat((x,action),dim=1)
        x = torch.relu(self.ff_1(x))

        out = self.ff_2(x)

        return out.view(batch_size, num_tau, 1), taus

    def get_qvalues(self, inputs, action):
        quantiles, _ = self.forward(inputs, action, self.N)
        actions = quantiles.mean(dim=1)
        return actions

class DeepIQN(nn.Module):
    def __init__(self, state_size, action_size, layer_size, seed, N, device="cuda:0"):
        super(DeepIQN, self).__init__()
        self.seed = torch.manual_seed(seed)
        self.input_shape = state_size
        self.action_size = action_size
        self.input_dim = action_size+state_size+layer_size
        self.N = N
        self.n_cos = 64
        self.layer_size = layer_size
        self.pis = torch.FloatTensor([np.pi*i for i in range(1,self.n_cos+1)]).view(1,1,self.n_cos).to(device) # Starting from 0 as in the paper
        self.device = device

        # Network Architecture

        self.head = nn.Linear(self.action_size+self.input_shape, layer_size)
        self.ff_1 = nn.Linear(self.input_dim, layer_size)
        self.ff_2 = nn.Linear(self.input_dim, layer_size)
        self.cos_embedding = nn.Linear(self.n_cos, layer_size)
        self.ff_3 = nn.Linear(self.input_dim, layer_size)
        self.ff_4 = nn.Linear(self.layer_size, 1)
        #weight_init([self.head_1, self.ff_1])

    def calc_input_layer(self):
        x = torch.zeros(self.input_shape).unsqueeze(0)
        x = self.head(x)
        return x.flatten().shape[0]

    def calc_cos(self, batch_size, n_tau=32):
        """
        Calculating the cosinus values depending on the number of tau samples
        """
        taus = torch.rand(batch_size, n_tau).unsqueeze(-1).to(self.device) #(batch_size, n_tau, 1)  .to(self.device)
        cos = torch.cos(taus*self.pis)

        assert cos.shape == (batch_size,n_tau,self.n_cos), "cos shape is incorrect"
        return cos, taus

    def forward(self, input, action, num_tau=32):
        """
        Quantile Calculation depending on the number of tau

        Return:
        quantiles [ shape of (batch_size, num_tau, action_size)]
        taus [shape of ((batch_size, num_tau, 1))]

        """
        batch_size = input.shape[0]
        xs = torch.cat((input, action), dim=1)
        x = torch.relu(self.head(xs))
        x = torch.cat((x, xs), dim=1)
        x = torch.relu(self.ff_1(x))
        x = torch.cat((x, xs), dim=1)
        x = torch.relu(self.ff_2(x))

        cos, taus = self.calc_cos(batch_size, num_tau) # cos shape (batch, num_tau, layer_size)
        cos = cos.view(batch_size*num_tau, self.n_cos)
        cos_x = torch.relu(self.cos_embedding(cos)).view(batch_size, num_tau, self.layer_size) # (batch, n_tau, layer)

        # x has shape (batch, layer_size) for multiplication –> reshape to (batch, 1, layer)
        x = (x.unsqueeze(1)*cos_x).view(batch_size*num_tau, self.layer_size)  #batch_size*num_tau, self.cos_layer_out
        # Following reshape and transpose is done to bring the action in the same shape as batch*tau:
        # first 32 entries are tau for each action -> thats why each action one needs to be repeated 32 times
        # x = [[tau1   action = [[a1
        #       tau1              a1
        #        ..               ..
        #       tau2              a2
        #       tau2              a2
        #       ..]]              ..]]
        action = action.repeat(num_tau,1).reshape(num_tau,batch_size*self.action_size).transpose(0,1).reshape(batch_size*num_tau,self.action_size)
        state = input.repeat(num_tau,1).reshape(num_tau,batch_size*self.input_shape).transpose(0,1).reshape(batch_size*num_tau,self.input_shape)

        x = torch.cat((x,action,state),dim=1)
        x = torch.relu(self.ff_3(x))

        out = self.ff_4(x)

        return out.view(batch_size, num_tau, 1), taus

    def get_qvalues(self, inputs, action):
        quantiles, _ = self.forward(inputs, action, self.N)
        actions = quantiles.mean(dim=1)
        return actions


text:
"""
The river class code.
"""

import pandas as pd
import urllib
import xml.etree.ElementTree as ET
from datetime import datetime
#import numpy as np

class River:

    def __init__(self, riv_name, rs):
        # initialize all expected fields
        self.att_list = ['name', 'usgs_code',
            'nws_code', 'ec_code', 'scale_factor',
            'got_data', 'memo', 'flow_units']
        for att in self.att_list:
            setattr(self, att, '')
        self.name = riv_name
        self.rs = rs
        # default values
        self.Q = []
        self.T = []
        self.qt = pd.Series(self.Q, index=self.T)
        try:
            self.scale_factor = self.rs['ratio']
            self.usgs_code = self.rs['usgs']
            self.nws_code = self.rs['nws']
            self.ec_code = self.rs['ec']
        except KeyError:
            # needed for analytical rivers
            pass
        self.got_data = False
        self.memo = 'no message'

    def print_info(self):
        print(50*'-')
        for att in self.att_list:
            print(att + ' = ' + str(getattr(self,att)))
        print(50*'-')

    def get_usgs_data(self, days):
        # this deals with rivers that need gauge combinations
        if self.name == 'skokomish':
            print('+++ combining to form Skokomish River +++')
            self.usgs_code = 12060500
            self.scale_factor = 1.4417
            self.get_usgs_data_sub(days)
            qt1 = self.qt.copy()
            self.usgs_code = 12059500
            self.scale_factor = 1.0
            self.get_usgs_data_sub(days)
            qt2 = self.qt.copy()
            self.qt = qt1 + qt2
        elif self.name == 'hamma':
            print('+++ combining to form Hamma Hamma River +++')
            self.usgs_code = 12060500
            self.scale_factor = 1.4417
            self.get_usgs_data_sub(days)
            qt1 = self.qt.copy()
            self.usgs_code = 12059500
            self.scale_factor = 1.0
            self.get_usgs_data_sub(days)
            qt2 = self.qt.copy()
            self.qt = 0.4125 * (qt1 + qt2)
        else:
            self.get_usgs_data_sub(days)

    def get_usgs_data_sub(self, days):
        # This gets USGS data for a past time period specfied by
        # the tuple of datetimes "days".  If "days" is empty
        # then we get the most recent 6 days.

        # Set the time period to get.
        if len(days) == 2:
            time_str = ('&startDT=' + days[0].strftime('%Y-%m-%d')
                +'&endDT=' + days[1].strftime('%Y-%m-%d'))
        else:
            # This gets the most recent 6 days (daily intervals?)
            time_str = '&period=P6D'
        gage = int(self.usgs_code)
        # Form the url.
        url_str = ('http://waterservices.usgs.gov/nwis/dv/'
            + '?format=waterml,1.1&sites=' + str(gage)
            + time_str + '&parameterCd=00060')
        try:
            # Get the XML.
            file = urllib.request.urlopen(url_str, timeout=10)
            tree = ET.parse(file)
            root = tree.getroot()
            # Extract the data from the XML.
            flag = True
            # aa = '{http://www.cuahsi.org/waterML/1.1/}'
            rt = root.tag
            aa = rt[rt.find('{'): rt.find('}') + 1]
            for e0 in root.findall(".//"):
                if e0.tag == aa+'value':
                    self.Q.append(float(e0.text))
                    self.T.append(pd.to_datetime(e0.get('dateTime')))
                if e0.tag == aa+'unitCode' and flag:
                    self.flow_units = e0.text
                    flag = False
            self.qt = pd.Series(self.Q, index=self.T)
            self.fix_units()
            self.qt = float(self.scale_factor) * self.qt
            # Note: this resampling of daily data just moves the timestamp to noon
            # of the day it started at.  Data is unchanged.
            self.qt = self.qt.resample('D', label='right', loffset='-12h').mean()
            self.got_data = True
            self.memo = 'success'
        except:
            self.memo = 'problem downloading or parsing data from XML'
        self.memo = (self.memo + ' USGS')

    def get_nws_data(self):
        # This gets NWS forecast data.
        url_str = ('http://www.nwrfc.noaa.gov/xml/xml.cgi?id=' +
                   self.nws_code +
                   '&pe=HG&dtype=b&numdays=10')
        try:
            # get the XML
            file = urllib.request.urlopen(url_str, timeout=10)
            tree = ET.parse(file)
            root = tree.getroot()
            flag = True
            # NOTE: you find this tag by looking at any instance of e0.tag
            #aa = '{http://www.nwrfc.noaa.gov/xml/schemas/2004/03/hydromet_data}'
            rt = root.tag
            aa = rt[rt.find('{'): rt.find('}') + 1]
            for e0 in root.findall(".//"):
                if e0.tag == aa+'observedData' or e0.tag == aa+'forecastData':
                    for e in e0:
                        if e.tag == aa+'observedValue' or e.tag == aa+'forecastValue':
                            for ee in e:
                                if ee.tag == aa+'discharge':
                                    self.Q.append(float(ee.text))
                                    if flag:
                                        self.flow_units = ee.get('units')
                                        flag = False
                                if ee.tag == aa+'dataDateTime':
                                    self.T.append(pd.to_datetime(ee.text))
            self.qt = pd.Series(self.Q, index=self.T)
            self.fix_units()
            self.qt = float(self.scale_factor) * self.qt
            self.qt = self.qt.resample('D', label='right', loffset='-12h').mean()
            self.got_data = True
            self.memo = 'success'
        except:
            self.memo = 'problem downloading or parsing data from XML'
        self.memo = (self.memo + ' NWS')

    def get_ec_data(self, days, timeout=10):
        # gets Environment Canada data, using code cribbed from:
        #https://bitbucket.org/douglatornell/ecget/src/
        # NOTE: this will get data up through today, but can only go back
        # 18 months into the past.
        # To get longer records use get_ec_data_historical.
        import requests
        import bs4
        try:
            PARAM_IDS = {'discharge': 47,}

            params = {
                'mode': 'Table',
                'type': 'realTime',
                'prm1': PARAM_IDS['discharge'],
                'prm2': -1,
                'stn': self.ec_code,
                'startDate': days[0].strftime('%Y-%m-%d'),
                'endDate': days[1].strftime('%Y-%m-%d'),
            }
            DATA_URL = 'http://wateroffice.ec.gc.ca/report/real_time_e.html'
            DISCLAIMER_COOKIE = {'disclaimer': 'agree'}
            response = requests.get(DATA_URL, params=params,
                                    cookies=DISCLAIMER_COOKIE)
            soup = bs4.BeautifulSoup(response.content, 'lxml')
            table = soup.find('table')
            table_body = table.find('tbody')
            rows = table_body.find_all('tr')
            data = []
            for row in rows:
                cols = row.find_all('td')
                cols = [ele.text.strip() for ele in cols]
                data.append([ele for ele in cols if ele]) # Get rid of empty values
            d_dict = dict()
            for item in data:
                d_dict[pd.to_datetime(item[0])] = float(item[1].replace(',',''))
            self.qt = pd.Series(d_dict)
            self.flow_units = '$m^{3}s^{-1}$'
            #self.fix_units() # not needed
            self.qt = float(self.scale_factor) * self.qt
            self.qt = self.qt.resample('D', label='right', loffset='-12h').mean()
            # NEW 2019.03.20 to deal with the problem that when you request date from
            # before 18 months ago it gives the most recent data instead.
            dt0_actual = self.qt.index[0]
            dt0_requested = days[0]
            import numpy as np
            if np.abs((dt0_actual - dt0_requested).days) >= 1:
                memo = 'That date range was not available'
                qt = ''
            else:
                self.got_data = True
                self.memo = 'success'
        except:
            self.memo = 'problem parsing data from soup'
        self.memo = (self.memo + ' EC')

    def get_ec_data_historical(self, year):
        # gets Environment Canada data, using code cribbed from:
        #https://bitbucket.org/douglatornell/ecget/src/
        # NOTE: this will get data up through the end of 2016.
        import requests
        import bs4
        try:
            params = {
                'mode': 'Table',
                'type': 'h2oArc',
                'stn': self.ec_code,
                'dataType': 'Daily',
                'parameterType': 'Flow',
                'year': str(year),
                'y1Max': '1',
                'y1Min': '1',
            }
            DATA_URL = 'http://wateroffice.ec.gc.ca/report/historical_e.html'
            DISCLAIMER_COOKIE = {'disclaimer': 'agree'}
            response = requests.get(DATA_URL, params=params,
                                    cookies=DISCLAIMER_COOKIE)
            soup = bs4.BeautifulSoup(response.content, 'lxml')
            if (str(year) + ' Daily Discharge') in soup.text:
                # we do the test above because the website will return the
                # table for the most recent year if the requested year is
                # missing
                table = soup.find('table')
                table_body = table.find('tbody')
                rows = table_body.find_all('tr')
                d_dict = dict()
                for row in rows:
                    # what day is it
                    for ele in row.find_all('th'):
                        day = int(ele.text.strip())
                    cols = row.find_all('td')
                    cols = [ele.text.strip() for ele in cols] # a list of strings
                    imo = 1
                    for item in cols:
                        if len(item) == 0:
                            pass
                        else:
                            this_data = float(item.split()[0].replace(',',''))
                            # the split call is to remove trailing letters
                            # that occasionally appear after the data
                            d_dict[datetime(year, imo, day, 12, 0, 0)] = this_data
                        imo += 1
                self.qt = pd.Series(d_dict)
                self.flow_units = '$m^{3}s^{-1}$'
                #self.fix_units() # not needed
                self.qt = float(self.scale_factor) * self.qt
                self.got_data = True
                self.memo = 'success'
            else:
                self.memo = 'That year was not available'
        except:
            self.memo = 'problem parsing data from soup'
        self.memo = (self.memo + ' EC')

    def fix_units(self):
        # fix units
        if self.flow_units == 'kcfs':
            self.qt = self.qt*28.3168466
            self.flow_units = '$m^{3}s^{-1}$'
        elif (self.flow_units == 'cubic feet per second'
            or self.flow_units == 'ft3/s'):
            self.qt = self.qt*0.0283168466
            self.flow_units = '$m^{3}s^{-1}$'
        else:
            pass

max_stars_repo_name: parkermac/LiveOcean

max_issues_repo_name: parkermac/LiveOcean

max_forks_repo_name: parkermac/LiveOcean


text:
from enum import Enum

import cadquery as cq
import numpy as np

from cq_cam.utils.utils import wire_to_ordered_edges

class Transition(Enum):
    NORMAL = 0
    TAB = 1
    IGNORE = 2

class Tabs:
    def __init__(self, *, width: float, height: float):
        self.width = width
        self.height = height

    def load_ordered_edges(self, ordered_edges):
        raise NotImplementedError()

    def edge_tab_transitions(self, edge_index):
        raise NotImplementedError()

class NoTabs(Tabs):
    def __init__(self):
        super().__init__(width=0, height=0)

    def load_ordered_edges(self, ordered_edges):
        pass

    def edge_tab_transitions(self, edge_index):
        return [
            (0, Transition.NORMAL),
            (1, Transition.NORMAL)
        ]

class EdgeTabs(Tabs):
    def __init__(self, *, spacing: float, width: float, height: float, only=None):
        super().__init__(width=width, height=height)
        self.edges_ds = None
        self.spacing = spacing
        self.only = only

    def load_ordered_edges(self, ordered_edges):
        self.edges_ds = []
        half = self.width / 2
        for edge in ordered_edges:
            length = edge.Length()
            edge_ds = []
            self.edges_ds.append(edge_ds)
            geom_type = edge.geomType()
            if self.only:
                if geom_type != self.only:
                    continue
            count = int(length // self.spacing)
            if count <= 0:
                continue

            half_d = 1 / length * half
            for d in np.linspace(0, 1, count + 2, endpoint=True)[1:-1]:
                edge_ds.append((d - half_d, d + half_d))

    def edge_tab_transitions(self, edge_index):
        edge_ds = self.edges_ds[edge_index]
        transitions = [(0, Transition.NORMAL)]
        for tab_start_d, tab_end_d in edge_ds:
            transitions.append((tab_start_d, Transition.TAB))
            transitions.append((tab_end_d, Transition.NORMAL))
        transitions.append((1, Transition.NORMAL))
        return transitions

class WireTabs(Tabs):
    def __init__(self, *, count: int, width: float, height: float):
        super().__init__(width=width, height=height)
        self.edge_lengths = None
        self.edge_ranges = None
        self.wire_tab_ds = None
        self.count = count
        self.width = width

    def load_wire(self, wire):
        length = wire.Length()
        half = self.width / 2
        half_d = 1 / length * half

        wire_tab_ds = []

        for d in np.linspace(0, 1, self.count, endpoint=False):
            range = (d - half_d, d + half_d)
            wire_tab_ds.append(range)

            # Handle edge cases of ranges crossing 0 or 1 by adding extra range(s)
            if range[0] < 0:
                wire_tab_ds.append((range[0] + 1, range[1] + 1))
            if range[1] > 1:
                wire_tab_ds.append((range[0] - 1, range[1] - 1))

        self.wire_tab_ds = wire_tab_ds

    def load_ordered_edges(self, ordered_edges):
        edge_lengths = [edge.Length() for edge in ordered_edges]
        wire_length = sum(edge_lengths)
        edge_ranges = []
        previous_range_end = 0
        for edge_length in edge_lengths:
            edge_fraction = edge_length / wire_length
            edge_range_end = previous_range_end + edge_fraction
            edge_ranges.append((previous_range_end, min(edge_range_end, 1)))
            previous_range_end = edge_range_end

        self.edge_lengths = edge_lengths
        self.edge_ranges = edge_ranges

    def edge_tab_transitions(self, edge_index):
        edge_range = self.edge_ranges[edge_index]
        edge_start_d, edge_end_d = edge_range
        transitions = []
        for tab_start_d, tab_end_d in self.wire_tab_ds:
            # Tab crosses start
            if tab_start_d <= edge_start_d < tab_end_d:
                transitions.append((0, Transition.TAB))

            # Tab starts inside edge
            if edge_start_d < tab_start_d <= edge_end_d:
                edge_d = min(max(self.wire_d_to_edge_d(tab_start_d, edge_range), 0), 1)
                transitions.append((edge_d, Transition.TAB))

            # Tab ends inside edge
            if edge_start_d < tab_end_d <= edge_end_d:
                edge_d = min(max(self.wire_d_to_edge_d(tab_end_d, edge_range), 0), 1)
                transitions.append((edge_d, Transition.NORMAL))

        # Add helper transitions for later processing
        if not transitions or transitions[0][0] != 0:
            transitions.insert(0, (0, Transition.NORMAL))
        if transitions[-1][0] != 1:
            transitions.append((1, Transition.IGNORE))
        return transitions

    @staticmethod
    def wire_edge_d_ranges(wire: cq.Wire):
        wire_length = wire.Length()
        edges = wire_to_ordered_edges(wire)
        edge_lengths = [edge.Length() for edge in edges]
        edge_ranges = []
        previous_range_end = 0
        for edge_length in edge_lengths:
            edge_fraction = edge_length / wire_length
            edge_range_end = previous_range_end + edge_fraction
            edge_ranges.append((previous_range_end, min(edge_range_end, 1)))
            previous_range_end = edge_range_end

        return edges, edge_ranges

    @staticmethod
    def wire_d_to_edge_d(wire_d, edge_range):
        return (wire_d - edge_range[0]) / (edge_range[1] - edge_range[0])

# TODO WireTabs EdgeTabs+ .,n
if 'show_object' in locals() or __name__ == '__main__':
    box = cq.Workplane().box(5, 5, 5)
    bottom = box.wires('<Z')
    tabs = Tabs(1, 1, 4)
    tabs.process(bottom.objects[0])
    show_object(box, 'box')

max_stars_repo_name: cyan-at/cq-cam

max_issues_repo_name: cyan-at/cq-cam

max_forks_repo_name: cyan-at/cq-cam


text:
'''
CBOR reader and writing functionality used by the `Artifact` class.

Exported definitions:
    PersistentArray (`numpy.memmap` subclass): A `memmap` backed by a CBOR file.
    PersistentList (`list` subclass): A `list` backed by a CBOR file.
    read_cbor_file (function): Read a CBOR file.
    write_object_as_cbor (function): Write an object to a CBOR file.
'''

from __future__ import annotations

import sys
from contextlib import contextmanager
from io import BufferedRandom
from itertools import chain
from os import SEEK_END, SEEK_SET
from pathlib import Path
from typing import Any, Iterable, Iterator, Sequence, Tuple, cast
from typing_extensions import Annotated

try:
    from fcntl import LOCK_EX, LOCK_SH, LOCK_UN, lockf
    locking_is_supported = True
except ImportError:
    locking_is_supported = False

import cbor2
import numpy as np

from ._namespaces import dictify, namespacify

__all__ = [
    'PersistentArray', 'PersistentList',
    'read_cbor_file', 'write_object_as_cbor']

#-- CBOR primitives ------------------------------------------------------------

MAJOR_TYPE_UINT = 0 << 5
MAJOR_TYPE_BYTE_STRING = 2 << 5
MAJOR_TYPE_ARRAY = 4 << 5
MAJOR_TYPE_TAG = 6 << 5

TAG_MULTIDIM_ARRAY = 40

INFO_NEXT_BYTE = 24
INFO_NEXT_2_BYTES = 25
INFO_NEXT_4_BYTES = 26
INFO_NEXT_8_BYTES = 27

dtypes_by_tag = {
    64: np.dtype('u1'),
    65: np.dtype('>u2'),
    66: np.dtype('>u4'),
    67: np.dtype('>u8'),
    68: np.dtype('u1'),
    69: np.dtype('<u2'),
    70: np.dtype('<u4'),
    71: np.dtype('<u8'),
    72: np.dtype('i1'),
    73: np.dtype('>i2'),
    74: np.dtype('>i4'),
    75: np.dtype('>i8'),
    77: np.dtype('<i2'),
    78: np.dtype('<i4'),
    79: np.dtype('<i8'),
    80: np.dtype('>f2'),
    81: np.dtype('>f4'),
    82: np.dtype('>f8'),
    84: np.dtype('<f2'),
    85: np.dtype('<f4'),
    86: np.dtype('<f8')}

tags_by_dtype = {
    dtype: tag
    for tag, dtype
    in dtypes_by_tag.items()}

#-- Persistent collections -----------------------------------------------------

class PersistentList(list):
    '''
    A `list` backed by a CBOR file.

    For performance, a `PersistentList` is invalidated when another object,
    including another `PersistentList`, writes to its backing file. An
    invalidated `PersistentList` is a potentially out-of-date read-only view
    into the file, and calling `append` or `extend` on it will corrupt the file.
    '''
    def __init__(self, file_: BufferedRandom, length: int) -> None:
        # Read the file into a buffer.
        file_.seek(0, SEEK_END)
        buf = bytearray(file_.tell())
        file_.seek(0, SEEK_SET)
        file_.readinto(buf)

        # Overwrite the header, in case it is currently being written to.
        header = list_header(length)
        buf[:len(header)] = header

        # Parse the buffer's contents as list items.
        super().__init__(namespacify(cbor2.loads(buf)))

        # Store the file pointer for `extend` calls.
        self._file = file_

    def __setitem__(self, index: object, value: object) -> None:
        raise TypeError('`PersistentList`s do not support item assignment')

    def __delitem__(self, index: object) -> None:
        raise TypeError('`PersistentList`s do not support item deletion')

    def extend(self, items: Iterable[object]) -> None:
        '''
        Extend list by appending elements from the iterable.
        '''
        # Coerce the collection of items to add into a sequence.
        items = items if isinstance(items, Sequence) else list(items)

        # Append the items, CBOR-encoded, to the backing file.
        data = b''.join(map(cbor2.dumps, dictify(items)))
        self._file.seek(0, SEEK_END)
        self._file.write(data)
        self._file.flush()

        # Update the header with the new list length.
        header = list_header(len(self) + len(items))
        with locking_header(self._file, LOCK_EX):
            self._file.seek(0, SEEK_SET)
            self._file.write(header)
            self._file.flush()

        # Add the items to `self`.
        super().extend(items)

    def append(self, item: object) -> None:
        '''
        Append object to the end of the list.
        '''
        self.extend([item])

class PersistentArray(np.memmap):
    '''
    A `numpy.memmap` backed by a CBOR file.

    The file must contain a row-major multidimensional array as defined in IETF
    RFC 8746. For performance, a `PersistentArray` is invalidated when another
    object, including another `PersistentArray`, writes to its backing file. An
    invalidated `PersistentArray` is a potentially out-of-date read-only view
    into the file, and calling `append` or `extend` on it will corrupt the file.

    Due to NumPy issue 4198 (https://github.com/numpy/numpy/issues/4198),
    `PersistenArray` extends `np.memmap` by proxy, meaning that it delegates
    attribute accesses and method calls to an internal `np.memmap` object
    instead of using Python's native subclassing mechanism.
    '''
    def extend(self, items: object) -> None:
        '''
        Extend the array by appending elements from `items`.
        '''

    def append(self, item: object) -> None:
        '''
        Append `item` to the array.
        '''

class PersistentArrayImpl:
    __name__ = 'PersistentArray'
    __qualname__ = 'PersistentArray'
    __doc__ = PersistentArray.__doc__

    def __init__(self,
                 file_: BufferedRandom,
                 shape: Tuple[int, ...],
                 dtype: np.dtype) -> None:
        self._file = file_
        self._memmap = np.memmap(
            file_, dtype, 'r+',
            data_offset(len(shape)),
            shape)

    def __array__(self) -> np.memmap:
        return self._memmap

    def extend(self, items: object) -> None:
        '''
        Extend the array by appending elements from `items`.
        '''
        # Convert `items` to a NumPy array.
        item_array = np.require(items, self._memmap.dtype, ['C_CONTIGUOUS'])

        # Raise an error if the arrays' shapes are not compatible.
        if self._memmap.ndim == 0:
            raise ValueError('scalars cannot be extended')
        if item_array.ndim == 0:
            raise ValueError('`items` must be a sequence')
        if item_array.shape[1:] != self._memmap.shape[1:]:
            raise ValueError('container and item shapes do not match')

        # Write data.
        self._file.seek(0, SEEK_END)
        self._file.write(item_array)
        self._file.flush()

        # Expand the memory-mapped array.
        dtype = self._memmap.dtype
        offset = data_offset(self._memmap.ndim)
        shape = (len(self._memmap) + len(item_array), *self._memmap.shape[1:])
        self._memmap = np.memmap(self._file, dtype, 'r+', offset, shape)

        # Overwrite the header.
        self._file.seek(0, SEEK_SET)
        with locking_header(self._file, LOCK_EX):
            self._file.write(ndarray_header(
                self._memmap.shape, self._memmap.dtype))
            self._file.flush()

    def append(self, item: object) -> None:
        '''
        Append `item` to the array.
        '''
        self.extend(np.asanyarray(item, self._memmap.dtype)[None])

class MemMapForwardingAttr:
    '''
    A descriptor that returns `obj._memmap.{key1}` when accessed.
    '''
    def __init__(self, key: str) -> None:
        self._key = key

    def __get__(self, obj: object, type_: type = None) -> Any:
        return getattr(getattr(obj, '_memmap'), self._key)

if 'sphinx' not in sys.modules:
    # Replace `PersistentArray` with `PersistentArrayImpl`
    # and add `np.memmap` methods and attribute-accessors.
    globals()['PersistentArray'] = PersistentArrayImpl
    for key in set(dir(np.memmap)) - set(dir(PersistentArrayImpl)):
        wrapper = MemMapForwardingAttr(key)
        wrapper.__doc__ = getattr(np.memmap, key).__doc__
        setattr(PersistentArrayImpl, key, wrapper)

#-- Reading --------------------------------------------------------------------

def read_cbor_file(path: Annotated[Path, '.cbor']) -> Any:
    '''
    Read a CBOR file.

    If the file encodes an indefinite-length array, a `PersistentList` will be
    returned.

    If the file encodes a 0–12-dimensional row-major array as specified in IETF
    RFC 8746, and the shape elements and byte string length are encoded as
    8-byte unsigned integers, a `PersistentArray` will be returned.

    Otherwise, a JSON-like object will be returned.
    '''
    # Defer to other readers if the path does not correspond to a CBOR file.
    if path.suffix != '.cbor':
        raise ValueError()

    # Open the specified file and read the first 128 bytes.
    f = cast(BufferedRandom, open(path, 'rb+'))
    with locking_header(f, LOCK_SH):
        header = cast(bytes, f.read(128))
    f.seek(0)

    # Try parsing the file as a `PersistentList`.
    try: return PersistentList(f, parse_list(header))
    except (ValueError, IndexError): pass

    # Try parsing the file as a `PersistentArray`.
    try: return PersistentArray(f, *parse_ndarray(header))
    except (ValueError, IndexError): pass

    # Parse the file using `cbor2`.
    return namespacify(cbor2.loads(f.read()))

def parse_list(buf: bytes) -> int:
    '''
    Parse the given buffer as the header of a `PersistentList` and return the
    number of items in the list.

    A `ValueError` is raised if an unexpected token is encountered and an
    `IndexError` is raised if the end of the buffer was reached while parsing.
    '''
    pos, size = parse_token(buf, 0, MAJOR_TYPE_ARRAY)
    fail_if(pos != 9)
    return size

def parse_ndarray(buf: bytes) -> Tuple[Tuple[int, ...], np.dtype]:
    '''
    Parse the given buffer as the header of a `PersistentArray` and return its
    shape and data type.

    A `ValueError` is raised if an unexpected token is encountered and an
    `IndexError` is raised if the end of the buffer was reached while parsing.
    '''
    # Check for a "multidimensional array" tag.
    pos, root_tag = parse_token(buf, 0, MAJOR_TYPE_TAG)
    fail_if(pos != 2 or root_tag != TAG_MULTIDIM_ARRAY)

    # Check whether the payload is a length-2 array.
    pos, root_len = parse_token(buf, pos, MAJOR_TYPE_ARRAY)
    fail_if(pos != 3 or root_len != 2)

    # Check for a shape array with up to 12 entries.
    pos, ndim = parse_token(buf, pos, MAJOR_TYPE_ARRAY)
    fail_if(pos != 4 or ndim > 12)

    # Read the shape array.
    shape = ndim * [0]
    for i in range(ndim):
        pos, shape[i] = parse_token(buf, pos, MAJOR_TYPE_UINT)
        fail_if(pos != 4 + 9 * (i + 1))

    # Check whether the shape array is followed by a typed data array.
    pos, dtype_tag = parse_token(buf, pos, MAJOR_TYPE_TAG)
    fail_if(pos != 6 + 9 * ndim or dtype_tag not in dtypes_by_tag)
    dtype = dtypes_by_tag[dtype_tag]

    # Check whether the data array is a byte string with an 8-byte size.
    pos, nbytes = parse_token(buf, pos, MAJOR_TYPE_BYTE_STRING)
    fail_if(pos != 15 + 9 * ndim or nbytes != np.prod(shape) * dtype.itemsize)

    # Return metadata if parsing succeeded.
    return tuple(shape), dtype

def parse_token(buf: bytes, pos: int,
                expected_major_type: int) -> Tuple[int, int]:
    '''
    Parse the CBOR token starting at `buf[pos]` and return the position of the
    next token in the buffer and the token's value.

    A `ValueError` is raised if the major type of the token does not match
    `expected_major_type`.
    '''
    major_type = buf[pos] & 0b1110_0000
    extra_info = buf[pos] & 0b0001_1111

    if major_type != expected_major_type:
        raise ValueError('CBOR parsing failed.')
    elif extra_info < INFO_NEXT_BYTE:
        return pos + 1, int(extra_info)
    elif extra_info == INFO_NEXT_BYTE:
        return pos + 2, int.from_bytes(buf[pos+1:pos+2], 'big')
    elif extra_info == INFO_NEXT_2_BYTES:
        return pos + 3, int.from_bytes(buf[pos+1:pos+3], 'big')
    elif extra_info == INFO_NEXT_4_BYTES:
        return pos + 5, int.from_bytes(buf[pos+1:pos+5], 'big')
    elif extra_info == INFO_NEXT_8_BYTES:
        return pos + 9, int.from_bytes(buf[pos+1:pos+9], 'big')
    else:
        raise ValueError('CBOR parsing failed.')

def fail_if(condition: bool) -> None:
    '''
    Raise a `ValueError` if the given condition is not true.
    '''
    if condition:
        raise ValueError('CBOR parsing failed')

#-- Writing --------------------------------------------------------------------

def write_object_as_cbor(path: Path, val: object) -> str:
    '''
    Write a JSON-encodable object or a NumPy array to a CBOR file.
    '''
    if isinstance(val, np.ndarray):
        write_ndarray(path, val)
    elif hasattr(val, '__array__'):
        write_ndarray(path, val.__array__()) # type: ignore
    elif isinstance(val, list):
        write_list(path, val)
    else:
        with open(path, 'wb') as f:
            cbor2.dump(dictify(val), f)
    return '.cbor'

def write_list(path: Path, list_: list) -> None:
    '''
    Write a list as a CBOR file containing an array.
    '''
    with open(path, 'wb') as f:
        f.write(list_header(len(list_)))
        for elem in list_:
            cbor2.dump(dictify(elem), f)
        f.flush()

def write_ndarray(path: Path, array: np.ndarray) -> None:
    '''
    Write an array as a CBOR file containing an row-major multidimensional
    array, as specified in IETF RFC 8746.

    The shape elements and the size of the byte string will be encoded as 8-byte
    unsigned integers.
    '''
    with open(path, 'wb') as f:
        f.write(ndarray_header(array.shape, array.dtype))
        f.write(np.ascontiguousarray(array).data)
        f.flush()

def list_header(length: int) -> bytes:
    '''
    Return the CBOR header for a list.
    '''
    return bytes((
        MAJOR_TYPE_ARRAY | INFO_NEXT_8_BYTES,
        *length.to_bytes(8, 'big')))

def ndarray_header(shape: Tuple[int, ...], dtype: np.dtype) -> bytes:
    '''
    Return the CBOR header for a multidimensional array.
    '''
    return bytes((
        MAJOR_TYPE_TAG | INFO_NEXT_BYTE,
        TAG_MULTIDIM_ARRAY,
        MAJOR_TYPE_ARRAY | 2,
        MAJOR_TYPE_ARRAY | len(shape),
        *chain.from_iterable(
            (MAJOR_TYPE_UINT | INFO_NEXT_8_BYTES,
             *n.to_bytes(8, 'big'))
            for n in shape),
        MAJOR_TYPE_TAG | INFO_NEXT_BYTE,
        tags_by_dtype[dtype],
        MAJOR_TYPE_BYTE_STRING | INFO_NEXT_8_BYTES,
        *int(np.prod(shape) * dtype.itemsize).to_bytes(8, 'big')))

def data_offset(ndim: int) -> int:
    '''
    Return the byte offset corresponding to the start of an `ndarray`'s data in
    a CBOR file.
    '''
    return 15 + 9 * ndim

@contextmanager
def locking_header(file_: BufferedRandom, mode: int) -> Iterator[None]:
    '''
    Return a context manager that acquires a lock on a CBOR file's header.
    '''
    if locking_is_supported:
        lockf(file_, mode, 128)
        yield
        lockf(file_, LOCK_UN)
    else:
        yield

max_stars_repo_name: MasonMcGill/artisan

max_issues_repo_name: MasonMcGill/artisan

max_forks_repo_name: MasonMcGill/artisan