Frattalita' in misure di spessore (2)

Aggiornamento bibliografico (che mi ero completamente saltato all'Universita')

Sediment Accumulation Rates and the Completeness of Stratigraphic Sections

 Peter M. Sadler

The Journal of Geology, Vol. 89, No. 5 (Sep., 1981), pp. 569-584

Estimation of completeness of stratigraphical sections using
empirical data and theoretical models  

Peter M.Sadler & David J. Strauss 

Journal of the Geological Society, London, Vol. 147, 1990,pp. 471-485

 

The nature of stratigraphical record (3 ed) 

D. Ager (1993) 

Wiley and Sons 

capitolo 3 More gaps than records pag.43 

 

Fractals models in Earth Sciences

Korvin Gabor (1992)

Elsevier Science Ltd

cap.2 Fractals in Flatland : a romance in <2 dimensions 

 ----------------------------------------------------------------------------- 

 

Seconda puntata dell'esperimento iniziato qui

 Su suggerimento di uno degli autori ho provato ad utilizzare i dati di questo articolo in cui sono presenti dettagliate stratigrafie del Scaglia Toscana

The Scaglia Toscana Formation of the Monti del Chianti:
new lithostratigraphic and biostratigraphic data
Enrico Pandeli (*, **), Milvio Fazzuoli (*), Fabio Sandrelli (***), Roberto Mazzei (***),
Simonetta Monechi (*), Marisa Nocchi (****), Ivan Martini (***) & Gigliola Valleri (*) Ital. J. Geosci., Vol. 137, No. 1 (2018), pp. 38-61, 12 figs., 1 tab., 1 plate, 2 extrapl. (doi: 10.3301/IJG.2017.16)
© Società Geologica Italiana, Roma 2018

 

 Sono state considerati considerati gli affioramenti Sugame (48 dati Middle Eocene-Lower Oligocene) , Castello di Brolio (251 dati Upper Cretaceous-Lower Oligocene) , Cintoia West (150 dati Upper Cretaceous-Middle Upper Eocene ) , Montegrossi (208 dati) e Lucolena North (Upper Cretaceous--Lower Oligocene 63 dati). Come si vede non tutte le serie ricoprono esattamente lo stesso periodo temporale e le stess formazioni. Si evidenzia che alcune serie non presentano una serie deposizionale continua

In generale tutti i grafici mostrano una curva spezzata con valori di D>1 sulle scale di maggiori dimensioni e valori nettamente inferiori a 1 per le scale di dettaglio. Alle scale maggiori emerge quindi una tendenza alla frattalita' ed autosomiglianza che si perde alle piccola scala

Il motivo di questo andamento e' perche' i cicli tettonici e climatici influenzano solo le scali maggiori (ma potrebbe essere anche un problema di sotto campionamento degli strati meno potenti)

Questo e' il codice gentilmente fornito da Gemini 

 

import numpy as np

import matplotlib.pyplot as plt

import os

def box_counting_dimension_1d(data, max_box_size=None, plot=True):

    """

    Calcola la dimensione frattale di una serie di dati 1D utilizzando il metodo del box-counting.

    Args:

        data (np.array): La serie di dati 1D (spessore di strato) come array NumPy.

        max_box_size (int, optional): La dimensione massima delle scatole da testare.

        plot (bool): Se True, visualizza il log-log plot.

    Returns:

        float: La dimensione frattale calcolata (pendenza del log-log plot).

    """

    # 1. Preparazione dei dati (la logica rimane la stessa)

    data = np.asarray(data)

   

    # Rimuovi eventuali valori non finiti (NaN, Inf) che potrebbero derivare dal caricamento

    data = data[np.isfinite(data)]

    if len(data) < 2:

        print("Errore: Il set di dati è troppo piccolo o vuoto dopo la pulizia.")

        return np.nan

    min_val = np.min(data)

    max_val = np.max(data)

   

    if max_val == min_val:

        return 1.0

    L = len(data)

   

    if max_box_size is None:

        max_box_size = L // 4

   

    # Raccogli le dimensioni delle scatole (epsilon) da testare

    r_values = 2**np.arange(1, int(np.log2(max_box_size)) + 1)

   

    # Scalatura dell'asse Y per mappare i punti in uno spazio 2D approssimativamente quadrato

    y_scaled = (data - min_val) / (max_val - min_val) * (L - 1)

    points = np.column_stack((np.arange(L), y_scaled))

    log_r = []

    log_N = []

   

    # 3. Box-Counting

    for r in r_values:

        # Calcola le coordinate delle scatole

        box_coords = np.floor(points / r).astype(int)

       

        # Conta le scatole uniche

        N_r = len(np.unique(box_coords, axis=0))

       

        if N_r > 0:

            log_r.append(np.log(1/r))

            log_N.append(np.log(N_r))

           

    log_r = np.array(log_r)

    log_N = np.array(log_N)

   

    if len(log_r) < 2:

        print("Errore: Non abbastanza punti per calcolare la pendenza. Prova con una serie di dati più lunga.")

        return np.nan

       

# Nuova sezione: Selezione dei primi 4 punti

    num_punti_fit = 3

    if len(log_r) < num_punti_fit:

        print(f"Attenzione: Trovati solo {len(log_r)} punti. Eseguo la regressione su tutti i punti disponibili.")

        log_r_fit = log_r

        log_N_fit = log_N

    else:

        # Seleziona i primi 'num_punti_fit' punti (corrispondenti alle scale più piccole/fini)

        #log_r_fit = log_r[-num_punti_fit:]

        #log_N_fit = log_N[-num_punti_fit:]

        log_r_fit = log_r[:num_punti_fit]

        log_N_fit = log_N[:num_punti_fit]

   

    # Esegue la regressione lineare solo sui punti selezionati

    coeffs, residuals, rank, singular_values, rcond = np.polyfit(log_r_fit, log_N_fit, 1, full=True)

    fractal_dimension = coeffs[0]

    # ----------------------------------------------------

   

    # 5. Plot (Opzionale)

    if plot:

        plt.figure(figsize=(10, 6))

        # Plot di TUTTI i punti

        plt.plot(log_r, log_N, 'bo', label='Punti Dati Totali')

       

        # Plot della linea di regressione utilizzando solo l'intervallo di punti selezionato

        p = np.poly1d(coeffs)

       

        # Plot della linea SOLO sul range dei punti usati per il fit (opzionale, per chiarezza)

        plt.plot(log_r_fit, p(log_r_fit), 'r-',

                 label=f'Fit (Primi {len(log_r_fit)} punti): D = {fractal_dimension:.4f}')

        # Evidenzia i punti usati per il fit

        plt.plot(log_r_fit, log_N_fit, 'ro', markerfacecolor='none', markersize=10,

                 label=f'Punti usati per la regressione')

        plt.title('Sugame')

        plt.xlabel(r'$\log(1/\epsilon)$')

        plt.ylabel(r'$\log(N(\epsilon))$')

        plt.legend()

        plt.grid(True, which="both", linestyle="--")

        plt.show()

    return fractal_dimension

# --- Esempio di Utilizzo AGGIORNATO ---

NOME_FILE_CSV = "sugame.csv"

print("--- Calcolo della Dimensione Frattale ---")

if not os.path.exists(NOME_FILE_CSV):

    print(f"ERRORE: File '{NOME_FILE_CSV}' non trovato.")

    print("Assicurati che il file esista nella stessa directory dello script.")

else:

    try:

        # Carica i dati dal file CSV

        # usecols=0 specifica di usare solo la prima colonna

        # dtype=int assicura che i dati vengano letti come interi

        spessore_strato_dati = np.genfromtxt(

            NOME_FILE_CSV,

            delimiter=',',

            dtype=int,

            skip_header=0,

            filling_values=np.nan # Sostituisce i valori mancanti con NaN (che puliremo dopo)

        )

       

        if spessore_strato_dati.ndim > 1:

            # Se per qualche motivo numpy legge più colonne, forziamo solo la prima

            spessore_strato_dati = spessore_strato_dati[:, 0]

       

        print(f"Dati caricati con successo da '{NOME_FILE_CSV}'.")

        print(f"Numero di punti dati: {len(spessore_strato_dati)}")

        # Calcola e stampa il risultato

        D = box_counting_dimension_1d(spessore_strato_dati, plot=True)

        print(f"\nRisultato:")

        print(f"Dimensione Frattale (D) calcolata: {D:.4f}")

    except Exception as e:

        print(f"Si è verificato un errore durante il caricamento o l'analisi dei dati: {e}")

 Ben sapendo che non e' propriamente corretto ho provato a vedere cosa succedeva sommando tutti i dati di tutte le serie

 

 

 

Esponente di Hurst su torbidita' di acquiferi apuani

In questo post vengono analizzati i valori di torbidita' della sorgente del Frigido (periodo 13/02/2020 01/12/2025 dati orari ) e Carbonera (26/07/2006 21/06/2013 dati ogni 3 minuti ) legate all'acquiferi carsico delle Alpi Apuane

 

Torbidita' del tempo Frigido

 

Il parametro Esponente di Hurst e' stato introdotto Hurst, H. E. (1951). Long-term storage capacity of reservoirs. Transactions of the American Society of Civil Engineers, 116(1), 770-808. 

E' una misura dell'inerzia del sistema idrogeologico (ovvero capacita' di stoccaggio) e della complessita' della rete carsica. Per i valore di 0.82 e 0.85 trovati si ha che il sistema non e' casuale (H=0.5), non e' deterministicamente semplice ed il comportamento e' correlato su differenti scale di tempo ed e' un sistema che ha "memoria" del proprio passato

Non e' possibile prevedere in maniera completamente deterministica il valore esatto del parametro indagato  

Frigido

Con modelli deterministici si avra' un discreto successo su scale di tempo brevi ma non sulla previsione a lunga scadenza. Il passaggio da regime deterministico a frattale e' indicato come tempo di crossover e si individua dai punti di flesso sull'andamento dei punti. Per Carbonera (grafico sottostante) il punto di cross over e'  tra 1000 e 10000 unita' che tradotto nel passo di campionamento si arriva ad intervallo da le 5 e le 50 ore

Carbonera

Valori di H vicini ad 1 indicano persistenza nei valori (sia che siano massimi che minimi) e questo puo' incidere anche sul valore massimo e non solo sulle code della distribuzione

La memoria non si esaurisce mai ma ha una legge di decadimento come una potenze (in questo caso 2*0.84 -2 = -0.32 

 

con questa legge e l'esponente calcolate la correlazione ha un peso superiore al  10% fino a 1000 intervalli temporali (rho = 0.10^(1/-0.32)

Il fatto che il calcolo di H risulti sostanzialmente identico su serie con diverso passo di campionamento e' un indice della frattalita' del fenomeno osservato (e' indipendente dalla scala del tempo di osservazione)

 

il valore di Hurst prossimo ad 1 indica che e' probabile che la torbidita' aumenti nel prossimo lag (passo di campionamento) e' fortemente dipendente da cio' che accade nel lag precedente . Cio' sarebbe in contrasto con il concetto di massimo relativo della curva di scarica ma la statistica deve essere considerata valida al mantenersi delle condizioni al contorno e delle relative forzanti. Con un valore di 0.82 la coda di esaurimento sara' lunga

 

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

# --- 1. FUNZIONE DEFINIZIONE per l'ANALISI R/S MANUALE ---

def compute_rs_data(series):

    """

    Calcola la statistica Rango Riscalato (R/S) per una serie temporale.

    Restituisce gli intervalli temporali (n) e i valori medi di R/S.

    """

    N = len(series)

    # Definisce gli intervalli temporali (n) in progressione geometrica fino a N/4

    n_values = np.unique(np.logspace(1, np.log10(N // 4), num=50, dtype=int))

    n_values = n_values[n_values >= 1]

   

    rs_values = []

   

    for n in n_values:

        if n == 0: continue

       

        num_blocks = N // n

        rs_per_block = []

       

        for i in range(num_blocks):

            block = series[i * n : (i + 1) * n]

           

            mean = np.mean(block)

            std_dev = np.std(block)

           

            if std_dev == 0:

                continue

            X = np.cumsum(block - mean)

            R = np.max(X) - np.min(X)

           

            rs_per_block.append(R / std_dev)

       

        if rs_per_block:

            rs_values.append(np.mean(rs_per_block))

        else:

            break

    return n_values[:len(rs_values)], np.array(rs_values)

# --- 2. CARICAMENTO E PREPARAZIONE DEI DATI ADATTATI ---

FILE_NAME = "frigido.csv"

COLUMN_NAME = "Torbidita NTU"

FILTER_N_THRESHOLD = 100

try:

    # ATTENZIONE: Caricamento con separatore PUNTO E VIRGOLA (;) e decimale VIRGOLA (,)

    df = pd.read_csv(FILE_NAME, sep=';', decimal=',')

   

    # Seleziona la colonna 'Torbidita NTU'

    discharge_series = df[COLUMN_NAME]

   

    # Conversione in numerico e gestione dei NaN (valori mancanti o non validi)

    discharge_series = pd.to_numeric(discharge_series, errors='coerce').dropna()

   

    # Converti in array NumPy

    discharge_data = discharge_series.values

   

except FileNotFoundError:

    print(f"Errore: Il file '{FILE_NAME}' non è stato trovato. Assicurati che sia nella directory corretta.")

    exit()

except KeyError:

    print(f"Errore: La colonna '{COLUMN_NAME}' non è stata trovata nel file. Controlla il nome esatto.")

    print(f"Colonne disponibili: {df.columns.tolist()}")

    exit()

if len(discharge_data) < 100:

     print("Errore: Non ci sono abbastanza punti dati validi per eseguire l'analisi R/S.")

     exit()

# --- 3. CALCOLO R/S ---

rs_n, rs_rs = compute_rs_data(discharge_data)

# Calcolo di H (Esponente di Hurst) per l'INTERVALLO COMPLETO

log_n_full = np.log10(rs_n)

log_rs_full = np.log10(rs_rs)

coefficients_full = np.polyfit(log_n_full, log_rs_full, 1)

H_full = coefficients_full[0]

c_full = 10**coefficients_full[1]

# --- 4. PLOT: ANALISI R/S INTERVALLO COMPLETO ---

plt.figure(figsize=(8, 5))

plt.plot(rs_n, c_full * rs_n**H_full, color="deepskyblue", label=f'Pendenza Frattale (H={H_full:.2f})')

plt.scatter(rs_n, rs_rs, color="navy", alpha=0.5, label='Dati R/S empirici')

plt.xscale('log')

plt.yscale('log')

plt.xlabel('Intervallo Temporale (n)')

plt.ylabel('Rango Riscalato R/S')

plt.title(f'Scaling Frattale della Torbidità (Full Range)')

plt.legend()

plt.savefig('Frigido_Full_R_S_Plot.png')

# --- 5. PLOT: INTERVALLO FILTRATO (n > 100) ---

filter_mask = rs_n > FILTER_N_THRESHOLD

X_filtered = rs_n[filter_mask]

Y_filtered = rs_rs[filter_mask]

# Calcolo di H per il solo intervallo filtrato

if len(X_filtered) > 1:

    log_n_filtered = np.log10(X_filtered)

    log_rs_filtered = np.log10(Y_filtered)

   

    coefficients_filtered = np.polyfit(log_n_filtered, log_rs_filtered, 1)

    H_filtered = coefficients_filtered[0]

    c_filtered = 10**coefficients_filtered[1]

   

    plt.figure(figsize=(8, 5))

    plt.scatter(X_filtered, Y_filtered, color="darkgreen", alpha=0.7, label=f'Dati R/S empirici (n > {FILTER_N_THRESHOLD})')

    plt.plot(X_filtered, c_filtered * X_filtered**H_filtered, color="red", linestyle='--',

             label=f'Pendenza Frattale Filtrata (H={H_filtered:.2f})')

    plt.xscale('log')

    plt.yscale('log')

    plt.xlabel('Intervallo Temporale (n)')

    plt.ylabel('Rango Riscalato R/S')

    plt.title(f'Scaling Frattale della Torbidità (Intervallo Temporale > {FILTER_N_THRESHOLD})')

    plt.legend()

    plt.savefig('Frigido_Filtered_R_S_Plot_gt_100.png')

   

    print("\n--- Risultati del Calcolo ---")

    print(f"Hurst Esponent (Torbidità - Intervallo Completo): {H_full:.4f}")

    print(f"Hurst Esponent (Torbidità - Filtrato n>{FILTER_N_THRESHOLD}): {H_filtered:.4f}")

else:

    print(f"Nota: Non ci sono abbastanza punti dati rimasti dopo il filtraggio (n > {FILTER_N_THRESHOLD}) per generare il grafico filtrato.")

 

 

Montagne frattali

 Negli anni 90 quando i miei amici matematici leggevano (Ri)creazioni al calcolatore di Dewdney (vedi anche pag. 107 del libro The Magic Machine a handbook of computer sorcery  A.K.Dewdney 1990 oppure Le Scienze n 222 1987 pag 98 e succ)  andava di moda programmare il proprio algoritmo di morphing e di  creazioni di paesaggi frattali

L'algoritmo che mi avevano spiegato era 

1) partire da un triangolo

2) cercare un punto casuale interno al triangolo 

3) uniire i vertici del triangolo con il punto casuale interno ed iterare sui nuovi tre triangoli

 All'epoca non ero riuscito a fare il mio codice...oggi ho deciso di chiedera all'AI. la cosa interessante del codice proposto e' l'algoritmo con cui viene cercato il punto casuale interno nella funzione random_point_inside. Vengono generati due numeri casuali compresi tra 0 e 1. Se la somma e' maggiore di 1 viene ridotta perche' il sistema esplora in realta' rettangoli, se la somma e' maggiore di 1 si e' selezionato un punto non nel triangolo

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.tri import Triangulation

# ------------------------------------------------------

# 1. SCEGLI UN TRIANGOLO DI PARTENZA

# ------------------------------------------------------

# Triangolo grande in coordinate 2D

points = [

    np.array([0.0, 0.0]),

    np.array([1.0, 0.0]),

    np.array([0.5, 0.9])

]

triangles = [(0,1,2)]     # un singolo triangolo iniziale

points = list(points)

# ------------------------------------------------------

# 2. FUNZIONE: scegli un punto interno casuale (barycentric random)

# ------------------------------------------------------

def random_point_inside_triangle(a,b,c):

    r1 = np.random.rand()

    r2 = np.random.rand()

    # riflessione per distribuzione uniforme

    if r1 + r2 > 1:

        r1 = 1 - r1

        r2 = 1 - r2

    return a + r1*(b-a) + r2*(c-a)

# ------------------------------------------------------

# 3. ITERAZIONE DI SUDDIVISIONE

# ------------------------------------------------------

def subdivide_random(triangles, points, iterations=500):

    for _ in range(iterations):

        # scegli casualmente un triangolo da suddividere

        idx = np.random.randint(len(triangles))

        i,j,k = triangles[idx]

        a,b,c = points[i], points[j], points[k]

        # nuovo punto interno

        p = random_point_inside_triangle(a,b,c)

        p_index = len(points)

        points.append(p)

        # crea tre nuovi triangoli

        new_tris = [

            (i, j, p_index),

            (j, k, p_index),

            (k, i, p_index)

        ]

        # sostituisci il triangolo originale con i 3 nuovi

        triangles[idx] = new_tris[0]

        triangles.extend(new_tris[1:])

    return triangles, points

triangles, points = subdivide_random(triangles, points, iterations=1500)

# ------------------------------------------------------

# 4. PLOT FINALE

# ------------------------------------------------------

pts = np.array(points)

tri = np.array(triangles)

plt.figure(figsize=(7,7))

plt.tripcolor(pts[:,0], pts[:,1], tri, edgecolors='k', linewidth=0.1, alpha=0.9)

plt.title("Random Triangle Subdivision – Paesaggio frattale 2D")

plt.axis("equal")

plt.show()

 

Un altro algoritmo (diamond square) e' basato su una griglia regolare lavorando sui vertici. Il modello e' frattale peri il metodo con cui viene calcolata la altezza  

 

 

 

 

 import numpy as np
import matplotlib.pyplot as plt

def midpoint_displacement(size, roughness=0.7):
    """
    Genera un terreno frattale 2D (matrice di elevazione Z)
    utilizzando l'algoritmo Diamond-Square semplificato (Midpoint Displacement).
    
    :param size: Il lato della griglia, deve essere (2^n) + 1.
    :param roughness: Controlla la "frattalità" (0=liscio, alto=irregolare).
    :return: Una matrice NumPy di elevazioni (Z).
    """
    # Inizializza la griglia Z con zeri
    Z = np.zeros((size, size), dtype=float)
    
    # 1. Inizializzazione degli angoli
    # Assegna valori casuali iniziali ai quattro angoli per dare il via al processo
    Z[0, 0] = np.random.rand() * 10
    Z[0, size - 1] = np.random.rand() * 10
    Z[size - 1, 0] = np.random.rand() * 10
    Z[size - 1, size - 1] = np.random.rand() * 10
    
    step = size - 1
    rand_range = 10.0  # Range iniziale della casualità
    
    # Processo iterativo
    while step > 1:
        half = step // 2
        
        # Quadrati (Diamond Step)
        # Calcola il punto medio di ogni quadrato (i, j)
        for i in range(0, size - 1, step):
            for j in range(0, size - 1, step):
                
                # Media dei quattro angoli del quadrato
                avg = (Z[i, j] + Z[i + step, j] + Z[i, j + step] + Z[i + step, j + step]) / 4.0
                
                # Aggiunge casualità al centro del quadrato
                Z[i + half, j + half] = avg + (np.random.rand() - 0.5) * 2 * rand_range
        
        # Diamanti (Square Step)
        # Calcola i punti medi dei bordi (punti 'diamante')
        for i in range(0, size - 1, half):
            for j in range((i + half) % step, size - 1, step):
                
                # Punti da mediare (gestione bordi per non sforare l'array)
                punti = [
                    Z[(i - half + size) % size, j], # Sopra
                    Z[(i + half) % size, j],         # Sotto
                    Z[i, (j - half + size) % size], # Sinistra
                    Z[i, (j + half) % size]          # Destra
                ]
                
                # Media dei quattro vicini (quelli che esistono)
                avg = np.mean(punti)
                
                # Aggiunge casualità al punto medio del bordo
                Z[i, j] = avg + (np.random.rand() - 0.5) * 2 * rand_range
                
                # Questo codice semplificato è per i bordi interni; un'implementazione
                # completa dovrebbe gestire meglio gli indici di confine
        
        # Riduce l'ampiezza di casualità per la prossima iterazione (l'algoritmo frattale)
        rand_range *= 2**(-roughness)
        step = half
        
    return Z

def visualize_wireframe(Z):
    """
    Visualizza la matrice di elevazione Z come un wireframe 3D.
    """
    size = Z.shape[0]
    
    # Creazione delle coordinate X e Y
    x = np.linspace(0, 1, size)
    y = np.linspace(0, 1, size)
    X, Y = np.meshgrid(x, y)

    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')

    # Plot del wireframe
    ax.plot_wireframe(X, Y, Z, color='brown', linewidth=0.3)

    # Impostazioni estetiche
    ax.set_title(f'Montagna Frattale Wireframe (Roughness={ROUGHNESS})')
    ax.set_xlabel('Asse X')
    ax.set_ylabel('Asse Y')
    ax.set_zlabel('Elevazione Z')
    
    # Regola gli assi Z per centrare meglio la montagna
    z_min, z_max = np.min(Z), np.max(Z)
    z_range = z_max - z_min
    ax.set_zlim(z_min - z_range * 0.1, z_max + z_range * 0.1)

    plt.show()

# --- Parametri di Esecuzione ---
# La dimensione (size) DEVE essere (2^n) + 1, ad esempio: 33, 65, 129, 257.
# 65 è un buon compromesso tra dettaglio e velocità.
SIZE = 65 
ROUGHNESS = 0.85 # Valori più alti (es. 0.9) danno montagne più "spigolose"

# 2. Generazione del Frattale e Visualizzazione
elevation_data = midpoint_displacement(SIZE, ROUGHNESS)
visualize_wireframe(elevation_data)