prominences, left_bases, right_bases = signal.peak_prominences(x, peaks)

# === Rééchantillonnage ===

# Rééchantillonnage à nombre de points différent
num_samples_new = 500
x_resampled = signal.resample(x, num_samples_new)

# Avec temps associé
x_resampled, t_resampled = signal.resample(x, num_samples_new, t=t)

# Rééchantillonnage à fréquence différente
up = 3  # Facteur suréchantillonnage
down = 2  # Facteur sous-échantillonnage
x_resampled_poly = signal.resample_poly(x, up, down)

# Décimation (sous-échantillonnage avec anti-aliasing)
q = 4  # Facteur
x_decimated = signal.decimate(x, q)

# === Fenêtrage ===

# Appliquer fenêtre
windowed = x * window_hann[:len(x)]

# Détrending (élimination tendance)
x_detrended = signal.detrend(x, type='linear')  # Tendance linéaire
x_detrended_constant = signal.detrend(x, type='constant')  # Moyenne

# === Systèmes LTI (Linear Time-Invariant) ===

# Créer système depuis coefficients
system = signal.TransferFunction(b, a)

# Créer système depuis zéros, pôles, gain
zeros = [1, 2]
poles = [-1, -2, -3]
gain = 10
system_zpk = signal.ZerosPolesGain(zeros, poles, gain)

# Conversion vers autres formes
system_ba = system_zpk.to_tf()  # Vers b, a
system_ss = system_zpk.to_ss()  # Vers state-space

# State-space
A = np.array([[0, 1], [-1, -2]])
B = np.array([[0], [1]])
C = np.array([[1, 0]])
D = np.array([[0]])
system_ss = signal.StateSpace(A, B, C, D)

# Simulation système
t_sim = np.linspace(0, 5, 500)
u = np.sin(2 * np.pi * t_sim)  # Entrée
tout, yout, xout = signal.lsim(system, u, t_sim)

# Réponse impulsionnelle
t_impulse, y_impulse = signal.impulse(system)

# Réponse indicielle
t_step, y_step = signal.step(system)

# === Transformations bilinéaires ===

# Analogique vers digital (bilinear transform)
b_analog = [1]
a_analog = [1, 1]
b_digital, a_digital = signal.bilinear(b_analog, a_analog, fs=fs)

# Inverse
b_analog_back, a_analog_back = signal.bilinear_zpk([0], [-1], 1, fs=fs)

# === Corrélation et convolution ===

# Autocorrélation
autocorr = signal.correlate(x, x, mode='full')

# Corrélation croisée
crosscorr = signal.correlate(x, y, mode='full')

# Convolution
conv = signal.convolve(x, y, mode='full')
# mode: 'full', 'same', 'valid'

# Convolution 2D
kernel = np.ones((3, 3))
image = np.random.rand(100, 100)
conv2d = signal.convolve2d(image, kernel, mode='same')

# Corrélation 2D
corr2d = signal.correlate2d(image, kernel, mode='same')

# FFT convolution (plus rapide grandes données)
conv_fft = signal.fftconvolve(x, y, mode='same')

# Ordre arbitraire convolution
conv_oa = signal.oaconvolve(x, y, mode='same')

# === Hilbert transform ===

# Transformée de Hilbert (signal analytique)
analytic_signal = signal.hilbert(x)
amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase) / 
                          (2.0*np.pi) * fs)

# Hilbert 2D
hilbert2d = signal.hilbert2(image)

# === Enveloppe ===

# Enveloppe supérieure et inférieure
upper, lower = signal.hilbert(x), -signal.hilbert(x)

# === Lissage ===

# Moyennage mobile
window_size = 10
smoothed = np.convolve(x, np.ones(window_size)/window_size, mode='same')

# Lissage polynomial Savitzky-Golay
smoothed_sg = signal.savgol_filter(x, window_length=51, polyorder=3)

# === Détection d'événements ===

# Détection de fronts
threshold = 0.5
rising_edges = np.where(np.diff(x > threshold) > 0)[0]
falling_edges = np.where(np.diff(x > threshold) < 0)[0]

# === Conception filtre analogique ===

# Filtre prototype (normalized lowpass)
b_proto, a_proto = signal.butter(4, 1, analog=True)

# Transformation lowpass to highpass
b_hp, a_hp = signal.lp2hp(b_proto, a_proto, wo=10)

# Transformation lowpass to bandpass
b_bp, a_bp = signal.lp2bp(b_proto, a_proto, wo=10, bw=5)

# Transformation lowpass to bandstop
b_bs, a_bs = signal.lp2bs(b_proto, a_proto, wo=10, bw=5)

# Transformation lowpass to lowpass (mise à l'échelle)
b_lp, a_lp = signal.lp2lp(b_proto, a_proto, wo=10)

# === B-splines ===

# B-spline basis functions
n = 3  # Ordre
x_spline = np.linspace(0, 1, 100)
knots = [0, 0, 0, 0, 0.5, 1, 1, 1, 1]
basis = signal.bspline(x_spline, n)

# Interpolation B-spline
coeffs = [1, 2, 3, 4, 5]
y_spline = signal.cspline1d(coeffs)

# === Ordre de filtrage ===

# Ordre minimum pour specs données
N, Wn = signal.buttord(10, 20, 3, 40, fs=fs)

# === Analyse stabilité ===

# Vérifier stabilité (pôles dans cercle unité)
_, p, _ = signal.tf2zpk(b, a)
is_stable = np.all(np.abs(p) < 1)


[OK] SCIPY.NDIMAGE - TRAITEMENT D'IMAGES N-DIMENSIONNELLES


from scipy import ndimage
import numpy as np

# === Filtres de lissage ===

# Image test
image = np.random.rand(100, 100)

# Filtre gaussien
sigma = 2
smoothed = ndimage.gaussian_filter(image, sigma=sigma)

# Gaussien avec sigma différent par dimension
smoothed_aniso = ndimage.gaussian_filter(image, sigma=[1, 3])

# Filtre uniforme (moyennage)
size = 5
uniform = ndimage.uniform_filter(image, size=size)

# Filtre médian (élimine bruit impulsif)
median = ndimage.median_filter(image, size=5)

# Filtre percentile
percentile = ndimage.percentile_filter(image, percentile=50, size=5)

# Filtre minimum/maximum
minimum = ndimage.minimum_filter(image, size=5)
maximum = ndimage.maximum_filter(image, size=5)

# Filtre rang
rank = ndimage.rank_filter(image, rank=10, size=5)

# === Dérivées et gradients ===

# Gradient (Sobel)
gradient_x = ndimage.sobel(image, axis=0)
gradient_y = ndimage.sobel(image, axis=1)
gradient_magnitude = np.hypot(gradient_x, gradient_y)

# Prewitt
prewitt_x = ndimage.prewitt(image, axis=0)
prewitt_y = ndimage.prewitt(image, axis=1)

# Laplacien
laplacian = ndimage.laplace(image)

# Dérivée gaussienne
gaussian_gradient_x = ndimage.gaussian_gradient_magnitude(image, sigma=2)

# Laplacien de gaussienne (LoG)
log = ndimage.gaussian_laplace(image, sigma=2)

# Dérivée générique
deriv_x = ndimage.correlate1d(image, [-1, 0, 1], axis=0)

# === Morphologie mathématique ===

# Image binaire
binary = image > 0.5

# Érosion
eroded = ndimage.binary_erosion(binary)
eroded_iterations = ndimage.binary_erosion(binary, iterations=3)

# Dilatation
dilated = ndimage.binary_dilation(binary)

# Ouverture (érosion puis dilatation)
opened = ndimage.binary_opening(binary)

# Fermeture (dilatation puis érosion)
closed = ndimage.binary_closing(binary)

# Hit-or-miss
structure1 = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
structure2 = np.array([[1, 0, 1], [0, 0, 0], [1, 0, 1]])
hitmiss = ndimage.binary_hit_or_miss(binary, structure1, structure2)

# Remplissage trous
filled = ndimage.binary_fill_holes(binary)

# Propagation (flood fill)
mask = np.zeros_like(binary)
mask[50, 50] = 1
propagated = ndimage.binary_propagation(binary, mask=mask)

# Élément structurant personnalisé
struct = ndimage.generate_binary_structure(2, 1)  # 2D, connectivité 1
eroded_custom = ndimage.binary_erosion(binary, structure=struct)

# Morphologie niveaux de gris
gray_erosion = ndimage.grey_erosion(image, size=(5, 5))
gray_dilation = ndimage.grey_dilation(image, size=(5, 5))
gray_opening = ndimage.grey_opening(image, size=(5, 5))
gray_closing = ndimage.grey_closing(image, size=(5, 5))

# Gradient morphologique
morphological_gradient = ndimage.morphological_gradient(image, size=(3, 3))

# Top-hat blanc (différence avec ouverture)
white_tophat = ndimage.white_tophat(image, size=(5, 5))

# Top-hat noir (différence avec fermeture)
black_tophat = ndimage.black_tophat(image, size=(5, 5))

# === Segmentation ===

# Étiquetage composantes connexes
labeled, num_features = ndimage.label(binary)

# Avec structure personnalisée
struct = ndimage.generate_binary_structure(2, 2)  # Connectivité 2
labeled_8, num_features_8 = ndimage.label(binary, structure=struct)

# Trouver objets
objects = ndimage.find_objects(labeled)
# Retourne slices pour chaque objet

# Taille des objets
sizes = ndimage.sum(binary, labeled, range(num_features + 1))

# === Mesures sur régions ===

# Centre de masse
com = ndimage.center_of_mass(image)

# Centre de masse par région
coms = ndimage.center_of_mass(image, labeled, range(1, num_features + 1))

# Histogramme par région
hist = ndimage.histogram(image, 0, 1, 10, labels=labeled, index=1)

# Somme par région
sums = ndimage.sum(image, labeled, range(num_features + 1))

# Moyenne par région
means = ndimage.mean(image, labeled, range(1, num_features + 1))

# Variance
variances = ndimage.variance(image, labeled, range(1, num_features + 1))

# Écart-type
stds = ndimage.standard_deviation(image, labeled, range(1, num_features + 1))

# Minimum/Maximum
mins = ndimage.minimum(image, labeled, range(1, num_features + 1))
maxs = ndimage.maximum(image, labeled, range(1, num_features + 1))

# Position minimum/maximum
min_pos = ndimage.minimum_position(image, labeled, range(1, num_features + 1))
max_pos = ndimage.maximum_position(image, labeled, range(1, num_features + 1))

# Extrêmes
extrema = ndimage.extrema(image, labeled, range(1, num_features + 1))
# Retourne (min, max, min_pos, max_pos)

# === Transformations géométriques ===

# Rotation
rotated = ndimage.rotate(image, angle=45, reshape=True)
rotated_same = ndimage.rotate(image, angle=45, reshape=False)

# Avec interpolation spécifique
rotated_nearest = ndimage.rotate(image, 45, order=0)  # Nearest
rotated_bilinear = ndimage.rotate(image, 45, order=1)  # Bilinear
rotated_cubic = ndimage.rotate(image, 45, order=3)  # Cubic

# Zoom
zoomed = ndimage.zoom(image, zoom=2)  # Agrandir 2x
zoomed_aniso = ndimage.zoom(image, zoom=[2, 1])  # Anisotrope

# Shift (translation)
shifted = ndimage.shift(image, shift=[10, 5])

# Transformation affine générale
matrix = np.array([[1.2, 0.1], [0.1, 0.8]])
offset = [5, 10]
affined = ndimage.affine_transform(image, matrix, offset=offset)

# Coordonnées de sortie spécifiques
output_shape = (150, 150)
affined_large = ndimage.affine_transform(image, matrix, 
                                        output_shape=output_shape)

# === Interpolation ===

# Grille de coordonnées
coords = np.meshgrid(np.arange(50), np.arange(50))

# Map coordinates (warping arbitraire)
warped = ndimage.map_coordinates(image, coords, order=3)

# Spline filter (prétraitement pour interpolation)
spline_filtered = ndimage.spline_filter(image, order=3)
warped_spline = ndimage.map_coordinates(spline_filtered, coords, 
                                       prefilter=False)

# === Convolution et corrélation ===

# Kernel personnalisé
kernel = np.array([[0, -1, 0], [-1, 5, -1], [0, -1, 0]])

# Convolution
convolved = ndimage.convolve(image, kernel)

# Corrélation
correlated = ndimage.correlate(image, kernel)

# Avec mode spécifique
# mode: 'reflect', 'constant', 'nearest', 'mirror', 'wrap'
convolved_reflect = ndimage.convolve(image, kernel, mode='reflect')

# Avec valeur constante
convolved_const = ndimage.convolve(image, kernel, mode='constant', cval=0)

# Convolution 1D
kernel_1d = np.array([1, 2, 1]) / 4
conv_1d = ndimage.convolve1d(image, kernel_1d, axis=0)

# Corrélation 1D
corr_1d = ndimage.correlate1d(image, kernel_1d, axis=0)

# === Filtre générique ===

# Fonction personnalisée
def custom_filter(values):
    return np.sum(values ** 2)

filtered_custom = ndimage.generic_filter(image, custom_filter, size=3)

# Filtre générique 1D
filtered_1d = ndimage.generic_filter1d(image, custom_filter, 
                                      filter_size=5, axis=0)

# === Distance transform ===

# Transform distance euclidienne
distances = ndimage.distance_transform_edt(binary)

# Avec indices des plus proches
distances, indices = ndimage.distance_transform_edt(binary, 
                                                   return_indices=True)

# Distance Manhattan
distances_manhattan = ndimage.distance_transform_cdt(binary, metric='taxicab')

# Distance Chebyshev
distances_chebyshev = ndimage.distance_transform_cdt(binary, 
                                                    metric='chessboard')

# Distance avec fonction personnalisée
distances_bf = ndimage.distance_transform_bf(binary, metric='euclidean')

# === Watershed ===

# Marqueurs
markers = np.zeros_like(image, dtype=int)
markers[20, 20] = 1
markers[80, 80] = 2

# Watershed segmentation
from scipy.ndimage import watershed_ift
segmented = watershed_ift(image.astype(np.uint8), markers)

# === Ligne de partage des eaux (morphologique) ===

# Maxima locaux comme marqueurs
local_max = ndimage.maximum_filter(image, size=20) == image
markers, num_markers = ndimage.label(local_max)

# === Fourier ===

# FFT avec décalage
fourier = ndimage.fourier_shift(image, shift=[5, 10])

# Filtre gaussien dans domaine de Fourier
fourier_gaussian = ndimage.fourier_gaussian(image, sigma=2)

# Filtre uniforme dans Fourier
fourier_uniform = ndimage.fourier_uniform(image, size=5)

# Filtre ellipsoïde dans Fourier
fourier_ellipsoid = ndimage.fourier_ellipsoid(image, size=5)

# === Interpolation de voxels ===

# Pour images 3D
volume = np.random.rand(50, 50, 50)

# Rotation 3D
rotated_3d = ndimage.rotate(volume, angle=45, axes=(0, 1))

# Zoom 3D
zoomed_3d = ndimage.zoom(volume, zoom=2)

# === Utilitaires ===

# Générer grille de coordonnées
coords = np.indices(image.shape)

# Élément structurant
struct_2d_4 = ndimage.generate_binary_structure(2, 1)  # 4-connectivité
struct_2d_8 = ndimage.generate_binary_structure(2, 2)  # 8-connectivité
struct_3d_6 = ndimage.generate_binary_structure(3, 1)  # 6-connectivité

# Itérer structure
struct_iter = ndimage.iterate_structure(struct_2d_4, iterations=3)


[OK] SCIPY.SPATIAL - GÉOMÉTRIE SPATIALE ET ALGORITHMES


from scipy import spatial
import numpy as np

# === Distance entre points ===

# Points 2D
points = np.random.rand(10, 2)

# Matrice de distances euclidiennes
dist_matrix = spatial.distance_matrix(points, points)

# Distance entre deux points
p1 = np.array([0, 0])
p2 = np.array([3, 4])
dist = spatial.distance.euclidean(p1, p2)  # 5.0

# Autres distances
dist_manhattan = spatial.distance.cityblock(p1, p2)
dist_chebyshev = spatial.distance.chebyshev(p1, p2)
dist_minkowski = spatial.distance.minkowski(p1, p2, p=3)
dist_cosine = spatial.distance.cosine(p1, p2)
dist_correlation = spatial.distance.correlation(p1, p2)
dist_hamming = spatial.distance.hamming([1, 0, 1], [1, 1, 0])
dist_jaccard = spatial.distance.jaccard([1, 0, 1], [1, 1, 0])

# Distances condensées (forme triangulaire)
dists_condensed = spatial.distance.pdist(points, metric='euclidean')
# Métriques: 'euclidean', 'cityblock', 'cosine', 'correlation', etc.

# Convertir vers forme carrée
dists_square = spatial.distance.squareform(dists_condensed)

# Distance avec fonction personnalisée
def custom_metric(u, v):
    return np.sum(np.abs(u - v) ** 1.5)

dists_custom = spatial.distance.pdist(points, metric=custom_metric)

# Distances entre deux ensembles
set1 = np.random.rand(5, 2)
set2 = np.random.rand(8, 2)
dists_cdist = spatial.distance.cdist(set1, set2, metric='euclidean')

# === KDTree - Recherche de voisins ===

# Construction
points = np.random.rand(1000, 3)
tree = spatial.KDTree(points)

# Plus proche voisin
query_point = np.array([0.5, 0.5, 0.5])
dist, idx = tree.query(query_point)

# K plus proches voisins
k = 5
dists, indices = tree.query(query_point, k=k)

# Voisins dans rayon
radius = 0.1
indices_radius = tree.query_ball_point(query_point, radius)

# Paires de points à distance < r
pairs = tree.query_pairs(r=0.05)

# Compter voisins dans rayon
counts = tree.query_ball_tree(tree, r=0.05)

# === cKDTree - Version optimisée (C++) ===

tree_fast = spatial.cKDTree(points)

# Mêmes méthodes mais plus rapides
dists, indices = tree_fast.query(query_point, k=k)

# Requêtes multiples (parallélisées)
query_points = np.random.rand(100, 3)
dists, indices = tree_fast.query(query_points, k=k, workers=-1)

# === Triangulation de Delaunay ===

# Points 2D
points_2d = np.random.rand(30, 2)

# Triangulation
tri = spatial.Delaunay(points_2d)

# Simplexes (triangles)
triangles = tri.simplices  # Indices des sommets

# Trouver simplex contenant point
point_test = np.array([0.5, 0.5])
simplex_idx = tri.find_simplex(point_test)

# Voisins de chaque simplex
neighbors = tri.neighbors

# Équations des hyperplans
equations = tri.equations

# Transformation barycentrique
transform = tri.transform

# === Diagramme de Voronoï ===

vor = spatial.Voronoi(points_2d)

# Sommets des régions
vertices = vor.vertices

# Régions (indices des sommets)
regions = vor.regions

# Points par région
point_region = vor.point_region

# Arêtes (ridge)
ridge_points = vor.ridge_points
ridge_vertices = vor.ridge_vertices

# === Enveloppe convexe (Convex Hull) ===

hull = spatial.ConvexHull(points_2d)

# Simplexes de l'enveloppe
hull_simplices = hull.simplices

# Points sur l'enveloppe
hull_points = points_2d[hull.vertices]

# Volume (aire en 2D)
area = hull.volume

# Périmètre (2D) ou surface (3D)
perimeter = hull.area

# Équations des facettes
equations_hull = hull.equations

# === Half-space intersection ===

# Demi-espaces: Ax <= b
halfspaces = np.array([[1, 0, 1],   # x <= 1
                       [-1, 0, 1],  # -x <= 1  =>  x >= -1
                       [0, 1, 1],   # y <= 1
                       [0, -1, 1]]) # y >= -1

feasible_point = np.array([0, 0])
hs = spatial.HalfspaceIntersection(halfspaces, feasible_point)

# Sommets de l'intersection
intersection_vertices = hs.intersections

# === Triangulation sphérique ===

# Points 3D sur sphère
phi = np.random.uniform(0, 2*np.pi, 50)
theta = np.random.uniform(0, np.pi, 50)
x = np.sin(theta) * np.cos(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(theta)
points_sphere = np.column_stack([x, y, z])

sphere_tri = spatial.SphericalVoronoi(points_sphere)

# === Procrustes analysis ===

# Alignement optimal de deux ensembles de points
points_A = np.random.rand(10, 2)
points_B = points_A + np.random.randn(10, 2) * 0.1

mtx1, mtx2, disparity = spatial.procrustes(points_A, points_B)
# mtx1, mtx2: données alignées
# disparity: mesure de dissimilarité

# === Transformations géométriques ===

# Rotation
from spatial.transform import Rotation

# Depuis angles d'Euler
r = Rotation.from_euler('xyz', [90, 0, 0], degrees=True)
rotation_matrix = r.as_matrix()

# Depuis quaternion
quat = [0, 0, np.sin(np.pi/4), np.cos(np.pi/4)]
r_quat = Rotation.from_quat(quat)

# Depuis vecteur rotation (axis-angle)
rotvec = [0, 0, np.pi/2]
r_rotvec = Rotation.from_rotvec(rotvec)

# Depuis matrice
r_matrix = Rotation.from_matrix(rotation_matrix)

# Appliquer rotation
points_3d = np.random.rand(10, 3)
rotated_points = r.apply(points_3d)

# Composition de rotations
r1 = Rotation.from_euler('z', 45, degrees=True)
r2 = Rotation.from_euler('x', 30, degrees=True)
r_combined = r2 * r1  # Ordre important!

# Rotation inverse
r_inv = r.inv()

# Interpolation de rotations (Slerp)
key_times = [0, 1, 2]
key_rots = Rotation.from_euler('z', [0, 45, 90], degrees=True)
from scipy.spatial.transform import Slerp
slerp = Slerp(key_times, key_rots)
interp_times = np.linspace(0, 2, 10)
interp_rots = slerp(interp_times)

# === Géométrie computationnelle ===

# Test si point dans enveloppe convexe
in_hull = hull.find_simplex(point_test) >= 0

# Plus proche point sur enveloppe convexe
# (nécessite calculs supplémentaires)

# === QHull options ===

# Options avancées pour triangulation
tri_qhull = spatial.Delaunay(points_2d, qhull_options='QJ')  # Joggle input
hull_qhull = spatial.ConvexHull(points_2d, qhull_options='QJ')


[OK] SCIPY.CLUSTER - CLUSTERING ET ANALYSE DE GROUPES


from scipy.cluster import vq, hierarchy
from scipy.spatial.distance import pdist
import numpy as np

# === K-means (Vector Quantization) ===

# Données
data = np.random.rand(100, 2)

# K-means standard
k = 3
centroids, distortion = vq.kmeans(data, k)
# centroids: centres des clusters
# distortion: distortion moyenne

# Avec itérations min
centroids, dist = vq.kmeans(data, k, iter=50)

# Assigner points aux clusters
idx, dists = vq.vq(data, centroids)
# idx: index du centroid le plus proche pour chaque point
# dists: distances correspondantes

# K-means++  initialization
centroids_pp = vq.kmeans2(data, k, minit='points')[0]

# K-means2 (plus d'options)
centroids2, labels = vq.kmeans2(data, k, iter=100, minit='++')

# Whitening (normalisation)
whitened = vq.whiten(data)
centroids_w, _ = vq.kmeans(whitened, k)

# === Clustering hiérarchique ===

# Matrice de distances
distances = pdist(data, metric='euclidean')

# Linkage (construction hiérarchie)
Z = hierarchy.linkage(distances, method='ward')
# Méthodes: 'single', 'complete', 'average', 'weighted', 'centroid',
#          'median', 'ward'

# Alternative: linkage direct sur données
Z_direct = hierarchy.linkage(data, method='ward', metric='euclidean')

# Dendrogramme
dendrogram = hierarchy.dendrogram(Z)

# Avec options
dendrogram_fancy = hierarchy.dendrogram(
    Z,
    truncate_mode='lastp',  # Montrer seulement p derniers merged clusters
    p=12,
    leaf_rotation=90,
    leaf_font_size=10,
    show_contracted=True
)

# Couper dendrogramme (former clusters)
clusters = hierarchy.fcluster(Z, t=3, criterion='maxclust')  # 3 clusters
clusters_dist = hierarchy.fcluster(Z, t=1.5, criterion='distance')  # Seuil distance

# Critères disponibles:
# 'inconsistent', 'distance', 'maxclust', 'monocrit', 'maxclust_monocrit'

# Coefficient de cophenetic correlation
c, coph_dists = hierarchy.cophenet(Z, distances)
# Mesure qualité du clustering

# Matrice de distances cophenétiques
coph_matrix = hierarchy.cophenet(Z)

# Leader clustering (single-pass)
T = hierarchy.leaders(Z, clusters)

# Vérifier si linkage est valide
is_valid = hierarchy.is_valid_linkage(Z)

# Vérifier monotonie
is_monotonic = hierarchy.is_monotonic(Z)

# Calculer inconsistency
inconsistency = hierarchy.inconsistent(Z, d=2)

# Maxdists et maxinconsts
maxdists = hierarchy.maxdists(Z)
maxinconsts = hierarchy.maxinconsts(Z)

# Nombre de clusters à chaque étape
num_clusters = hierarchy.fclusterdata(data, t=3, criterion='maxclust', 
                                     method='ward')

# === Méthodes de linkage ===

# Single linkage (nearest neighbor)
Z_single = hierarchy.single(distances)

# Complete linkage (farthest neighbor)
Z_complete = hierarchy.complete(distances)

# Average linkage
Z_average = hierarchy.average(distances)

# Weighted average
Z_weighted = hierarchy.weighted(distances)

# Centroid linkage
Z_centroid = hierarchy.centroid(distances)

# Median linkage
Z_median = hierarchy.median(distances)

# Ward (minimize variance)
Z_ward = hierarchy.ward(distances)

# === Conversion et manipulation ===

# Linkage vers tree
rootnode, nodelist = hierarchy.to_tree(Z, rd=True)

# Optimal leaf ordering (minimize distance between adjacent leaves)
Z_optimal = hierarchy.optimal_leaf_ordering(Z, distances)

# Extraction de clusters à différentes hauteurs
clusters_multi = hierarchy.cut_tree(Z, n_clusters=[2, 3, 4])

# === Distances entre clusters ===

# Correspondance entre clusters
from scipy.cluster.hierarchy import correspond
# Comparer deux partitions

# === Validation ===

# Silhouette score (nécessite sklearn)
# from sklearn.metrics import silhouette_score
# score = silhouette_score(data, clusters)


[OK] SCIPY.FFT - TRANSFORMÉE DE FOURIER RAPIDE


from scipy import fft
import numpy as np

# === FFT 1D ===

# Signal
N = 1000
T = 1.0 / 800.0  # Période échantillonnage
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)

# FFT
yf = fft.fft(y)

# Fréquences
xf = fft.fftfreq(N, T)

# Spectre de puissance
power = 2.0/N * np.abs(yf[0:N//2])

# === FFT inverse ===

y_reconstructed = fft.ifft(yf)

# Vérification
assert np.allclose(y, y_reconstructed.real)

# === FFT réelle (optimisée pour signaux réels) ===

# FFT réelle (retourne seulement fréquences positives)
yf_real = fft.rfft(y)
xf_real = fft.rfftfreq(N, T)

# IFFT réelle
y_reconstructed_real = fft.irfft(yf_real)

# === DCT - Discrete Cosine Transform ===

# DCT type II (défaut, utilisée en JPEG)
y_dct = fft.dct(y)

# Types disponibles: 1, 2, 3, 4
y_dct_1 = fft.dct(y, type=1)
y_dct_2 = fft.dct(y, type=2)  # Défaut
y_dct_3 = fft.dct(y, type=3)
y_dct_4 = fft.dct(y, type=4)

# IDCT
y_idct = fft.idct(y_dct)

# Normalisation
y_dct_norm = fft.dct(y, norm='ortho')

# === DST - Discrete Sine Transform ===

y_dst = fft.dst(y)
y_idst = fft.idst(y_dst)

# Types 1, 2, 3, 4
y_dst_2 = fft.dst(y, type=2)

# === FFT 2D ===

# Image
image = np.random.rand(100, 100)

# FFT 2D
image_fft = fft.fft2(image)

# IFFT 2D
image_reconstructed = fft.ifft2(image_fft)

# FFT réelle 2D
image_fft_real = fft.rfft2(image)
image_reconstructed_real = fft.irfft2(image_fft_real)

# Fréquences 2D
freqs_x = fft.fftfreq(image.shape[0])
freqs_y = fft.fftfreq(image.shape[1])

# === FFT N-dimensionnelle ===

# Volume 3D
volume = np.random.rand(50, 50, 50)

# FFT ND
volume_fft = fft.fftn(volume)

# IFFT ND
volume_reconstructed = fft.ifftn(volume_fft)

# FFT réelle ND
volume_fft_real = fft.rfftn(volume)
volume_reconstructed_real = fft.irfftn(volume_fft_real)

# Sur axes spécifiques
volume_fft_axes = fft.fftn(volume, axes=(0, 1))  # FFT sur 2 premiers axes

# === FFT avec zero-padding ===

# Padding pour puissance de 2 (plus rapide)
n_padded = 2048
yf_padded = fft.fft(y, n=n_padded)

# === Shift FFT (centrer fréquence zéro) ===

yf_shifted = fft.fftshift(yf)
xf_shifted = fft.fftshift(xf)

# Inverse shift
yf_unshifted = fft.ifftshift(yf_shifted)

# === DCT/DST 2D et ND ===

# DCT 2D
image_dct = fft.dctn(image)
image_idct = fft.idctn(image_dct)

# DST 2D
image_dst = fft.dstn(image)
image_idst = fft.idstn(image_dst)

# === Convolution via FFT ===

# Deux signaux
signal1 = np.random.rand(100)
signal2 = np.random.rand(100)

# Convolution via FFT
fft1 = fft.fft(signal1)
fft2 = fft.fft(signal2)
conv_fft = fft.ifft(fft1 * fft2).real

# === Filtrage fréquentiel ===

# Filtre passe-bas
cutoff = 20  # Hz
filter_mask = np.abs(xf) < cutoff
yf_filtered = yf * filter_mask
y_filtered = fft.ifft(yf_filtered).real

# === FFT multi-threading ===

# Planification pour optimisation
# (scipy.fft utilise automatiquement FFTW si disponible)

# === Fonctions utilitaires ===

# Fréquences pour FFT réelle
freqs_rfft = fft.rfftfreq(N, T)

# Next fast length (pour optimisation)
next_fast = fft.next_fast_len(N)

# === Normalisation ===

# Modes de normalisation
# norm=None (défaut): pas de normalisation
# norm='ortho': normalisation orthogonale
# norm='forward': normalisation sur FFT
# norm='backward': normalisation sur IFFT

yf_ortho = fft.fft(y, norm='ortho')
y_ortho = fft.ifft(yf_ortho, norm='ortho')


[OK] SCIPY.IO - ENTRÉES/SORTIES FICHIERS


from scipy import io
import numpy as np

# === Fichiers MATLAB ===

# Sauvegarder variables MATLAB
data = {'variable1': np.array([1, 2, 3]),
        'variable2': np.random.rand(5, 5),
        'text': 'Hello MATLAB'}
io.savemat('data.mat', data)

# Charger fichier MATLAB
mat_contents = io.loadmat('data.mat')
var1 = mat_contents['variable1']

# Options de sauvegarde
io.savemat('data.mat', data, 
           do_compression=True,    # Compression
           format='5',             # Version MATLAB (4, 5)
           oned_as='column')       # 1D arrays as column vectors

# Charger avec options
mat_contents = io.loadmat('data.mat',
                          squeeze_me=True,      # Squeeze unit dimensions
                          struct_as_record=False)  # Structs as objects

# Who mat (lister variables sans charger)
variables = io.whosmat('data.mat')

# === Fichiers NetCDF ===

# Lire NetCDF
from scipy.io import netcdf_file

# Ouvrir fichier
with netcdf_file('data.nc', 'r') as f:
    # Dimensions
    print(f.dimensions)
    
    # Variables
    var = f.variables['temperature']
    data = var[:].copy()
    
    # Attributs
    attrs = var._attributes

# Créer fichier NetCDF
with netcdf_file('output.nc', 'w') as f:
    # Créer dimensions
    f.createDimension('time', 10)
    f.createDimension('lat', 5)
    f.createDimension('lon', 5)
    
    # Créer variable
    temp = f.createVariable('temperature', 'f', ('time', 'lat', 'lon'))
    temp[:] = np.random.rand(10, 5, 5)
    temp.units = 'celsius'

# === Fichiers WAV (audio) ===

# Lire fichier WAV
from scipy.io import wavfile

samplerate, data = wavfile.read('audio.wav')

# Écrire fichier WAV
wavfile.write('output.wav', samplerate, data)

# Avec données float
data_float = np.random.rand(44100) * 2 - 1  # [-1, 1]
wavfile.write('output_float.wav', 44100, data_float.astype(np.float32))

# === Fichiers ARFF (Weka) ===

from scipy.io import arff

# Charger ARFF
with open('data.arff') as f:
    data, meta = arff.loadarff(f)

# data: numpy structured array
# meta: metadata (attributes, relation name)

# === Fichiers Matrix Market ===

# Format sparse matrix
from scipy.io import mmread, mmwrite

# Lire
sparse_matrix = mmread('matrix.mtx')

# Écrire
from scipy.sparse import csr_matrix
sparse = csr_matrix([[1, 0, 2], [0, 0, 3], [4, 5, 6]])
mmwrite('output.mtx', sparse)

# === Fichiers IDL ===

# Lire fichiers IDL save
from scipy.io import readsav
data = readsav('data.sav')

# === Fichiers Harwell-Boeing ===

from scipy.io.harwell_boeing import hb_read, hb_write

# Lire
sparse = hb_read('matrix.hb')

# Écrire
hb_write('output.hb', sparse)

# === Utilitaires ===

# Vérifier si fichier MATLAB est valide
is_mat = io.matlab.whosmat('data.mat')


[OK] SCIPY.MISC - UTILITAIRES DIVERS


from scipy import misc
import numpy as np

# Note: scipy.misc est deprecated pour images
# Utiliser imageio ou PIL/Pillow à la place

# === Images de test (deprecated, utiliser scipy.datasets) ===

# Ancienne méthode (deprecated)
# face = misc.face()  # Image test

# Nouvelle méthode
from scipy.datasets import face, ascent, electrocardiogram

# Image visage
face_img = face()  # Array (768, 1024, 3)

# Image escalier
ascent_img = ascent()  # Array (512, 512)

# ECG signal
ecg = electrocardiogram()

# === Dérivées numériques ===

# Central difference
from scipy.misc import central_diff_weights
weights = central_diff_weights(3)  # Weights pour dérivée ordre 1

# Derivative (utiliser scipy.optimize.approx_fprime plutôt)


[OK] SCIPY.SPARSE.LINALG - ALGÈBRE LINÉAIRE CREUSE (SUITE)


# Voir section SCIPY.SPARSE pour détails complets

from scipy.sparse import linalg as sparse_linalg
from scipy import sparse
import numpy as np

# === Solveurs itératifs avancés ===

# GCROTMK
A = sparse.random(100, 100, density=0.1, format='csr')
b = np.random.rand(100)
x, info = sparse_linalg.gcrotmk(A, b)

# QMR - Quasi-Minimal Residual
x, info = sparse_linalg.qmr(A, b)

# TFQMR - Transpose-Free QMR
x, info = sparse_linalg.tfqmr(A, b)

# === Eigenvalue problems ===

# LOBPCG - Locally Optimal Block Preconditioned Conjugate Gradient
X = np.random.rand(100, 3)  # Initial guess
eigenvalues, eigenvectors = sparse_linalg.lobpcg(A, X, largest=True)

# === Interface ARPACK ===

# Plus flexible que eigs/eigsh
from scipy.sparse.linalg import ArpackNoConvergence

try:
    eigenvalues, eigenvectors = sparse_linalg.eigs(A, k=6, which='LM')
except ArpackNoConvergence as e:
    eigenvalues = e.eigenvalues
    eigenvectors = e.eigenvectors


[OK] FONCTIONNALITÉS AVANCÉES ET BONNES PRATIQUES


# === Vectorisation et performance ===

# Toujours utiliser opérations vectorisées NumPy
x = np.linspace(0, 10, 1000000)

# Lent
# result = [np.sin(xi) for xi in x]

# Rapide
result = np.sin(x)

# === Gestion mémoire matrices creuses ===

# Choisir bon format
# - COO: construction
# - CSR: opérations ligne, produit matrice-vecteur
# - CSC: opérations colonne, solveurs
# - LIL/DOK: construction incrémentale

# Conversion explicite pour performance
lil = sparse.lil_matrix((1000, 1000))
# ... remplissage ...
csr = lil.tocsr()  # Convertir avant calculs

# === Parallélisation ===

# cKDTree avec workers
from scipy.spatial import cKDTree
tree = cKDTree(points)
dists, indices = tree.query(query_points, k=5, workers=-1)  # Tous les CPU

# FFT multithreaded (automatique avec FFTW)
from scipy import fft
result = fft.fft(large_signal)  # Utilise threads automatiquement

# === Précision numérique ===

# Utiliser forme SOS pour filtres IIR (plus stable)
from scipy import signal
sos = signal.butter(10, 0.1, output='sos')
y = signal.sosfilt(sos, x)

# Éviter coefficients b, a pour ordres élevés
# b, a = signal.butter(10, 0.1)  # Peut être instable

# === Tests statistiques - corrections multiples ===

from scipy import stats

# Correction Bonferroni
p_values = [0.01, 0.04, 0.03, 0.05]
alpha = 0.05
bonferroni_alpha = alpha / len(p_values)
significant = [p < bonferroni_alpha for p in p_values]

# FDR - False Discovery Rate (Benjamini-Hochberg)
# from statsmodels.stats.multitest import multipletests
# reject, pvals_corrected, _, _ = multipletests(p_values, method='fdr_bh')

# === Optimisation - Callbacks et monitoring ===

from scipy import optimize

history = {'x': [], 'fun': []}

def callback(xk):
    history['x'].append(xk.copy())
    history['fun'].append(objective(xk))
    print(f"Iteration {len(history['x'])}: f={history['fun'][-1]:.6f}")

def objective(x):
    return np.sum(x**2)

result = optimize.minimize(objective, [1, 1], callback=callback)

# === Intégration ODE - Gestion événements ===

from scipy import integrate

def pendulum(t, y, g, L):
    theta, omega = y
    return [omega, -g/L * np.sin(theta)]

def apex_event(t, y, g, L):
    return y[1]  # omega = 0

apex_event.terminal = False  # Continue après événement
apex_event.direction = -1    # Détecte seulement passages décroissants

sol = integrate.solve_ivp(
    pendulum, 
    [0, 10], 
    [np.pi/4, 0],
    args=(9.81, 1.0),
    events=apex_event,
    dense_output=True
)

# Temps des apex
apex_times = sol.t_events[0]

# === Interpolation - Extrapolation prudente ===

from scipy import interpolate

x = np.linspace(0, 10, 11)
y = np.sin(x)

# Avec avertissement hors limites
f = interpolate.interp1d(x, y, bounds_error=True)
# f(11)  # Erreur!

# Extrapolation linéaire
f_extrap = interpolate.interp1d(x, y, fill_value='extrapolate')
y_extrap = f_extrap(11)  # OK mais prudence

# === Sparse matrices - Préallocation ===

# Mauvais: construction incrémentale avec CSR
# A = sparse.csr_matrix((1000, 1000))
# for i in range(1000):
#     A[i, i] = 1  # Très lent!

# Bon: utiliser LIL ou COO puis convertir
A_lil = sparse.lil_matrix((1000, 1000))
for i in range(1000):
    A_lil[i, i] = 1
A = A_lil.tocsr()

# Meilleur: construction directe COO
row = np.arange(1000)
col = np.arange(1000)
data = np.ones(1000)
A = sparse.coo_matrix((data, (row, col)), shape=(1000, 1000)).tocsr()

# === Distribution fitting - AIC/BIC ===

from scipy import stats

data = np.random.normal(0, 1, 1000)

# Ajuster et calculer AIC
params = stats.norm.fit(data)
log_likelihood = np.sum(stats.norm.logpdf(data, *params))
k = len(params)
n = len(data)

aic = 2*k - 2*log_likelihood
bic = k*np.log(n) - 2*log_likelihood

# Comparer distributions
best_aic = np.inf
best_dist = None

for dist_name in ['norm', 'expon', 'gamma']:
    dist = getattr(stats, dist_name)
    params = dist.fit(data)
    ll = np.sum(dist.logpdf(data, *params))
    aic_current = 2*len(params) - 2*ll
    
    if aic_current < best_aic:
        best_aic = aic_current
        best_dist = dist_name

print(f"Meilleure distribution: {best_dist}")

# === Gestion erreurs numériques ===

import warnings
from scipy import linalg

# Supprimer warnings temporairement
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    result = linalg.inv(nearly_singular_matrix)

# Vérifier convergence
from scipy.optimize import OptimizeResult

result = optimize.minimize(objective, x0)
if not result.success:
    print(f"Échec: {result.message}")
    print(f"Nombre itérations: {result.nit}")


[OK] RESSOURCES ET RÉFÉRENCES


# Documentation officielle
# https://docs.scipy.org/doc/scipy/

# Tutoriels
# https://docs.scipy.org/doc/scipy/tutorial/

# API Reference
# https://docs.scipy.org/doc/scipy/reference/

# Cookbook (exemples)
# https://scipy-cookbook.readthedocs.io/

# GitHub
# https://github.com/scipy/scipy

# Mailing lists
# https://scipy.org/community/

# Livres recommandés:
# - "Python for Data Analysis" - Wes McKinney
# - "Elegant SciPy" - Juan Nunez-Iglesias et al.
# - "Scientific Computing with Python" - Various

# Citation SciPy:
# Virtanen, P. et al. (2020). SciPy 1.0: fundamental algorithms for 
# scientific computing in Python. Nature Methods, 17(3), 261-272.


[OK] VERSION ET COMPATIBILITÉ


import scipy
import numpy as np

print(f"SciPy version: {scipy.__version__}")
print(f"NumPy version: {np.__version__}")

# Vérifier disponibilité fonctionnalités
import sys
print(f"Python version: {sys.version}")

# Compiler info
scipy.show_config()

# Dépendances BLAS/LAPACK
np.show_config()


# === FIN DU CHEATSHEET SCIPY ===
# Dernière mise à jour: Compatible avec SciPy 1.11+
# Pour mises à jour: consulter https://docs.scipy.org/doc/scipy/release.html
maxA_op = sparse_linalg.LinearOperator((3, 3), matvec=matvec, rmatvec=rmatvec)

# Utilisation avec solveurs
x, info = sparse_linalg.gmres(A_op, b)

# Exponentielle matricielle sparse
expm_A = sparse_linalg.expm(A)

# Multiplication exp(A) @ v (efficace)
v = np.array([1, 1, 1])
result = sparse_linalg.expm_multiply(A, v)

# === Graphe de la matrice ===

# Élimination minimum degree
perm = sparse.csgraph.reverse_cuthill_mckee(A)
A_reordered = A[perm, :][:, perm]

# === Performances ===

# Élimination zéros
A_sparse.eliminate_zeros()

# Somme doublons
A_sparse.sum_duplicates()

# Trier indices
A_sparse.sort_indices()

# Vérifier format
A_sparse.check_format()

# === Sauvegarde/Chargement ===

# Sauvegarder
sparse.save_npz('matrix.npz', A_sparse)

# Charger
A_loaded = sparse.load_npz('matrix.npz')


[OK] SCIPY.STATS - STATISTIQUES ET DISTRIBUTIONS


from scipy import stats
import numpy as np

# === Distributions continues ===

# Normale (Gaussienne)
mu, sigma = 0, 1
normal = stats.norm(loc=mu, scale=sigma)

# PDF - Probabilité density function
x = np.linspace(-3, 3, 100)
pdf = normal.pdf(x)
pdf_point = stats.norm.pdf(0, loc=0, scale=1)

# CDF - Cumulative distribution function
cdf = normal.cdf(x)
cdf_point = stats.norm.cdf(1.96)  # ≈ 0.975

# Quantile (inverse CDF)
q = normal.ppf(0.975)  # ≈ 1.96

# Survival function (1 - CDF)
sf = normal.sf(x)

# Inverse survival function
isf = normal.isf(0.025)  # ≈ 1.96

# Génération aléatoire
samples = normal.rvs(size=1000)

# Moments
mean = normal.mean()
variance = normal.var()
std = normal.std()
skewness = normal.stats(moments='s')
kurtosis = normal.stats(moments='k')
mean, var, skew, kurt = normal.stats(moments='mvsk')

# Intervalle contenant x% de la probabilité
interval = normal.interval(0.95)  # (-1.96, 1.96)

# Support (domaine)
a, b = normal.support()

# Entropie
entropy = normal.entropy()

# Fonction de survie médiane
median = normal.median()

# === Autres distributions continues ===

# Uniforme
uniform = stats.uniform(loc=0, scale=10)  # [0, 10)

# Exponentielle
exponential = stats.expon(scale=1/2)  # λ=2

# Gamma
gamma = stats.gamma(a=2, scale=2)  # a=shape, scale=1/rate

# Beta
beta = stats.beta(a=2, b=5)

# Chi-carré
chi2 = stats.chi2(df=5)

# Student t
t_dist = stats.t(df=10)

# F de Fisher
f_dist = stats.f(dfn=5, dfd=10)

# Log-normale
lognormal = stats.lognorm(s=1, scale=np.exp(0))

# Weibull
weibull = stats.weibull_min(c=1.5)

# Cauchy
cauchy = stats.cauchy(loc=0, scale=1)

# Laplace
laplace = stats.laplace(loc=0, scale=1)

# Pareto
pareto = stats.pareto(b=2.5)

# Rayleigh
rayleigh = stats.rayleigh(scale=1)

# Gumbel
gumbel = stats.gumbel_r(loc=0, scale=1)

# Logistique
logistic = stats.logistic(loc=0, scale=1)

# Triangle
triangular = stats.triang(c=0.5, loc=0, scale=1)

# Multivariate normale
mean = [0, 0]
cov = [[1, 0.5], [0.5, 1]]
mvn = stats.multivariate_normal(mean=mean, cov=cov)
samples_mvn = mvn.rvs(size=100)
pdf_mvn = mvn.pdf([[0, 0], [1, 1]])

# === Distributions discrètes ===

# Binomiale
n, p = 10, 0.5
binomial = stats.binom(n=n, p=p)

# PMF - Probability mass function
pmf = binomial.pmf([0, 1, 2, 3])

# CDF
cdf = binomial.cdf(5)

# Échantillonnage
samples = binomial.rvs(size=1000)

# Poisson
lambda_param = 3
poisson = stats.poisson(mu=lambda_param)

# Géométrique
geometric = stats.geom(p=0.5)

# Binomiale négative
nbinom = stats.nbinom(n=5, p=0.5)

# Hypergéométrique
M, n, N = 20, 7, 12  # Population, succès, tirages
hypergeom = stats.hypergeom(M=M, n=n, N=N)

# Discrète uniforme
discrete_uniform = stats.randint(low=1, high=7)  # Dé

# Multinomiale
n = 10
p = [0.2, 0.3, 0.5]
multinomial = stats.multinomial(n=n, p=p)
samples_multi = multinomial.rvs(size=100)

# === Tests statistiques ===

# Test t de Student (un échantillon)
data = stats.norm.rvs(loc=5, scale=2, size=100)
t_stat, p_value = stats.ttest_1samp(data, popmean=0)
print(f"t={t_stat:.3f}, p={p_value:.3f}")

# Test t (deux échantillons indépendants)
data1 = stats.norm.rvs(loc=5, scale=2, size=100)
data2 = stats.norm.rvs(loc=4, scale=2, size=100)
t_stat, p_value = stats.ttest_ind(data1, data2)

# Sans hypothèse de variance égale (Welch)
t_stat, p_value = stats.ttest_ind(data1, data2, equal_var=False)

# Test t (échantillons appariés)
before = stats.norm.rvs(loc=5, scale=2, size=50)
after = before + stats.norm.rvs(loc=0.5, scale=1, size=50)
t_stat, p_value = stats.ttest_rel(before, after)

# Test de Kolmogorov-Smirnov (normalité)
ks_stat, p_value = stats.kstest(data, 'norm')

# K-S deux échantillons
ks_stat, p_value = stats.ks_2samp(data1, data2)

# Test de Shapiro-Wilk (normalité)
w_stat, p_value = stats.shapiro(data)

# Test de Anderson-Darling (normalité)
result = stats.anderson(data, dist='norm')
print(f"Statistique: {result.statistic}")
print(f"Valeurs critiques: {result.critical_values}")

# Test du Chi-carré (indépendance)
obs = np.array([[10, 10, 20], [20, 20, 20]])
chi2_stat, p_value, dof, expected = stats.chi2_contingency(obs)

# Test exact de Fisher (tables 2x2)
table = [[8, 2], [1, 5]]
odds_ratio, p_value = stats.fisher_exact(table)

# Test de Mann-Whitney U (non-paramétrique, deux échantillons)
u_stat, p_value = stats.mannwhitneyu(data1, data2)

# Test de Wilcoxon (non-paramétrique, appariés)
w_stat, p_value = stats.wilcoxon(before, after)

# Test de Kruskal-Wallis (ANOVA non-paramétrique)
data3 = stats.norm.rvs(loc=6, scale=2, size=100)
h_stat, p_value = stats.kruskal(data1, data2, data3)

# ANOVA à un facteur
f_stat, p_value = stats.f_oneway(data1, data2, data3)

# Test de Levene (homogénéité des variances)
w_stat, p_value = stats.levene(data1, data2, data3)

# Test de Bartlett (homogénéité des variances, normal)
t_stat, p_value = stats.bartlett(data1, data2, data3)

# Test de Friedman (ANOVA non-paramétrique répétée)
measurements = np.array([[5.2, 5.5, 5.8],
                         [4.9, 5.1, 5.3],
                         [5.4, 5.7, 6.0]])
statistic, p_value = stats.friedmanchisquare(*measurements.T)

# Test de corrélation de Pearson
x = np.random.randn(100)
y = x + np.random.randn(100) * 0.5
r, p_value = stats.pearsonr(x, y)

# Test de corrélation de Spearman (non-paramétrique)
rho, p_value = stats.spearmanr(x, y)

# Test de corrélation de Kendall tau
tau, p_value = stats.kendalltau(x, y)

# Test binomial
k = 7  # Succès
n = 10  # Essais
p_null = 0.5
p_value = stats.binom_test(k, n, p_null)

# Test de normalité multivariée (Jarque-Bera)
jb_stat, p_value = stats.jarque_bera(data)

# Test de signe
differences = after - before
statistic, p_value = stats.binomtest(np.sum(differences > 0), len(differences))

# === Statistiques descriptives ===

data = np.random.randn(1000)

# Statistiques de base
mean = np.mean(data)
median = np.median(data)
std = np.std(data, ddof=1)  # ddof=1 pour échantillon
variance = np.var(data, ddof=1)

# Avec scipy.stats
description = stats.describe(data)
# Retourne: (n, (min, max), mean, variance, skewness, kurtosis)

# Moments
mean = stats.tmean(data)  # Moyenne tronquée
variance = stats.tvar(data)
std = stats.tstd(data)
skewness = stats.skew(data)
kurtosis = stats.kurtosis(data)

# Moments géométriques et harmoniques
gmean = stats.gmean(data + 100)  # Moyenne géométrique (positifs)
hmean = stats.hmean(data + 100)  # Moyenne harmonique (positifs)

# Percentiles
q25, q50, q75 = np.percentile(data, [25, 50, 75])
iqr = stats.iqr(data)  # Interquartile range

# Mode
mode_result = stats.mode(data, keepdims=True)
mode_value = mode_result.mode[0]

# Étendue
data_range = np.ptp(data)  # Peak to peak

# Écart médian absolu
mad = stats.median_abs_deviation(data)

# Coefficient de variation
cv = stats.variation(data)

# Entropie
entropy = stats.entropy([0.1, 0.3, 0.6])

# Z-scores
z_scores = stats.zscore(data)

# Moments statistiques d'ordre supérieur
moment_3 = stats.moment(data, moment=3)

# === Ajustement de distribution ===

# Ajuster normale
mu_fit, sigma_fit = stats.norm.fit(data)
print(f"mu={mu_fit:.3f}, sigma={sigma_fit:.3f}")

# Ajuster gamma
params_gamma = stats.gamma.fit(np.abs(data))

# Avec valeurs fixées
params_fixed = stats.norm.fit(data, floc=0)  # Fixer loc=0

# Ajustement automatique (meilleure distribution)
# Tester plusieurs distributions
distributions = [stats.norm, stats.expon, stats.gamma]
best_dist = None
best_aic = np.inf

for dist in distributions:
    params = dist.fit(data)
    # Log-likelihood
    ll = np.sum(dist.logpdf(data, *params))
    # AIC
    k = len(params)
    aic = 2 * k - 2 * ll
    if aic < best_aic:
        best_aic = aic
        best_dist = dist

# Test de goodness-of-fit
D, p_value = stats.kstest(data, 'norm', args=(mu_fit, sigma_fit))

# === Régression linéaire ===

x = np.random.randn(100)
y = 2 * x + 1 + np.random.randn(100) * 0.5

# Régression linéaire simple
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
print(f"y = {slope:.3f}x + {intercept:.3f}")
print(f"R² = {r_value**2:.3f}")

# Prédiction
y_pred = slope * x + intercept

# Résidus
residuals = y - y_pred

# Test de régression de Theil-Sen (robuste)
result = stats.theilslopes(y, x)
slope_ts, intercept_ts = result.slope, result.intercept

# === Bootstrap ===

# Bootstrapping manuel
n_bootstrap = 1000
bootstrap_means = []
for _ in range(n_bootstrap):
    sample = np.random.choice(data, size=len(data), replace=True)
    bootstrap_means.append(np.mean(sample))

# Intervalle de confiance
ci_lower = np.percentile(bootstrap_means, 2.5)
ci_upper = np.percentile(bootstrap_means, 97.5)

# === Intervalles de confiance ===

# IC pour moyenne (normale)
confidence = 0.95
ci = stats.norm.interval(confidence, loc=mean, scale=std/np.sqrt(len(data)))

# IC pour moyenne (t de Student)
ci = stats.t.interval(confidence, len(data)-1, loc=mean, 
                     scale=std/np.sqrt(len(data)))

# IC pour proportion
n_success = 45
n_trials = 100
ci_proportion = stats.binom.interval(0.95, n_trials, n_success/n_trials)

# Wilson score interval (meilleur pour petits échantillons)
from statsmodels.stats.proportion import proportion_confint
ci_wilson = proportion_confint(n_success, n_trials, method='wilson')

# === Transformations ===

# Box-Cox (normalisation)
data_positive = np.abs(data) + 1
transformed, lambda_param = stats.boxcox(data_positive)

# Inverse Box-Cox
original = stats.inv_boxcox(transformed, lambda_param)

# Yeo-Johnson (accepte valeurs négatives)
transformed_yj, lambda_yj = stats.yeojohnson(data)

# === Contingence et association ===

# Table de contingence
obs = np.array([[10, 15, 5], [20, 10, 15]])

# Chi-carré
chi2, p, dof, expected = stats.chi2_contingency(obs)

# Coefficient de contingence de Cramér
n = obs.sum()
min_dim = min(obs.shape) - 1
cramer_v = np.sqrt(chi2 / (n * min_dim))

# === Échantillonnage ===

# Échantillonnage stratifié
from scipy.stats import qmc
sampler = qmc.LatinHypercube(d=2)  # 2 dimensions
sample = sampler.random(n=100)

# Sobol sequence (quasi-random)
sampler_sobol = qmc.Sobol(d=2)
sample_sobol = sampler_sobol.random(n=100)

# Halton sequence
sampler_halton = qmc.Halton(d=2)
sample_halton = sampler_halton.random(n=100)

# === Fonctions de ranking ===

# Rangs
ranks = stats.rankdata(data)

# Rangs moyens pour égalités
ranks_avg = stats.rankdata(data, method='average')

# Rangs avec méthodes différentes
# method: 'average', 'min', 'max', 'dense', 'ordinal'
ranks_min = stats.rankdata(data, method='min')

# === Power analysis ===

# Taille d'effet Cohen's d
group1 = np.random.randn(50)
group2 = np.random.randn(50) + 0.5
pooled_std = np.sqrt((np.var(group1) + np.var(group2)) / 2)
cohens_d = (np.mean(group2) - np.mean(group1)) / pooled_std

# === Distributions multivariées ===

# Dirichlet
alpha = [0.5, 0.5, 0.5]
dirichlet = stats.dirichlet(alpha)
samples_dir = dirichlet.rvs(size=100)

# Wishart (distribution de matrices)
df = 5
scale = np.eye(3)
wishart = stats.wishart(df=df, scale=scale)
sample_wishart = wishart.rvs()

# Multivariate t
mean = [0, 0]
cov = [[1, 0.5], [0.5, 1]]
df = 5
mvt = stats.multivariate_t(loc=mean, shape=cov, df=df)
samples_mvt = mvt.rvs(size=100)


[OK] SCIPY.SIGNAL - TRAITEMENT DU SIGNAL


from scipy import signal
import numpy as np

# === Génération de signaux ===

# Signal carré
t = np.linspace(0, 1, 1000)
square_wave = signal.square(2 * np.pi * 5 * t)  # 5 Hz

# Signal triangulaire
sawtooth = signal.sawtooth(2 * np.pi * 5 * t)
triangle = signal.sawtooth(2 * np.pi * 5 * t, width=0.5)

# Chirp (fréquence variable)
chirp_linear = signal.chirp(t, f0=5, f1=50, t1=1, method='linear')
chirp_quadratic = signal.chirp(t, f0=5, f1=50, t1=1, method='quadratic')
chirp_log = signal.chirp(t, f0=5, f1=50, t1=1, method='logarithmic')
chirp_hyperbolic = signal.chirp(t, f0=5, f1=50, t1=1, method='hyperbolic')

# Sweep cosine
sweep = signal.sweep_poly(t, [5, 50])

# Bruit blanc gaussien
noise = np.random.randn(len(t))

# === Fenêtres ===

# Fenêtre rectangulaire
window_rect = signal.windows.boxcar(51)

# Hann (Hanning)
window_hann = signal.windows.hann(51)

# Hamming
window_hamming = signal.windows.hamming(51)

# Blackman
window_blackman = signal.windows.blackman(51)

# Blackman-Harris
window_blackman_harris = signal.windows.blackmanharris(51)

# Kaiser
beta = 8.6
window_kaiser = signal.windows.kaiser(51, beta)

# Tukey
alpha = 0.5
window_tukey = signal.windows.tukey(51, alpha)

# Bartlett
window_bartlett = signal.windows.bartlett(51)

# Gaussian
std = 7
window_gaussian = signal.windows.gaussian(51, std)

# Fenêtre générale
window_general = signal.windows.general_hamming(51, alpha=0.54)

# Propriétés fenêtre
window = signal.windows.hann(51)
A, B = signal.windows.get_window('hann', 51, fftbins=True)

# === Filtres FIR (Finite Impulse Response) ===

# Filtre passe-bas FIR
numtaps = 51  # Ordre filtre
cutoff = 10  # Fréquence coupure (Hz)
fs = 100  # Fréquence échantillonnage
fir_lowpass = signal.firwin(numtaps, cutoff, fs=fs)

# Filtre passe-haut
fir_highpass = signal.firwin(numtaps, cutoff, fs=fs, pass_zero=False)

# Filtre passe-bande
lowcut = 10
highcut = 20
fir_bandpass = signal.firwin(numtaps, [lowcut, highcut], fs=fs, pass_zero=False)

# Filtre coupe-bande
fir_bandstop = signal.firwin(numtaps, [lowcut, highcut], fs=fs, pass_zero=True)

# Avec fenêtre spécifique
fir_custom = signal.firwin(numtaps, cutoff, fs=fs, window='hamming')

# Méthode firwin2 (réponse arbitraire)
freq = [0, 10, 15, 50]  # Points fréquence
gain = [1, 1, 0, 0]  # Gain désiré
fir_custom2 = signal.firwin2(numtaps, freq, gain, fs=fs)

# Filtre FIR optimal (Parks-McClellan / Remez)
bands = [0, 8, 12, 50]  # Bandes [passband_start, passband_end, stopband_start, stopband_end]
desired = [1, 0]  # Gain désiré dans chaque bande
fir_remez = signal.remez(numtaps, bands, desired, fs=fs)

# Filtre minimum phase
fir_minphase = signal.minimum_phase(fir_lowpass)

# === Filtres IIR (Infinite Impulse Response) ===

# Butterworth passe-bas
order = 4
Wn = 10  # Fréquence coupure normalisée (0-1, où 1 = Nyquist)
b, a = signal.butter(order, Wn, btype='low', fs=fs)

# Passe-haut
b, a = signal.butter(order, Wn, btype='high', fs=fs)

# Passe-bande
b, a = signal.butter(order, [10, 20], btype='band', fs=fs)

# Coupe-bande
b, a = signal.butter(order, [10, 20], btype='bandstop', fs=fs)

# Forme SOS (Second-Order Sections, plus stable)
sos = signal.butter(order, Wn, btype='low', fs=fs, output='sos')

# Chebyshev Type I (ondulation dans bande passante)
ripple_db = 0.5  # Ondulation maximale (dB)
b, a = signal.cheby1(order, ripple_db, Wn, fs=fs)
sos = signal.cheby1(order, ripple_db, Wn, fs=fs, output='sos')

# Chebyshev Type II (ondulation dans bande atténuée)
attenuation_db = 40
b, a = signal.cheby2(order, attenuation_db, Wn, fs=fs)

# Elliptique (Cauer) - ondulation dans les deux bandes
b, a = signal.ellip(order, ripple_db, attenuation_db, Wn, fs=fs)

# Bessel (phase linéaire)
b, a = signal.bessel(order, Wn, fs=fs)

# === Conception filtre IIR (specs analogiques) ===

# Butterworth par specs
wp = 10  # Passband edge
ws = 20  # Stopband edge
gpass = 3  # Passband ripple (dB)
gstop = 40  # Stopband attenuation (dB)
order, Wn = signal.buttord(wp, ws, gpass, gstop, fs=fs)
b, a = signal.butter(order, Wn, fs=fs)

# Chebyshev I par specs
order, Wn = signal.cheb1ord(wp, ws, gpass, gstop, fs=fs)
b, a = signal.cheby1(order, gpass, Wn, fs=fs)

# Chebyshev II par specs
order, Wn = signal.cheb2ord(wp, ws, gpass, gstop, fs=fs)
b, a = signal.cheby2(order, gstop, Wn, fs=fs)

# Elliptique par specs
order, Wn = signal.ellipord(wp, ws, gpass, gstop, fs=fs)
b, a = signal.ellip(order, gpass, gstop, Wn, fs=fs)

# === Filtrage (application) ===

# Signal test
x = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)

# Filtrage avec coefficients b, a
y = signal.lfilter(b, a, x)

# Filtrage avec SOS (recommandé)
y_sos = signal.sosfilt(sos, x)

# Filtrage avant-arrière (zero-phase)
y_filtfilt = signal.filtfilt(b, a, x)

# Avec SOS
y_filtfilt_sos = signal.sosfiltfilt(sos, x)

# Convolution (FIR)
y_fir = signal.convolve(x, fir_lowpass, mode='same')

# Corrélation
y_corr = signal.correlate(x, fir_lowpass, mode='same')

# Filtrage avec fenêtre glissante
y_savgol = signal.savgol_filter(x, window_length=51, polyorder=3)

# Dérivée avec Savitzky-Golay
y_deriv = signal.savgol_filter(x, window_length=51, polyorder=3, deriv=1)

# Médian filter (non-linéaire, élimine outliers)
y_median = signal.medfilt(x, kernel_size=5)

# Wiener filter (estimation optimale)
y_wiener = signal.wiener(x, mysize=5)

# === Analyse fréquentielle ===

# Réponse fréquentielle filtre
w, h = signal.freqz(b, a, worN=2000, fs=fs)
# w: fréquences, h: réponse complexe

# Magnitude et phase
magnitude_db = 20 * np.log10(abs(h))
phase = np.angle(h)

# Avec SOS
w_sos, h_sos = signal.sosfreqz(sos, worN=2000, fs=fs)

# Réponse impulsionnelle
impulse_response = signal.impulse((b, a), N=100)

# Réponse indicielle (step response)
step_response = signal.step((b, a), N=100)

# Délai de groupe (group delay)
w, gd = signal.group_delay((b, a), fs=fs)

# === Transformée de Fourier ===

# FFT
X = np.fft.fft(x)
freq = np.fft.fftfreq(len(x), 1/fs)

# Avec scipy
f, X_scipy = signal.periodogram(x, fs=fs)

# Welch (moyenne de périodogrammes)
f, Pxx = signal.welch(x, fs=fs, nperseg=256)

# Spectrogramme (STFT)
f, t_spec, Sxx = signal.spectrogram(x, fs=fs)

# STFT complète
f, t_stft, Zxx = signal.stft(x, fs=fs, nperseg=256)

# STFT inverse
t_reconstructed, x_reconstructed = signal.istft(Zxx, fs=fs)

# Densité spectrale croisée
f, Pxy = signal.csd(x, y, fs=fs)

# Cohérence
f, Cxy = signal.coherence(x, y, fs=fs)

# === Analyse temps-fréquence ===

# Transformée en ondelettes continue (CWT)
widths = np.arange(1, 31)
cwt_matrix = signal.cwt(x, signal.ricker, widths)

# Ondelettes disponibles: ricker (Mexican hat), morlet, etc

# === Détection de pics ===

# Trouver pics
peaks, properties = signal.find_peaks(x)

# Avec conditions
peaks, _ = signal.find_peaks(x, height=0.5)  # Hauteur minimale
peaks, _ = signal.find_peaks(x, threshold=0.1)  # Différence verticale min
peaks, _ = signal.find_peaks(x, distance=20)  # Distance minimale entre pics
peaks, _ = signal.find_peaks(x, prominence=0.5)  # Proéminence minimale
peaks, _ = signal.find_peaks(x, width=5)  # Largeur minimale
peaks, _ = signal.find_peaks(x, wlen=50)  # Fenêtre évaluation proéminence

# Propriétés des pics
peaks, properties = signal.find_peaks(x, prominence=0.5, width=5)
# properties: 'peak_heights', 'prominences', 'left_bases', 'right_bases',
#            'widths', 'width_heights', 'left_ips', 'right_ips'

# Largeur pic
widths, width_heights, left_ips, right_ips = signal.peak_widths(
    x, peaks, rel_height=0.5
)

# Proéminence pic
prominences, left_bases,# Fichier: python_cheats/cheatsheets/scipy.txt
# Cheatsheet SciPy - Guide Complet des Fonctionnalités Scientifiques


[OK] INSTALLATION & IMPORTS

# Installation
pip install scipy
pip install scipy numpy  # NumPy requis
conda install scipy      # Avec conda

# Import global (déconseillé)
import scipy

# Imports par module (recommandé)
import numpy as np
from scipy import stats, optimize, integrate, interpolate
from scipy import linalg, sparse, signal, ndimage, spatial
from scipy import special, fft, cluster, constants

# Vérifier version
import scipy
print(scipy.__version__)


[OK] SCIPY.CONSTANTS - CONSTANTES PHYSIQUES ET MATHÉMATIQUES


# === Constantes mathématiques ===
from scipy import constants as const

const.pi                    # 3.141592653589793
const.golden                # Nombre d'or: 1.618033988749895
const.golden_ratio          # Même chose

# === Constantes physiques fondamentales ===
const.c                     # Vitesse de la lumière: 299792458.0 m/s
const.speed_of_light        # Même chose
const.h                     # Constante de Planck: 6.62607015e-34 J⋅s
const.Planck                # Même chose
const.hbar                  # ℏ = h/(2π): 1.054571817e-34 J⋅s
const.G                     # Constante gravitationnelle: 6.6743e-11 m³/(kg⋅s²)
const.g                     # Accélération gravité standard: 9.80665 m/s²
const.e                     # Charge élémentaire: 1.602176634e-19 C
const.R                     # Constante des gaz parfaits: 8.314462618 J/(mol⋅K)
const.N_A                   # Nombre d'Avogadro: 6.02214076e23 mol⁻¹
const.k                     # Constante de Boltzmann: 1.380649e-23 J/K
const.sigma                 # Constante Stefan-Boltzmann: 5.670374419e-8 W/(m²⋅K⁴)

# === Masse des particules ===
const.m_e                   # Masse électron: 9.1093837015e-31 kg
const.m_p                   # Masse proton: 1.67262192369e-27 kg
const.m_n                   # Masse neutron: 1.67492749804e-27 kg
const.m_u                   # Unité de masse atomique: 1.66053906660e-27 kg

# === Constantes électromagnétiques ===
const.mu_0                  # Perméabilité vide: 1.25663706212e-6 N/A²
const.epsilon_0             # Permittivité vide: 8.8541878128e-12 F/m
const.e                     # Charge élémentaire

# === Conversions d'unités ===

# Distance
const.inch                  # Pouce vers mètre: 0.0254
const.foot                  # Pied vers mètre: 0.3048
const.mile                  # Mile vers mètre: 1609.344
const.nautical_mile         # Mile nautique: 1852.0
const.angstrom              # Angström: 1e-10 m
const.light_year            # Année-lumière: 9.4607304725808e15 m
const.parsec                # Parsec: 3.0856775814913673e16 m
const.au                    # Unité astronomique: 1.495978707e11 m

# Temps
const.minute                # 60.0 secondes
const.hour                  # 3600.0 secondes
const.day                   # 86400.0 secondes
const.week                  # 604800.0 secondes
const.year                  # 31536000.0 secondes (365 jours)
const.Julian_year           # 31557600.0 secondes (365.25 jours)

# Masse
const.gram                  # 0.001 kg
const.metric_ton            # 1000.0 kg
const.grain                 # 6.479891e-05 kg
const.pound                 # Livre: 0.45359237 kg
const.ounce                 # Once: 0.028349523125 kg
const.stone                 # 6.35029318 kg
const.ton_TNT               # 4.184e9 J

# Volume
const.liter                 # Litre: 0.001 m³
const.gallon                # Gallon US: 0.003785411784 m³
const.gallon_US             # Même chose
const.gallon_imp            # Gallon impérial: 0.00454609 m³
const.fluid_ounce           # 2.9573529562e-05 m³
const.barrel                # Baril: 0.158987294928 m³

# Pression
const.atm                   # Atmosphère: 101325.0 Pa
const.atmosphere            # Même chose
const.bar                   # 100000.0 Pa
const.torr                  # 133.32236842105263 Pa
const.mmHg                  # Millimètre mercure: 133.322387415 Pa
const.psi                   # Livre par pouce carré: 6894.757293168361 Pa

# Énergie
const.eV                    # Électron-volt: 1.602176634e-19 J
const.electron_volt         # Même chose
const.calorie               # Calorie: 4.184 J
const.calorie_th            # Calorie thermochimique: 4.184 J
const.calorie_IT            # Calorie internationale: 4.1868 J
const.erg                   # 1e-7 J
const.Btu                   # British thermal unit: 1055.05585262 J
const.Btu_IT                # BTU international: 1055.05585262 J
const.ton_TNT               # Tonne de TNT: 4.184e9 J

# Puissance
const.hp                    # Cheval-vapeur: 745.6998715822702 W
const.horsepower            # Même chose

# Température (conversions)
const.zero_Celsius          # 0°C en Kelvin: 273.15 K
const.degree_Fahrenheit     # 0.5555555555555556 K

# Angle
const.degree                # Degré vers radian: 0.017453292519943295
const.arcmin                # Minute d'arc: 0.0002908882086657216 rad
const.arcsec                # Seconde d'arc: 4.84813681109536e-06 rad

# === Préfixes SI ===
const.yotta                 # 1e24
const.zetta                 # 1e21
const.exa                   # 1e18
const.peta                  # 1e15
const.tera                  # 1e12
const.giga                  # 1e9
const.mega                  # 1e6
const.kilo                  # 1000.0
const.hecto                 # 100.0
const.deka                  # 10.0
const.deci                  # 0.1
const.centi                 # 0.01
const.milli                 # 0.001
const.micro                 # 1e-6
const.nano                  # 1e-9
const.pico                  # 1e-12
const.femto                 # 1e-15
const.atto                  # 1e-18
const.zepto                 # 1e-21

# === Préfixes binaires ===
const.kibi                  # 1024
const.mebi                  # 1048576
const.gibi                  # 1073741824
const.tebi                  # 1099511627776
const.pebi                  # 1125899906842624
const.exbi                  # 1152921504606846976

# === Recherche de constantes ===
const.find('light')         # Cherche toutes les constantes contenant 'light'
const.physical_constants    # Dict de toutes les constantes physiques
const.unit('light year')    # Détails sur une constante

# === Conversions de température ===
from scipy.constants import convert_temperature
convert_temperature(100, 'Celsius', 'Fahrenheit')      # 212.0
convert_temperature(32, 'Fahrenheit', 'Celsius')       # 0.0
convert_temperature(273.15, 'Kelvin', 'Celsius')       # 0.0
convert_temperature(0, 'Celsius', 'Kelvin')            # 273.15

# === Valeur, unité et incertitude ===
val, unit, uncertainty = const.physical_constants['proton mass']
# val: 1.67262192369e-27
# unit: 'kg'
# uncertainty: 5.1e-37


[OK] SCIPY.SPECIAL - FONCTIONS SPÉCIALES MATHÉMATIQUES


from scipy import special
import numpy as np

# === Fonctions de Bessel ===

# Bessel première espèce
special.jv(0, 2.4)          # J₀(2.4)
special.jv(1, [1, 2, 3])    # J₁ pour plusieurs valeurs
special.j0(2.4)             # J₀ optimisé
special.j1(2.4)             # J₁ optimisé

# Bessel deuxième espèce (Neumann)
special.yv(0, 2.4)          # Y₀(2.4)
special.y0(2.4)             # Y₀ optimisé
special.y1(2.4)             # Y₁ optimisé

# Bessel modifiées
special.iv(0, 2.4)          # I₀(2.4) - première espèce modifiée
special.kv(0, 2.4)          # K₀(2.4) - deuxième espèce modifiée

# Bessel sphériques
special.spherical_jn(2, 1.0)    # j₂(1.0)
special.spherical_yn(2, 1.0)    # y₂(1.0)

# Zéros des fonctions de Bessel
special.jn_zeros(0, 5)      # 5 premiers zéros de J₀
special.jnp_zeros(0, 5)     # 5 premiers zéros de J₀'

# === Fonctions gamma et associées ===

special.gamma(5)            # Γ(5) = 4! = 24
special.gammaln(100)        # ln(Γ(100)) - évite overflow
special.loggamma(100)       # Même chose
special.gammainc(2, 1)      # Gamma incomplète normalisée γ(a,x)/Γ(a)
special.gammaincc(2, 1)     # Gamma incomplète complémentaire
special.digamma(5)          # ψ(x) = Γ'(x)/Γ(x) - fonction digamma
special.polygamma(1, 5)     # Polygamma d'ordre 1
special.rgamma(5)           # 1/Γ(5)

# Fonction beta
special.beta(2, 3)          # B(2,3) = Γ(2)Γ(3)/Γ(5)
special.betaln(2, 3)        # ln(B(2,3))
special.betainc(2, 3, 0.5)  # Beta incomplète normalisée

# === Fonctions d'erreur ===

special.erf(1.0)            # Fonction d'erreur erf(x)
special.erfc(1.0)           # Fonction d'erreur complémentaire 1-erf(x)
special.erfcx(1.0)          # exp(x²)·erfc(x) - évite underflow
special.erfi(1.0)           # Fonction d'erreur imaginaire
special.erfinv(0.5)         # Inverse de erf
special.erfcinv(0.5)        # Inverse de erfc

# Intégrales de Fresnel
special.fresnel(1.0)        # Retourne (S(x), C(x))

# === Fonctions elliptiques ===

# Intégrales elliptiques complètes
special.ellipk(0.5)         # K(m) - première espèce
special.ellipe(0.5)         # E(m) - deuxième espèce
special.ellipkm1(0.5)       # K(1-m)

# Intégrales elliptiques incomplètes
special.ellipkinc(np.pi/4, 0.5)  # F(φ,m) - première espèce
special.ellipeinc(np.pi/4, 0.5)  # E(φ,m) - deuxième espèce

# Fonctions elliptiques de Jacobi
special.ellipj(0.5, 0.3)    # Retourne (sn, cn, dn, ph)

# === Polynômes orthogonaux ===

# Legendre
special.legendre(3)         # Retourne polynôme P₃(x)
x = np.linspace(-1, 1, 100)
P3 = special.eval_legendre(3, x)  # Évalue P₃(x)
special.lpmv(2, 3, 0.5)     # Legendre associé P₃²(0.5)

# Chebyshev
special.chebyt(3)           # Polynôme Chebyshev 1ère espèce T₃(x)
special.chebyu(3)           # Polynôme Chebyshev 2ème espèce U₃(x)
special.eval_chebyt(3, x)   # Évalue T₃(x)

# Hermite
special.hermite(3)          # Polynôme Hermite H₃(x) (physicien)
special.hermitenorm(3)      # Hermite normalisé (probabiliste)
special.eval_hermite(3, x)  # Évalue H₃(x)

# Laguerre
special.laguerre(3)         # Polynôme Laguerre L₃(x)
special.genlaguerre(3, 0.5) # Laguerre généralisé L₃⁰·⁵(x)
special.eval_laguerre(3, x) # Évalue L₃(x)

# Gegenbauer
special.gegenbauer(3, 0.5)  # Polynôme Gegenbauer C₃⁰·⁵(x)

# Jacobi
special.jacobi(3, 0.5, 0.5) # Polynôme Jacobi P₃⁽⁰·⁵'⁰·⁵⁾(x)

# === Fonctions hypergeométriques ===

special.hyp2f1(1, 2, 3, 0.5)    # ₂F₁(a,b;c;z) - Gauss
special.hyp1f1(1, 2, 0.5)       # ₁F₁(a;b;z) - Kummer
special.hyp0f1(2, 0.5)          # ₀F₁(b;z)
special.hyperu(1, 2, 0.5)       # U(a,b,x) - Tricomi

# === Fonctions de Airy ===

special.airy(1.0)           # Retourne (Ai, Ai', Bi, Bi')
special.ai_zeros(5)         # 5 premiers zéros de Ai(x)
special.bi_zeros(5)         # 5 premiers zéros de Bi(x)

# === Fonctions combinatoires ===

special.comb(10, 3)         # Combinaisons C(10,3) = 120
special.comb(10, 3, exact=True)  # Résultat exact (entier)
special.perm(10, 3)         # Permutations P(10,3) = 720
special.factorial(5)        # 5! = 120
special.factorial2(5)       # Double factorielle 5!! = 5·3·1 = 15
special.factorialk(5, 3)    # Factorielle généralisée

# === Fonctions exponentielles intégrales ===

special.expi(1.0)           # Ei(x) = -∫_{-x}^∞ e^{-t}/t dt
special.exp1(1.0)           # E₁(x) = ∫_x^∞ e^{-t}/t dt
special.expn(2, 1.0)        # Eₙ(x) - intégrale exponentielle généralisée

# Intégrales sinus et cosinus
special.sici(1.0)           # Retourne (Si(x), Ci(x))
special.shichi(1.0)         # Retourne (Shi(x), Chi(x))

# === Fonctions zeta et associées ===

special.zeta(2)             # ζ(2) = π²/6
special.zetac(2)            # ζ(x) - 1

# === Fonctions de Lambert W ===

special.lambertw(1)         # W(1) - branche principale
special.lambertw(1, k=-1)   # Branche -1

# === Fonctions logistiques ===

special.expit(0)            # 1/(1+exp(-x)) - sigmoïde: 0.5
special.logit(0.5)          # log(p/(1-p)) - inverse sigmoïde: 0.0
special.log_expit(0)        # log(expit(x)) - stable numériquement

# === Fonctions de Mathieu ===

special.mathieu_a(2, 5)     # Valeur caractéristique a
special.mathieu_b(2, 5)     # Valeur caractéristique b
special.mathieu_cem(2, 5, 0.5)  # Fonction Mathieu ce_m
special.mathieu_sem(2, 5, 0.5)  # Fonction Mathieu se_m

# === Fonctions de Kelvin ===

special.kelvin(2)           # Retourne (Be, Ke, Bep, Kep)
special.bei(1.0)            # bei(x)
special.ber(1.0)            # ber(x)
special.kei(1.0)            # kei(x)
special.ker(1.0)            # ker(x)

# === Fonctions de Struve ===

special.struve(0, 1.0)      # H₀(x)
special.modstruve(0, 1.0)   # L₀(x) - Struve modifiée

# === Fonctions de partition sphérique ===

special.sph_harm(1, 2, np.pi/4, np.pi/3)  # Y₂¹(θ,φ)

# === Fonctions diverses ===

special.sinc(np.pi)         # sin(πx)/(πx): retourne 0
special.sindg(90)           # sin en degrés: 1.0
special.cosdg(180)          # cos en degrés: -1.0
special.tandg(45)           # tan en degrés: 1.0

# Intégrale de Dawson
special.dawsn(1.0)          # F(x) = exp(-x²)∫₀ˣ exp(t²) dt

# Voigt profile
special.voigt_profile(1.0, 2.0, 0.5)  # Convolution Gauss-Lorentz

# Spence (dilogarithme)
special.spence(0.5)         # Li₂(x)

# Fonction de Hurwitz zeta
special.hurwitz(2, 1)       # ζ(s,a)

# === Conversions et utilitaires ===

# Décibels
special.db2pow(20)          # dB vers puissance: 100
special.pow2db(100)         # Puissance vers dB: 20

# Softmax
x = np.array([1, 2, 3])
special.softmax(x)          # exp(x)/sum(exp(x))
special.log_softmax(x)      # log(softmax(x)) - stable

# xlogy - x*log(y) stable pour x=0
special.xlogy(0, 0)         # Retourne 0 (pas NaN)
special.xlog1py(2, 0.1)     # x*log(1+y)

# === Racines et zéros ===

# Racines des polynômes orthogonaux
roots, weights = special.roots_legendre(5)  # Gauss-Legendre
roots, weights = special.roots_chebyt(5)    # Gauss-Chebyshev
roots, weights = special.roots_hermite(5)   # Gauss-Hermite
roots, weights = special.roots_laguerre(5)  # Gauss-Laguerre
roots, weights = special.roots_jacobi(5, 0.5, 0.5)  # Gauss-Jacobi

# Utilisées pour intégration numérique (quadrature de Gauss)


[OK] SCIPY.OPTIMIZE - OPTIMISATION ET RECHERCHE DE RACINES


from scipy import optimize
import numpy as np

# === Minimisation scalaire (1D) ===

# Fonction à minimiser
def f(x):
    return (x - 2) ** 2 + 3

# Méthode de Brent (bracketing)
result = optimize.minimize_scalar(f, method='brent')
print(result.x, result.fun)  # 2.0, 3.0

# Avec intervalle borné
result = optimize.minimize_scalar(f, bounds=(0, 5), method='bounded')

# Golden section search
result = optimize.minimize_scalar(f, bracket=(0, 5), method='golden')

# === Minimisation non contrainte (multidimensionnelle) ===

# Fonction de Rosenbrock (exemple classique)
def rosen(x):
    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)

# Point initial
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])

# Nelder-Mead (simplex) - sans gradient
result = optimize.minimize(rosen, x0, method='nelder-mead')
result = optimize.minimize(rosen, x0, method='Nelder-Mead',
                          options={'xatol': 1e-8, 'disp': True})

# Powell - sans gradient
result = optimize.minimize(rosen, x0, method='powell')

# BFGS - quasi-Newton avec gradient
result = optimize.minimize(rosen, x0, method='BFGS')

# Avec gradient fourni
def rosen_grad(x):
    xm = x[1:-1]
    xm_m1 = x[:-2]
    xm_p1 = x[2:]
    der = np.zeros_like(x)
    der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
    der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
    der[-1] = 200*(x[-1]-x[-2]**2)
    return der

result = optimize.minimize(rosen, x0, method='BFGS', jac=rosen_grad)

# L-BFGS-B - mémoire limitée, supporte bornes
bounds = [(0, None)] * len(x0)  # x >= 0
result = optimize.minimize(rosen, x0, method='L-BFGS-B', bounds=bounds)

# Conjugate Gradient
result = optimize.minimize(rosen, x0, method='CG', jac=rosen_grad)

# Newton-CG - nécessite hessienne
def rosen_hess(x):
    # Hessienne de Rosenbrock
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H

result = optimize.minimize(rosen, x0, method='Newton-CG',
                          jac=rosen_grad, hess=rosen_hess)

# Trust Region methods
result = optimize.minimize(rosen, x0, method='trust-ncg',
                          jac=rosen_grad, hess=rosen_hess)
result = optimize.minimize(rosen, x0, method='trust-exact',
                          jac=rosen_grad, hess=rosen_hess)
result = optimize.minimize(rosen, x0, method='trust-krylov',
                          jac=rosen_grad, hess=rosen_hess)

# dogleg
result = optimize.minimize(rosen, x0, method='dogleg',
                          jac=rosen_grad, hess=rosen_hess)

# === Minimisation avec contraintes ===

# Contraintes linéaires et non-linéaires
def objective(x):
    return x[0]**2 + x[1]**2

# Contrainte: x[0] + x[1] >= 1
def constraint1(x):
    return x[0] + x[1] - 1

# Contraintes sous forme de dict
cons = ({'type': 'ineq', 'fun': constraint1})

# SLSQP - Sequential Least Squares Programming
result = optimize.minimize(objective, [0, 0], method='SLSQP',
                          constraints=cons)

# Avec bornes et contraintes multiples
bounds = [(-10, 10), (-10, 10)]
cons = [
    {'type': 'ineq', 'fun': lambda x: x[0] + x[1] - 1},
    {'type': 'eq', 'fun': lambda x: x[0]**2 + x[1]**2 - 4}
]
result = optimize.minimize(objective, [1, 1], method='SLSQP',
                          bounds=bounds, constraints=cons)

# Trust-region constrained
result = optimize.minimize(objective, [1, 1], method='trust-constr',
                          bounds=bounds, constraints=cons)

# COBYLA - Constrained Optimization BY Linear Approximation
result = optimize.minimize(objective, [1, 1], method='COBYLA',
                          constraints=cons)

# === Programmation linéaire ===

# Minimiser c^T x sous contraintes A_ub x <= b_ub et A_eq x = b_eq
c = [-1, 4]  # Coefficients fonction objectif
A_ub = [[3, 1], [1, 2]]  # Contraintes inégalité
b_ub = [9, 8]
A_eq = [[1, 1]]  # Contraintes égalité
b_eq = [5]
bounds = [(0, None), (0, None)]  # x >= 0, y >= 0

result = optimize.linprog(c, A_ub=A_ub, b_ub=b_ub,
                         A_eq=A_eq, b_eq=b_eq,
                         bounds=bounds, method='highs')
print(result.x)  # Solution optimale
print(result.fun)  # Valeur optimale

# Méthodes disponibles
# method='highs' (défaut, recommandé)
# method='highs-ds' (dual simplex)
# method='highs-ipm' (interior-point)
# method='interior-point' (legacy)
# method='revised simplex' (legacy)
# method='simplex' (legacy)

# === Recherche de racines (scalaire) ===

def func(x):
    return x**2 - 2

# Méthode de Brent (bracketing)
root = optimize.brentq(func, 0, 3)  # √2 ≈ 1.414
root = optimize.brenth(func, 0, 3)  # Hyperbolic variant

# Bisection
root = optimize.bisect(func, 0, 3, xtol=1e-6)

# Ridder
root = optimize.ridder(func, 0, 3)

# Toms748
root = optimize.toms748(func, 0, 3)

# Newton-Raphson (nécessite dérivée)
def func_prime(x):
    return 2*x

root = optimize.newton(func, 1.0, fprime=func_prime)

# Secant method (approxime dérivée)
root = optimize.newton(func, 1.0)

# Fixed point iteration
def g(x):
    return np.sqrt(2)  # x = g(x)

root = optimize.fixed_point(g, 1.0)

# === Recherche de racines (multidimensionnelle) ===

# Système d'équations
def equations(x):
    return [x[0]**2 + x[1]**2 - 4,
            x[0] - x[1]**2]

# fsolve - hybride Powell
solution = optimize.fsolve(equations, [1, 1])

# root - interface unifiée
result = optimize.root(equations, [1, 1], method='hybr')
result = optimize.root(equations, [1, 1], method='lm')  # Levenberg-Marquardt
result = optimize.root(equations, [1, 1], method='broyden1')
result = optimize.root(equations, [1, 1], method='broyden2')
result = optimize.root(equations, [1, 1], method='anderson')
result = optimize.root(equations, [1, 1], method='krylov')

# Avec jacobienne
def jacobian(x):
    return [[2*x[0], 2*x[1]],
            [1, -2*x[1]]]

result = optimize.root(equations, [1, 1], jac=jacobian, method='hybr')

# === Moindres carrés (curve fitting) ===

# Données expérimentales
xdata = np.array([0, 1, 2, 3, 4])
ydata = np.array([1.1, 2.9, 5.2, 7.1, 8.9])

# Modèle: y = a*x + b
def model(x, a, b):
    return a * x + b

# Ajustement de courbe
params, covariance = optimize.curve_fit(model, xdata, ydata)
print(f"a={params[0]:.2f}, b={params[1]:.2f}")

# Avec estimation initiale
p0 = [1, 0]
params, cov = optimize.curve_fit(model, xdata, ydata, p0=p0)

# Avec bornes sur paramètres
bounds = ([0, -np.inf], [10, np.inf])  # 0 <= a <= 10, b sans limite
params, cov = optimize.curve_fit(model, xdata, ydata, bounds=bounds)

# Avec poids (incertitudes)
sigma = np.array([0.1, 0.2, 0.1, 0.15, 0.2])
params, cov = optimize.curve_fit(model, xdata, ydata, sigma=sigma)

# Erreurs sur paramètres
perr = np.sqrt(np.diag(cov))

# Moindres carrés linéaires
A = np.vstack([xdata, np.ones(len(xdata))]).T
result = optimize.lsq_linear(A, ydata)

# Moindres carrés non-linéaires
def residuals(params, x, y):
    return model(x, *params) - y

result = optimize.least_squares(residuals, [1, 0], args=(xdata, ydata))
result = optimize.least_squares(residuals, [1, 0], args=(xdata, ydata),
                                method='trf')  # Trust Region Reflective
result = optimize.least_squares(residuals, [1, 0], args=(xdata, ydata),
                                method='dogbox')  # dogleg
result = optimize.least_squares(residuals, [1, 0], args=(xdata, ydata),
                                method='lm')  # Levenberg-Marquardt

# Avec bornes
bounds = ([0, -np.inf], [10, np.inf])
result = optimize.least_squares(residuals, [1, 0], bounds=bounds,
                                args=(xdata, ydata))

# Loss functions robustes (outliers)
result = optimize.least_squares(residuals, [1, 0], loss='soft_l1',
                                args=(xdata, ydata))
# loss='linear' (défaut), 'soft_l1', 'huber', 'cauchy', 'arctan'

# === Optimisation globale ===

# Fonction avec multiples minima locaux
def ackley(x):
    return (-20 * np.exp(-0.2 * np.sqrt(0.5 * (x[0]**2 + x[1]**2))) -
            np.exp(0.5 * (np.cos(2*np.pi*x[0]) + np.cos(2*np.pi*x[1]))) +
            np.e + 20)

# Differential Evolution
bounds = [(-5, 5), (-5, 5)]
result = optimize.differential_evolution(ackley, bounds)
result = optimize.differential_evolution(ackley, bounds, 
                                        strategy='best1bin',
                                        maxiter=1000,
                                        popsize=15,
                                        tol=0.01,
                                        mutation=(0.5, 1),
                                        recombination=0.7)

# Basin-hopping (perturbations aléatoires + minimisation locale)
x0 = [1, 1]
minimizer_kwargs = {"method": "BFGS"}
result = optimize.basinhopping(ackley, x0, minimizer_kwargs=minimizer_kwargs,
                               niter=100)

# Simulated Annealing
result = optimize.dual_annealing(ackley, bounds)

# SHGO - Simplicial Homology Global Optimization
result = optimize.shgo(ackley, bounds)

# === Recherche par force brute ===

result = optimize.brute(ackley, bounds, Ns=50)  # Grille 50x50

# Avec finition par minimisation locale
result = optimize.brute(ackley, bounds, Ns=20, finish=optimize.fmin)

# === Optimisation linéaire assignation ===

# Matrice de coûts
cost = np.array([[4, 1, 3],
                 [2, 0, 5],
                 [3, 2, 2]])

# Hungarian algorithm
row_ind, col_ind = optimize.linear_sum_assignment(cost)
print(cost[row_ind, col_ind].sum())  # Coût total minimal

# === Optimisation quadratique ===

# Minimiser 0.5*x^T*P*x + q^T*x sous contraintes
from scipy.optimize import minimize

P = np.array([[1, 0], [0, 1]])
q = np.array([1, 1])

def quadratic(x):
    return 0.5 * x @ P @ x + q @ x

result = optimize.minimize(quadratic, [0, 0], method='trust-constr')

# === Utilitaires ===

# Approximation numérique du gradient
def f(x):
    return sum(x**2)

x = np.array([1, 2, 3])
grad = optimize.approx_fprime(x, f, epsilon=1e-8)

# Check gradient
def grad_f(x):
    return 2 * x

error = optimize.check_grad(f, grad_f, x)

# Hessienne approximée
hess = optimize.approx_hess(x, f)

# === Options avancées ===

# Callback pour monitoring
def callback(xk):
    print(f"Iteration: {xk}")

result = optimize.minimize(rosen, x0, callback=callback)

# Contraintes avec jacobien
def con_f(x):
    return x[0]**2 + x[1]**2 - 1

def con_jac(x):
    return [2*x[0], 2*x[1]]

cons = {'type': 'eq', 'fun': con_f, 'jac': con_jac}


[OK] SCIPY.INTEGRATE - INTÉGRATION NUMÉRIQUE


from scipy import integrate
import numpy as np

# === Intégration quadrature simple ===

# Fonction à intégrer
def f(x):
    return np.exp(-x**2)

# Intégrale générale
result, error = integrate.quad(f, 0, 1)  # ∫₀¹ e^(-x²) dx
print(f"Résultat: {result:.6f}, Erreur: {error:.2e}")

# Intégrale infinie
result, error = integrate.quad(f, 0, np.inf)  # ∫₀^∞ e^(-x²) dx

# Intégrale double infinie
result, error = integrate.quad(f, -np.inf, np.inf)

# Avec points singuliers
def g(x):
    return 1/np.sqrt(abs(x))

result, error = integrate.quad(g, -1, 1, points=[0])  # Singularité en 0

# Avec poids (weight functions)
result, error = integrate.quad(lambda x: x**2, 0, 1, weight='alg',
                               wvar=(0, 0))  # Poids algébrique

# Poids disponibles:
# 'cos': cos(w*x)
# 'sin': sin(w*x)
# 'alg': (x-a)^α * (b-x)^β
# 'alg-loga': (x-a)^α * (b-x)^β * log(x-a)
# 'alg-logb': (x-a)^α * (b-x)^β * log(b-x)
# 'alg-log': (x-a)^α * (b-x)^β * log(x-a) * log(b-x)
# 'cauchy': 1/(x-c)

# Oscillatoire (cosinus)
result, error = integrate.quad(lambda x: np.cos(10*x), 0, np.pi)

# Avec tolérance
result, error = integrate.quad(f, 0, 1, epsabs=1e-10, epsrel=1e-10)

# Limite de subdivisions
result, error = integrate.quad(f, 0, 1, limit=100)

# === Intégration double ===

# Limites constantes
def integrand(y, x):
    return x * y**2

result, error = integrate.dblquad(integrand, 0, 2, 0, 1)
# ∫₀² ∫₀¹ xy² dy dx

# Limites variables
def y_lower(x):
    return 0

def y_upper(x):
    return x

result, error = integrate.dblquad(integrand, 0, 1, y_lower, y_upper)
# ∫₀¹ ∫₀ˣ xy² dy dx

# === Intégration triple ===

def integrand3d(z, y, x):
    return x * y * z

def z_lower(x, y):
    return 0

def z_upper(x, y):
    return x + y

result, error = integrate.tplquad(integrand3d,
                                  0, 1,        # x: 0 à 1
                                  lambda x: 0, lambda x: 1,  # y: 0 à 1
                                  z_lower, z_upper)  # z: 0 à x+y

# === Intégration n-dimensionnelle ===

# nquad - généralisation arbitraire
def integrand_nd(*args):
    return sum(args)

# 3D avec limites constantes
ranges = [[0, 1], [0, 2], [0, 3]]
result, error = integrate.nquad(integrand_nd, ranges)

# Avec limites fonctions
def y_range(x):
    return [0, x]

def z_range(x, y):
    return [0, x + y]

result, error = integrate.nquad(integrand3d,
                                [z_range, y_range, [0, 1]])

# === Intégration sur échantillon fixe ===

# Trapèzes
x = np.linspace(0, 2, 100)
y = x**2
result = integrate.trapezoid(y, x)  # ou integrate.trapz (deprecated)

# Simpson
result = integrate.simpson(y, x)  # ou integrate.simps (deprecated)

# Romberg (nécessite 2^k+1 points)
x = np.linspace(0, 2, 33)
y = x**2
result = integrate.romb(y, dx=(x[1]-x[0]))

# Cumulative trapezoid (intégrale cumulative)
result = integrate.cumulative_trapezoid(y, x, initial=0)

# === Quadrature de Gauss ===

# Intégration avec points et poids optimaux
def f(x):
    return x**2

# Gauss-Legendre
deg = 5
x, w = np.polynomial.legendre.leggauss(deg)
# Transformation [-1,1] -> [a,b]
a, b = 0, 2
x_transformed = 0.5 * (b - a) * x + 0.5 * (b + a)
result = 0.5 * (b - a) * sum(w * f(x_transformed))

# Avec scipy.special
from scipy.special import roots_legendre
x, w = roots_legendre(deg)
# Même transformation

# Quadrature fixe
result = integrate.fixed_quad(lambda x: x**2, 0, 2, n=5)  # n points

# Quadrature gaussienne adaptative
result, error = integrate.quadrature(lambda x: x**2, 0, 2)

# === ODE - Équations différentielles ordinaires ===

# solve_ivp - Solveur moderne (Python 3.6+)

# dy/dt = -2y, y(0) = 1
def dydt(t, y):
    return -2 * y

t_span = (0, 5)  # Intervalle temps
y0 = [1]  # Condition initiale

# Méthode RK45 (Runge-Kutta 4-5, défaut)
sol = integrate.solve_ivp(dydt, t_span, y0)
print(sol.t)  # Points temps
print(sol.y)  # Solutions

# Évaluation à des temps spécifiques
t_eval = np.linspace(0, 5, 50)
sol = integrate.solve_ivp(dydt, t_span, y0, t_eval=t_eval)

# Méthodes disponibles
sol = integrate.solve_ivp(dydt, t_span, y0, method='RK45')  # Défaut
sol = integrate.solve_ivp(dydt, t_span, y0, method='RK23')  # Ordre 2-3
sol = integrate.solve_ivp(dydt, t_span, y0, method='DOP853')  # Ordre 8
sol = integrate.solve_ivp(dydt, t_span, y0, method='Radau')  # Implicit, stiff
sol = integrate.solve_ivp(dydt, t_span, y0, method='BDF')  # Backward diff, stiff
sol = integrate.solve_ivp(dydt, t_span, y0, method='LSODA')  # Auto-switch

# Système d'EDO: oscillateur harmonique
# d²x/dt² = -ω²x  =>  dx/dt = v, dv/dt = -ω²x
def harmonic(t, y):
    x, v = y
    omega = 2 * np.pi
    return [v, -omega**2 * x]

y0 = [1, 0]  # x(0)=1, v(0)=0
sol = integrate.solve_ivp(harmonic, (0, 10), y0, t_eval=np.linspace(0, 10, 200))

# Avec paramètres
def dydt_param(t, y, k):
    return -k * y

sol = integrate.solve_ivp(dydt_param, t_span, y0, args=(2,))

# Avec jacobien (pour méthodes stiff)
def jacobian(t, y):
    return [[-2]]

sol = integrate.solve_ivp(dydt, t_span, y0, method='Radau', jac=jacobian)

# Événements (détection zéros)
def event(t, y):
    return y[0] - 0.5  # Détecte y = 0.5

event.terminal = True  # Arrête intégration
event.direction = -1   # Seulement décroissant (-1), croissant (1), tous (0)

sol = integrate.solve_ivp(harmonic, (0, 10), y0, events=event)
print(sol.t_events)  # Temps où événement détecté

# Dense output (interpolation)
sol = integrate.solve_ivp(dydt, t_span, y0, dense_output=True)
t_dense = np.linspace(0, 5, 100)
y_dense = sol.sol(t_dense)  # Interpolation continue

# Tolérance
sol = integrate.solve_ivp(dydt, t_span, y0, rtol=1e-6, atol=1e-9)

# Pas maximal
sol = integrate.solve_ivp(dydt, t_span, y0, max_step=0.1)

# === odeint - Interface legacy (plus simple) ===

def dydt_old(y, t):  # Ordre arguments inversé!
    return -2 * y

t = np.linspace(0, 5, 50)
y0 = 1
sol = integrate.odeint(dydt_old, y0, t)

# Système d'EDO
def system(y, t):
    x, v = y
    return [v, -(2*np.pi)**2 * x]

y0 = [1, 0]
sol = integrate.odeint(system, y0, t)

# Avec jacobien
def jac(y, t):
    return [[0, 1], [-(2*np.pi)**2, 0]]

sol = integrate.odeint(system, y0, t, Dfun=jac)

# === Boundary Value Problems (BVP) ===

# solve_bvp - Problèmes aux limites

# Exemple: y'' + y = 0 avec y(0)=0, y(π)=0
def ode(x, y):
    return np.vstack((y[1], -y[0]))

def bc(ya, yb):
    return np.array([ya[0], yb[0]])  # y(0)=0, y(π)=0

x = np.linspace(0, np.pi, 5)
y_guess = np.zeros((2, x.size))  # [y, y']

sol = integrate.solve_bvp(ode, bc, x, y_guess)
print(sol.success)

# Évaluation sur grille fine
x_plot = np.linspace(0, np.pi, 100)
y_plot = sol.sol(x_plot)

# Avec paramètres inconnus
def ode_param(x, y, p):
    k = p[0]  # Paramètre à déterminer
    return np.vstack((y[1], -k*y[0]))

def bc_param(ya, yb, p):
    return np.array([ya[0], yb[0] - 1, ya[1]])  # 3 conditions pour 2+1 inconnues

p_guess = [1]
sol = integrate.solve_bvp(ode_param, bc_param, x, y_guess, p=p_guess)
print(sol.p)  # Paramètre trouvé

# === Intégration complexe ===

# Quadrature complexe
def f_complex(x):
    return np.exp(1j * x)

result = integrate.quad(f_complex, 0, 2*np.pi)
# result est tuple (real_part, imag_part), error

# === Intégration Monte Carlo (avec scipy.stats) ===

from scipy import stats

# Intégration par échantillonnage
def f(x):
    return x**2

n_samples = 100000
x = np.random.uniform(0, 2, n_samples)
result = (2 - 0) * np.mean(f(x))  # Monte Carlo estimate

# === Intégrales à plusieurs dimensions (quasi-MC) ===

# qmc_quad pour haute dimension (nécessite scipy >= 1.8)
# from scipy.integrate import qmc_quad
# def f_nd(x):
#     return np.prod(x, axis=0)
# result = qmc_quad(f_nd, [(0, 1)] * 5, n_points=10000)  # 5D


[OK] SCIPY.INTERPOLATE - INTERPOLATION ET APPROXIMATION


from scipy import interpolate
import numpy as np

# === Interpolation 1D ===

# Données
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

# Linéaire
f_linear = interpolate.interp1d(x, y)
x_new = np.linspace(0, 4, 20)
y_new = f_linear(x_new)

# Différents types d'interpolation
f_nearest = interpolate.interp1d(x, y, kind='nearest')
f_linear = interpolate.interp1d(x, y, kind='linear')
f_slinear = interpolate.interp1d(x, y, kind='slinear')  # spline degré 1
f_quadratic = interpolate.interp1d(x, y, kind='quadratic')
f_cubic = interpolate.interp1d(x, y, kind='cubic')

# Splines polynomiales d'ordre arbitraire
f_spline3 = interpolate.interp1d(x, y, kind=3)  # Ordre 3
f_spline5 = interpolate.interp1d(x, y, kind=5)  # Ordre 5

# Avec extrapolation
f_extrap = interpolate.interp1d(x, y, fill_value='extrapolate')
y_outside = f_extrap(5)  # Extrapole au-delà de x=4

# Remplissage personnalisé hors limites
f_fill = interpolate.interp1d(x, y, bounds_error=False, fill_value=-999)
y_outside = f_fill([-1, 5])  # Retourne -999

# Remplissage asymétrique
f_fill = interpolate.interp1d(x, y, bounds_error=False, 
                              fill_value=(-1, 100))  # Gauche, droite

# === Splines cubiques ===

# CubicSpline - spline cubique C²
cs = interpolate.CubicSpline(x, y)
y_new = cs(x_new)

# Conditions aux bords
cs_natural = interpolate.CubicSpline(x, y, bc_type='natural')  # y''=0 aux bords
cs_clamped = interpolate.CubicSpline(x, y, bc_type='clamped')  # y'=0 aux bords
cs_notaknot = interpolate.CubicSpline(x, y, bc_type='not-a-knot')  # Défaut
cs_periodic = interpolate.CubicSpline(x, y, bc_type='periodic')  # Périodique

# Conditions personnalisées
cs_custom = interpolate.CubicSpline(x, y, bc_type=((1, 0), (1, 0)))
# (order, value): (1, 0) = première dérivée = 0

# Dérivées de la spline
y_prime = cs(x_new, 1)  # Première dérivée
y_second = cs(x_new, 2)  # Deuxième dérivée
y_third = cs(x_new, 3)  # Troisième dérivée

# Anti-dérivée (intégrale)
y_integral = cs.antiderivative()(x_new)

# Racines
roots = cs.roots()

# === Splines d'Akima ===

# Akima1DInterpolator - robuste aux outliers
akima = interpolate.Akima1DInterpolator(x, y)
y_new = akima(x_new)

# Avec extrapolation
akima_extrap = interpolate.Akima1DInterpolator(x, y, extrapolate=True)

# === Splines PCHIP ===

# Piecewise Cubic Hermite Interpolating Polynomial
# Préserve monotonie
pchip = interpolate.PchipInterpolator(x, y)
y_new = pchip(x_new)

# Dérivées
y_prime = pchip.derivative()(x_new)
y_integral = pchip.antiderivative()(x_new)

# Racines
roots = pchip.roots()

# === B-splines univariées ===

# BSpline - contrôle total
t = [0, 0, 0, 0, 1, 2, 3, 4, 4, 4, 4]  # Knots (avec multiplicité)
c = [0, 1, 4, 9, 16, 25, 36]  # Coefficients
k = 3  # Degré

spl = interpolate.BSpline(t, c, k)
y_new = spl(x_new)

# Créer B-spline par interpolation
tck = interpolate.splrep(x, y, s=0)  # s=0: passe par tous les points
# tck = (t, c, k) - tuple knots, coefficients, degré
y_new = interpolate.splev(x_new, tck)

# Lissage (s > 0)
tck_smooth = interpolate.splrep(x, y, s=1.0)  # Facteur lissage
y_smooth = interpolate.splev(x_new, tck_smooth)

# Avec poids
w = np.array([1, 2, 1, 2, 1])  # Poids pour chaque point
tck = interpolate.splrep(x, y, w=w, s=0)

# Dérivées
y_prime = interpolate.splev(x_new, tck, der=1)
y_second = interpolate.splev(x_new, tck, der=2)

# Intégrale
integral = interpolate.splint(0, 4, tck)

# Racines
roots = interpolate.sproot(tck)

# === UnivariateSpline - interface OO ===

# Interpolation exacte
spl = interpolate.UnivariateSpline(x, y, s=0)
y_new = spl(x_new)

# Lissage
spl_smooth = interpolate.UnivariateSpline(x, y, s=1.0)

# InterpolatedUnivariateSpline - toujours s=0
spl_interp = interpolate.InterpolatedUnivariateSpline(x, y)

# LSQUnivariateSpline - knots spécifiés
knots_interior = [1.5, 2.5]  # Knots intérieurs uniquement
spl_lsq = interpolate.LSQUnivariateSpline(x, y, knots_interior)

# Dérivées et intégrales
y_prime = spl.derivative()(x_new)
y_prime2 = spl.derivative(n=2)(x_new)
y_integral = spl.antiderivative()(x_new)

# Racines
roots = spl.roots()

# === Interpolation 2D ===

# Grille régulière
x_grid = np.linspace(0, 4, 5)
y_grid = np.linspace(0, 4, 5)
X, Y = np.meshgrid(x_grid, y_grid)
Z = X**2 + Y**2

# RectBivariateSpline - grille rectangulaire
f_2d = interpolate.RectBivariateSpline(x_grid, y_grid, Z)

# Évaluation
x_new = np.linspace(0, 4, 20)
y_new = np.linspace(0, 4, 20)
Z_new = f_2d(x_new, y_new)  # Grille complète

# Point unique
z_point = f_2d(1.5, 2.5)

# Dérivées partielles
dz_dx = f_2d(x_new, y_new, dx=1, dy=0)  # ∂z/∂x
dz_dy = f_2d(x_new, y_new, dx=0, dy=1)  # ∂z/∂y
d2z_dxdy = f_2d(x_new, y_new, dx=1, dy=1)  # ∂²z/∂x∂y

# Intégrale
integral = f_2d.integral(0, 4, 0, 4)  # ∫∫ f(x,y) dx dy

# interp2d - interface simple (deprecated, utiliser RectBivariateSpline)
# f_2d = interpolate.interp2d(x_grid, y_grid, Z, kind='cubic')

# RegularGridInterpolator - N-dimensions, moderne
from scipy.interpolate import RegularGridInterpolator

points = (x_grid, y_grid)
values = Z
interp = RegularGridInterpolator(points, values, method='linear')

# Évaluation sur points quelconques
pts = np.array([[1.5, 2.5], [2.0, 3.0]])
z_interp = interp(pts)

# Méthodes disponibles
interp_linear = RegularGridInterpolator(points, values, method='linear')
interp_nearest = RegularGridInterpolator(points, values, method='nearest')
interp_slinear = RegularGridInterpolator(points, values, method='slinear')
interp_cubic = RegularGridInterpolator(points, values, method='cubic')
interp_quintic = RegularGridInterpolator(points, values, method='quintic')

# === Interpolation points irréguliers (scattered data) ===

# Points aléatoires
np.random.seed(42)
n_points = 100
x_random = np.random.rand(n_points) * 4
y_random = np.random.rand(n_points) * 4
z_random = x_random**2 + y_random**2

# griddata - interpolation simple
X_new, Y_new = np.meshgrid(np.linspace(0, 4, 50), np.linspace(0, 4, 50))
points_scattered = np.column_stack([x_random, y_random])
points_regular = np.column_stack([X_new.ravel(), Y_new.ravel()])

Z_interp = interpolate.griddata(points_scattered, z_random, 
                                (X_new, Y_new), method='linear')
# method='nearest', 'linear', 'cubic'

# LinearNDInterpolator - interpolation linéaire Delaunay
interp_linear = interpolate.LinearNDInterpolator(points_scattered, z_random)
Z_linear = interp_linear(X_new, Y_new)

# Avec valeur par défaut hors convex hull
interp_linear = interpolate.LinearNDInterpolator(points_scattered, z_random,
                                                 fill_value=0)

# NearestNDInterpolator - plus proche voisin
interp_nearest = interpolate.NearestNDInterpolator(points_scattered, z_random)
Z_nearest = interp_nearest(X_new, Y_new)

# CloughTocher2DInterpolator - C1 continu, piecewise cubic
interp_ct = interpolate.CloughTocher2DInterpolator(points_scattered, z_random)
Z_ct = interp_ct(X_new, Y_new)

# Avec rescaling
interp_ct = interpolate.CloughTocher2DInterpolator(points_scattered, z_random,
                                                    rescale=True)

# === Radial Basis Function (RBF) ===

# RBFInterpolator - moderne (scipy >= 1.7)
from scipy.interpolate import RBFInterpolator

rbf = RBFInterpolator(points_scattered, z_random, kernel='thin_plate_spline')
Z_rbf = rbf(points_regular).reshape(X_new.shape)

# Kernels disponibles
# 'linear', 'thin_plate_spline', 'cubic', 'quintic',
# 'multiquadric', 'inverse_multiquadric', 'inverse_quadratic', 'gaussian'

rbf_gaussian = RBFInterpolator(points_scattered, z_random, 
                               kernel='gaussian', epsilon=1)

# Lissage
rbf_smooth = RBFInterpolator(points_scattered, z_random, smoothing=0.1)

# Rbf - ancienne interface (deprecated)
# rbf_old = interpolate.Rbf(x_random, y_random, z_random, function='thin_plate')

# === Interpolation paramétrique (courbes) ===

# Courbe 2D paramétrée
t = np.linspace(0, 2*np.pi, 10)
x = np.cos(t)
y = np.sin(t)

# splprep/splev - spline paramétrique
tck, u = interpolate.splprep([x, y], s=0)  # u = paramètre
u_new = np.linspace(0, 1, 100)
x_new, y_new = interpolate.splev(u_new, tck)

# Courbe 3D
z = t
tck_3d, u = interpolate.splprep([x, y, z], s=0)
x_new, y_new, z_new = interpolate.splev(u_new, tck_3d)

# Courbe fermée (périodique)
tck_periodic, u = interpolate.splprep([x, y], s=0, per=True)

# Avec lissage
tck_smooth, u = interpolate.splprep([x, y], s=1.0)

# Dérivées
dx_du, dy_du = interpolate.splev(u_new, tck, der=1)

# === BarycentricInterpolator - interpolation de Lagrange ===

# Interpolation barycentrique (stable numériquement)
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

bary = interpolate.BarycentricInterpolator(x, y)
y_new = bary(x_new)

# Ajout de points
bary.add_xi([5], [25])

# Poids personnalisés
weights = np.array([1, 1, 1, 1, 1])
bary = interpolate.BarycentricInterpolator(x, y, wi=weights)

# === KroghInterpolator - polynôme de Hermite ===

# Interpolation avec dérivées
krogh = interpolate.KroghInterpolator(x, y)
y_new = krogh(x_new)

# Dérivées
y_prime = krogh.derivative(x_new, der=1)

# === Approximation polynomiale ===

# Ajustement polynôme moindres carrés
degree = 3
coeffs = np.polyfit(x, y, degree)
poly = np.poly1d(coeffs)
y_fitted = poly(x_new)

# Avec poids
weights = np.array([1, 2, 1, 2, 1])
coeffs = np.polyfit(x, y, degree, w=weights)

# Évaluation avec np.polyval
y_fitted = np.polyval(coeffs, x_new)

# === Approximation rationnelle ===

# pade - approximation de Padé (ratio polynômes)
from scipy.interpolate import pade

# Coefficients série Taylor
coeffs = [1, 1, 0.5, 1/6, 1/24]  # exp(x) ≈ 1 + x + x²/2 + ...
p, q = pade(coeffs, 2)  # Padé [2/2]
# p, q sont des poly1d (numérateur, dénominateur)

# Évaluation
x_eval = 0.5
y_approx = p(x_eval) / q(x_eval)


[OK] SCIPY.LINALG - ALGÈBRE LINÉAIRE


from scipy import linalg
import numpy as np

# === Opérations basiques ===

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Inverse
A_inv = linalg.inv(A)

# Déterminant
det_A = linalg.det(A)

# Résolution système linéaire Ax = b
b = np.array([5, 11])
x = linalg.solve(A, b)

# Vérification
assert np.allclose(A @ x, b)

# Système multiple (plusieurs seconds membres)
B_multi = np.array([[5, 1], [11, 2]])
X = linalg.solve(A, B_multi)

# Avec matrice singulière (pseudo-inverse)
x_lstsq, residuals, rank, s = linalg.lstsq(A, b)

# === Décompositions matricielles ===

# LU decomposition
P, L, U = linalg.lu(A)
# PA = LU

# Sans permutation (retourne seulement L et U)
lu, piv = linalg.lu_factor(A)
x = linalg.lu_solve((lu, piv), b)  # Plus efficace pour multiples b

# QR decomposition
Q, R = linalg.qr(A)
# A = QR, Q orthogonale, R triangulaire supérieure

# Mode économique
Q, R = linalg.qr(A, mode='economic')

# Pivoting
Q, R, P = linalg.qr(A, pivoting=True)

# Cholesky (matrice symétrique définie positive)
A_sym = np.array([[4, 2], [2, 3]])
L = linalg.cholesky(A_sym, lower=True)  # A = LL^T
U = linalg.cholesky(A_sym, lower=False)  # A = U^TU

# Résolution avec Cholesky (plus rapide)
c, lower = linalg.cho_factor(A_sym)
x = linalg.cho_solve((c, lower), b)

# SVD - Singular Value Decomposition
U, s, Vt = linalg.svd(A)
# A = U @ diag(s) @ Vt

# Reconstruction
A_reconstructed = U @ np.diag(s) @ Vt

# SVD tronquée (économique)
U, s, Vt = linalg.svd(A, full_matrices=False)

# Valeurs singulières seulement
s = linalg.svdvals(A)

# Schur decomposition
T, Z = linalg.schur(A)
# A = Z @ T @ Z^H, T quasi-triangulaire

# Forme réelle de Schur
T_real, Z_real = linalg.schur(A, output='real')

# Hessenberg form
H, Q = linalg.hessenberg(A, calc_q=True)
# A = Q @ H @ Q^H, H Hessenberg

# === Valeurs et vecteurs propres ===

# Valeurs propres
eigenvalues = linalg.eigvals(A)

# Valeurs et vecteurs propres
eigenvalues, eigenvectors = linalg.eig(A)
# A @ v = λ @ v

# Vérification
lambda_0 = eigenvalues[0]
v_0 = eigenvectors[:, 0]
assert np.allclose(A @ v_0, lambda_0 * v_0)

# Valeurs propres seulement (plus rapide)
eigenvalues = linalg.eigvals(A)

# Matrice symétrique/hermitienne (plus rapide et stable)
A_sym = np.array([[1, 2], [2, 1]])
eigenvalues_sym, eigenvectors_sym = linalg.eigh(A_sym)

# Sous-ensemble de valeurs propres
eigenvalues_subset = linalg.eigh(A_sym, eigvals=(0, 0))  # Plus petite seulement

# Valeurs propres dans intervalle
eigenvalues_range, eigenvectors_range = linalg.eigh(
    A_sym, subset_by_value=(-10, 10)
)

# Indices spécifiques
eigenvalues_indices, eigenvectors_indices = linalg.eigh(
    A_sym, subset_by_index=[0, 1]  # 2 plus petites
)

# Problème aux valeurs propres généralisé: Av = λBv
B = np.array([[2, 1], [1, 2]])
eigenvalues_gen, eigenvectors_gen = linalg.eig(A, B)

# Symétrique généralisé
eigenvalues_gen_sym, eigenvectors_gen_sym = linalg.eigh(A_sym, B)

# === Normes matricielles ===

# Norme de Frobenius
norm_fro = linalg.norm(A, 'fro')

# Norme nucléaire (somme valeurs singulières)
norm_nuc = linalg.norm(A, 'nuc')

# Norme induite
norm_2 = linalg.norm(A, 2)  # Norme spectrale (max valeur singulière)
norm_1 = linalg.norm(A, 1)  # Max somme colonnes
norm_inf = linalg.norm(A, np.inf)  # Max somme lignes
norm_neg_inf = linalg.norm(A, -np.inf)  # Min somme lignes
norm_neg_1 = linalg.norm(A, -1)  # Min somme colonnes

# Condition number
cond = linalg.cond(A)  # κ(A) = ||A|| ||A^(-1)||
cond_1 = linalg.cond(A, 1)
cond_inf = linalg.cond(A, np.inf)

# === Matrices spéciales ===

# Matrice identité
I = np.eye(3)

# Matrice nulle
Z = np.zeros((3, 3))

# Matrice diagonale
d = np.array([1, 2, 3])
D = np.diag(d)

# Matrice triangulaire
A_tri = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
L_tri = np.tril(A_tri)  # Triangulaire inférieure
U_tri = np.triu(A_tri)  # Triangulaire supérieure

# Bloc diagonal
from scipy.linalg import block_diag
A1 = np.array([[1, 2], [3, 4]])
A2 = np.array([[5]])
A3 = np.array([[6, 7], [8, 9]])
BD = block_diag(A1, A2, A3)

# Matrice compagnon
coeffs = [1, -6, 11, -6]  # x³ - 6x² + 11x - 6
C = linalg.companion(coeffs)

# Matrice de Hadamard
H = linalg.hadamard(4)  # Taille doit être 2^k

# Matrice de Hilbert
from scipy.linalg import hilbert, invhilbert
H_hilbert = hilbert(5)  # H_ij = 1/(i+j-1)
H_inv = invhilbert(5)  # Inverse exacte

# Matrice de Toeplitz
from scipy.linalg import toeplitz
c = [1, 2, 3, 4]
r = [1, 5, 6, 7]
T = toeplitz(c, r)  # Première colonne, première ligne

# Circulante (cas spécial Toeplitz)
from scipy.linalg import circulant
circ = circulant([1, 2, 3, 4])

# Matrice de Hankel
from scipy.linalg import hankel
c = [1, 2, 3, 4]
r = [4, 5, 6]
H_hankel = hankel(c, r)

# Matrice de Leslie (démographie)
from scipy.linalg import leslie
f = [0, 2, 1]  # Fécondité
s = [0.5, 0.3]  # Survie
L = leslie(f, s)

# === Fonctions de matrices ===

# Exponentielle matricielle
exp_A = linalg.expm(A)

# Logarithme matriciel
log_A = linalg.logm(A)

# Puissance matricielle
A_power = linalg.fractional_matrix_power(A, 0.5)  # Racine carrée
sqrt_A = linalg.sqrtm(A)  # Équivalent

# Sinus, cosinus
sin_A = linalg.sinm(A)
cos_A = linalg.cosm(A)
tan_A = linalg.tanm(A)

# Hyperboliques
sinh_A = linalg.sinhm(A)
cosh_A = linalg.coshm(A)
tanh_A = linalg.tanhm(A)

# Fonction matricielle générale
def f(x):
    return x**2 + 2*x + 1

f_A = linalg.funm(A, f)

# Exponentielle avec multiplication par vecteur (efficace)
v = np.array([1, 1])
result = linalg.expm_multiply(A, v)  # exp(A) @ v

# Évolution temporelle: exp(t*A) @ v
t = 0.5
result_t = linalg.expm_multiply(A, v, start=0, stop=t, num=10)

# === Équations matricielles ===

# Sylvester: AX + XB = Q
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
Q = np.array([[1, 1], [1, 1]])
X = linalg.solve_sylvester(A, B, Q)

# Lyapunov: AX + XA^T = Q (cas spécial Sylvester)
X_lyap = linalg.solve_lyapunov(A, Q)

# Équation algébrique de Riccati continue
A = np.array([[0, 1], [-1, -1]])
B = np.array([[0], [1]])
Q = np.eye(2)
R = np.array([[1]])
X_riccati = linalg.solve_continuous_are(A, B, Q, R)

# Équation algébrique de Riccati discrète
X_riccati_d = linalg.solve_discrete_are(A, B, Q, R)

# Lyapunov continue: AX + XA^T + Q = 0
X_lyap_c = linalg.solve_continuous_lyapunov(A, Q)

# Lyapunov discrète: AXA^T - X + Q = 0
X_lyap_d = linalg.solve_discrete_lyapunov(A, Q)

# === Pseudo-inverse et projections ===

# Pseudo-inverse de Moore-Penrose
A_pinv = linalg.pinv(A)

# Avec tolérance
A_pinv = linalg.pinv(A, rcond=1e-10)

# Avec méthode spécifique
A_pinv = linalg.pinv(A, rtol=1e-10)

# Projection orthogonale sur colonne-space
orth_basis = linalg.orth(A)

# Null space (noyau)
null_basis = linalg.null_space(A)

# === Matrices bandes ===

# Matrice bande (stockage efficace)
ab = np.array([[0, 0, 1, 2, 3],   # Diagonale supérieure 2
               [0, 4, 5, 6, 7],   # Diagonale supérieure 1
               [8, 9, 10, 11, 12]])  # Diagonale principale

# Résolution système bande
l = 0  # Sous-diagonales
u = 2  # Sur-diagonales
b = np.array([1, 2, 3, 4, 5])
x = linalg.solve_banded((l, u), ab, b)

# Décomposition LU bande
lu_band, piv = linalg.lu_factor_banded(ab, l, u)

# Cholesky bande (symétrique définie positive)
ab_sym = np.array([[4, 5, 6, 7],  # Diagonale principale
                   [1, 2, 3, 0]])  # Sous-diagonale
c_band = linalg.cholesky_banded(ab_sym, lower=True)

# === Matrices triangulaires ===

# Résolution système triangulaire (plus rapide)
L = np.tril(A)
x = linalg.solve_triangular(L, b, lower=True)

U = np.triu(A)
x = linalg.solve_triangular(U, b, lower=False)

# Inverse triangulaire
L_inv = linalg.inv_triangular(L, lower=True)

# === Matrices de rotation ===

# Rotation 2D
theta = np.pi / 4
R_2d = np.array([[np.cos(theta), -np.sin(theta)],
                 [np.sin(theta), np.cos(theta)]])

# Rotation 3D autour axe
from scipy.spatial.transform import Rotation as R
r = R.from_euler('z', 45, degrees=True)
R_3d = r.as_matrix()

# === Interpolation matricielle ===

# Racine carrée matricielle par Schur
sqrt_A_schur = linalg.sqrtm(A)

# Logarithme matriciel par eigenvalues
log_A_eig = linalg.logm(A)

# === Matrices creuses (voir scipy.sparse) ===

# Conversion dense -> sparse
from scipy import sparse
A_sparse = sparse.csr_matrix(A)

# === Utilitaires ===

# Rang matriciel
rank = linalg.matrix_rank(A)

# Trace
trace = np.trace(A)

# Vérifier définition positive
try:
    linalg.cholesky(A_sym)
    is_positive_definite = True
except linalg.LinAlgError:
    is_positive_definite = False

# Vérifier symétrie
is_symmetric = np.allclose(A, A.T)

# Orthonormaliser vecteurs (Gram-Schmidt)
vectors = np.random.rand(5, 3)  # 3 vecteurs de dimension 5
ortho = linalg.orth(vectors)


[OK] SCIPY.SPARSE - MATRICES CREUSES


from scipy import sparse
import numpy as np

# === Formats de matrices creuses ===

# COO - Coordinate format (construction)
row = np.array([0, 0, 1, 2, 2])
col = np.array([0, 2, 1, 0, 2])
data = np.array([1, 2, 3, 4, 5])
coo = sparse.coo_matrix((data, (row, col)), shape=(3, 3))

# CSR - Compressed Sparse Row (calculs par ligne)
csr = sparse.csr_matrix((data, (row, col)), shape=(3, 3))

# CSC - Compressed Sparse Column (calculs par colonne)
csc = sparse.csc_matrix((data, (row, col)), shape=(3, 3))

# Depuis matrice dense
dense = np.array([[1, 0, 2], [0, 3, 0], [4, 0, 5]])
csr_from_dense = sparse.csr_matrix(dense)

# LIL - List of Lists (construction incrémentale)
lil = sparse.lil_matrix((3, 3))
lil[0, 0] = 1
lil[0, 2] = 2
lil[1, 1] = 3
lil[2, 0] = 4
lil[2, 2] = 5

# DOK - Dictionary of Keys (construction incrémentale)
dok = sparse.dok_matrix((3, 3))
dok[0, 0] = 1
dok[0, 2] = 2

# DIA - Diagonal (matrices bandes)
diagonals = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
offsets = [-1, 0, 1]
dia = sparse.diags(diagonals, offsets, shape=(3, 3))

# BSR - Block Sparse Row (matrices blocs)
indptr = np.array([0, 2, 3, 6])
indices = np.array([0, 2, 2, 0, 1, 2])
data = np.array([1, 2, 3, 4, 5, 6]).repeat(4).reshape(6, 2, 2)
bsr = sparse.bsr_matrix((data, indices, indptr), shape=(6, 6))

# === Conversions ===

# Conversion entre formats
csr_to_csc = csr.tocsc()
csr_to_coo = csr.tocoo()
csr_to_lil = csr.tolil()
csr_to_dok = csr.todok()
csr_to_dia = csr.todia()
csr_to_bsr = csr.tobsr()

# Vers dense
dense_from_sparse = csr.toarray()
dense_matrix = csr.todense()  # Retourne np.matrix (deprecated)

# === Opérations basiques ===

A_sparse = sparse.csr_matrix([[1, 0, 2], [0, 3, 0], [4, 0, 5]])
B_sparse = sparse.csr_matrix([[2, 0, 1], [0, 4, 0], [1, 0, 3]])

# Addition
C = A_sparse + B_sparse

# Soustraction
C = A_sparse - B_sparse

# Multiplication scalaire
C = 2 * A_sparse

# Multiplication matricielle
C = A_sparse @ B_sparse
C = A_sparse.dot(B_sparse)

# Multiplication élément par élément
C = A_sparse.multiply(B_sparse)

# Puissance élément par élément
C = A_sparse.power(2)

# Transposée
A_T = A_sparse.T
A_T = A_sparse.transpose()

# Conjuguée transposée
A_H = A_sparse.H
A_H = A_sparse.conjugate().transpose()

# === Propriétés ===

# Dimensions
print(A_sparse.shape)
print(A_sparse.ndim)

# Nombre d'éléments non-nuls
nnz = A_sparse.nnz

# Densité
density = A_sparse.nnz / (A_sparse.shape[0] * A_sparse.shape[1])

# Indices éléments non-nuls
rows, cols = A_sparse.nonzero()
values = A_sparse.data

# === Accès éléments ===

# Lecture
value = A_sparse[1, 1]

# Écriture (LIL ou DOK seulement, efficace)
lil_matrix = A_sparse.tolil()
lil_matrix[1, 1] = 10
A_sparse = lil_matrix.tocsr()

# Slicing
submatrix = A_sparse[0:2, 1:3]

# === Construction matrices spéciales ===

# Identité
I = sparse.eye(5)
I_csr = sparse.eye(5, format='csr')

# Diagonale
d = np.array([1, 2, 3, 4, 5])
D = sparse.diags(d)

# Diagonales multiples
diag_values = [[1, 2, 3], [4, 5, 6, 7], [8, 9, 10]]
offsets = [-1, 0, 1]
M = sparse.diags(diag_values, offsets, shape=(4, 4))

# Matrice aléatoire creuse
random_sparse = sparse.random(100, 100, density=0.05, format='csr')

# Avec distribution spécifique
from scipy.stats import norm
random_normal = sparse.random(100, 100, density=0.05, 
                              data_rvs=norm().rvs, format='csr')

# Bloc diagonal
from scipy.sparse import block_diag
A1 = sparse.csr_matrix([[1, 2], [3, 4]])
A2 = sparse.csr_matrix([[5]])
BD = block_diag([A1, A2])

# Matrice Kronecker
kron = sparse.kron(A_sparse, B_sparse)

# Matrice verticale/horizontale
vstack = sparse.vstack([A_sparse, B_sparse])
hstack = sparse.hstack([A_sparse, B_sparse])

# Matrices bandes (tridiagonale)
n = 5
diagonals = [np.ones(n-1), -2*np.ones(n), np.ones(n-1)]
tridiag = sparse.diags(diagonals, [-1, 0, 1], format='csr')

# === Algèbre linéaire creuse ===

from scipy.sparse import linalg as sparse_linalg

# Résolution système linéaire Ax = b
A = sparse.csr_matrix([[3, 0, 1], [0, 2, 0], [1, 0, 3]])
b = np.array([1, 2, 3])

# Factorisation LU directe (petites matrices)
x = sparse_linalg.spsolve(A, b)

# Solveur itératif (grandes matrices)
x, info = sparse_linalg.cg(A, b)  # Gradient conjugué (symétrique)
x, info = sparse_linalg.gmres(A, b)  # GMRES (général)
x, info = sparse_linalg.bicg(A, b)  # BiCG
x, info = sparse_linalg.bicgstab(A, b)  # BiCGSTAB (plus stable)
x, info = sparse_linalg.minres(A, b)  # MINRES (symétrique)
x, info = sparse_linalg.lgmres(A, b)  # LGMRES

# Avec préconditionnement
M = sparse_linalg.spilu(A.tocsc())  # ILU preconditioner
M_op = sparse_linalg.LinearOperator(A.shape, M.solve)
x, info = sparse_linalg.gmres(A, b, M=M_op)

# Avec callback
def callback(xk):
    print(f"Itération: résidu = {np.linalg.norm(A @ xk - b)}")

x, info = sparse_linalg.gmres(A, b, callback=callback)

# Factorisation LU creuse
lu = sparse_linalg.splu(A.tocsc())
x = lu.solve(b)

# Factorisation ILU (incomplete LU)
ilu = sparse_linalg.spilu(A.tocsc())
x = ilu.solve(b)

# Factorisation Cholesky (symétrique définie positive)
# Nécessite sksparse ou cholmod
# from sksparse.cholmod import cholesky
# factor = cholesky(A)
# x = factor(b)

# === Valeurs propres creuses ===

# Quelques valeurs propres
n_eigenvalues = 3

# Plus grandes valeurs propres (module)
eigenvalues, eigenvectors = sparse_linalg.eigs(A, k=n_eigenvalues)

# Plus petites valeurs propres (module)
eigenvalues, eigenvectors = sparse_linalg.eigs(A, k=n_eigenvalues, 
                                               which='SM')

# which: 'LM' (largest magnitude), 'SM' (smallest magnitude),
#        'LR' (largest real), 'SR' (smallest real),
#        'LI' (largest imaginary), 'SI' (smallest imaginary)

# Avec valeur cible (proches de sigma)
eigenvalues, eigenvectors = sparse_linalg.eigs(A, k=n_eigenvalues, 
                                               sigma=1.0, which='LM')

# Matrices symétriques/hermitiennes (plus rapide)
A_sym = sparse.csr_matrix([[2, 1, 0], [1, 2, 1], [0, 1, 2]])
eigenvalues, eigenvectors = sparse_linalg.eigsh(A_sym, k=n_eigenvalues)

# Intervalle de valeurs propres
eigenvalues, eigenvectors = sparse_linalg.eigsh(A_sym, k=n_eigenvalues,
                                                which='SA')  # Smallest algebraic
# which: 'LA' (largest algebraic), 'SA' (smallest algebraic),
#        'BE' (both ends)

# Valeurs singulières
u, s, vt = sparse_linalg.svds(A, k=n_eigenvalues)

# Plus grandes valeurs singulières
u, s, vt = sparse_linalg.svds(A, k=n_eigenvalues, which='LM')

# === Normes creuses ===

# Norme de Frobenius
norm_fro = sparse_linalg.norm(A, 'fro')

# Norme infinie
norm_inf = sparse_linalg.norm(A, np.inf)

# === Opérateurs linéaires ===

# LinearOperator - matrice implicite
def matvec(v):
    return A @ v

def rmatvec(v):
    return A.T @ v

A_op = sparse_linalg.LinearOperator((3, 3), matvec=matvec, rmatvec=rmatvec)