#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 28 11:34:19 2025

@author: jean-mariefurbringer
"""

import numpy as np

def doehlert(k: int) -> np.ndarray:
    """
    Génère une matrice de Doehlert pour k facteurs.

    Paramètres
    ----------
    k : int
        Nombre de facteurs.

    Retour
    ------
    Mexp : np.ndarray
        Matrice des expériences Doehlert (N x k)
    """
    # Calcul du nombre d’expériences
    # N = k * (k + 1) + 1

    # Préparation de la matrice Z (taille (k+1) x (k+1))
    Z = np.eye(k + 1)
    Z = Z - np.hstack((np.ones((k + 1, 1)), np.zeros((k + 1, k))))

    # Préparation de C : nombre de lignes = (k+1)*k - k = k^2 ; colonnes = k+1
    C = np.zeros((k * k, k + 1))
    i = 0

    # Boucles traduites fidèlement en considérant les indices MATLAB 1-based
    for l1 in range(1, k + 2):          # l1 = 1 .. k+1
        s = 1
        for m in range(2, k + 2):       # m = 2 .. k+1
            ll = l1 + m - s
            if ll > (k + 1):
                # correspond exactement au code MATLAB
                ll = 2
                s = l1 + m - 2
            if l1 != ll:
                # conversion vers 0-based pour indexer Z
                C[i, :] = Z[l1 - 1, :] - Z[ll - 1, :]
                i += 1

    # concatène Z et C
    Z = np.vstack((Z, C))

    # Préparation de la base orthogonale B (taille (k+1) x (k+1))
    B = np.zeros((k + 1, k + 1))

    # sous-diagonale -1 (équivalent de -diag(ones(k,1), -1) en MATLAB)
    for idx in range(1, k + 1):
        B[idx, idx - 1] = -1.0

    # remplissage B(l,ll) = 1/ll pour l=1..k, ll=l..k (attention 1-based)
    for l1 in range(1, k + 1):
        for ll in range(l1, k + 1):
            B[l1 - 1, ll - 1] = 1.0 / ll

    # dernière colonne = ones(k+1,1)
    B[:, k] = 1.0

    # Calcul de la matrice d'expérience X puis mise à l'échelle des colonnes
    X = Z @ B

    # mise à l'échelle : pour l = 1..k appliquer sqrt(l / (2*(l+1)))
    for l1 in range(1, k + 1):
        scale = np.sqrt(l1 / (2.0 * (l1 + 1)))
        X[:, l1 - 1] *= scale

    Mexp = X[:, :k]
    return Mexp

M = doehlert(5)
print(M.shape)   # doit être ((3*(3+1)+1) x 3) => (13, 3)
print(np.round(M, 3))