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)

Connettore N (disavventure di un radioamatore)

Pensavo che le antenne per radioamatori avessero connettori SMA, BNC e SO 239 (connettore HF femmina) e PL259 (connettore HF maschio)

Comprando una antenna Diamond x50n ho scoperto che esistono i connettori N

 

 

 

Frattalita' in spessore di strato

Sempre leggendo i libri dei post precedenti ho visto che e' stata trovata una correlazione frattale nello spessore degli strati turbiditici arenacei. Avendo fatto la tesi sulle Arenarie del Cervarola mi sono incuriosito ma ho misure di sezioni o affioramenti

Gli unici dati che ho trovato in merito a spessori di strato sono relativi agli stratotipo ed in particolare dell'affioramento di Punta Piccola per il Piacenziano (The Global Standard Stratotype-section and Point (GSSP) of the Piacenzian Stage (Middle Pliocene) by D. Castradori, D. Rio, F. J. Hilgen, and L. J. Lourens https://doi.org/10.18814/epiiugs/1998/v21i2/003.Non si tratta di una formazione turbiditica ma non avendo altri dati ci ho provato 

Questa e' la distribuzione di frequenza degli spessori..un limite e' che la distribuzione non e' oltre un ordine di grandezza 

 

 plottando log-log lo spessore degli strati con il numero di strati con spessore superiore alla soglia si ha che

 

usando lo spettro di potenza 

il valore β≈0.47 e' compreso tra il rumore bianco ed ed il fattore di scala 1/f dove β=1 (Rothman, D. H., Grotzinger, J. P., & Richter, F. M. (1994). Scaling in the Sedimentary Rock Record. Nature, 372(6505), 337–342. suggerisce che il valore atteso per le sequenze stratigrafiche e' β=1 altrimenti detto Pink Noise )

Cerchero' dati piu' accurati

 

Frattalita' in misure geofisiche di pozzo (2)

Nella conclusione del post precedente era che i vari tentativi di ricondurre i dati sperimentali all'analisi frattale davano risultati contradittori...tutto cio' fino sono arrivato a leggere Fractal Behaviour of the Earth System V.P.Dimitri  Springer 2005 pag 63-81 Regularity Analysis Applied to Well Log Data (M.Fedi, D.Fiore, Mauro La Mann) dove parlando proprio dei dati di KST viene introdotto il concetto di multifrattalita' (per essere precisi di multifrattalita' aveva gia' parlato Mandelbrot ed in D. Turcotte  Fractals and Chaos in Geology and Geophysics si trova applicazione dei multifrattali applicata alle geoscienze a pag)

Prima di tutto si puo' effettuare un test se i dati del pozzo possono essere inquadrati all'interno della multifrattalita' 

Questo puo' essere fatto l'analisi di regolarita' locale tramite analisi di fluttuazione a momenti multipli (Multifractal Detrended Fluctuation Analysis, MF-DFA) 

Vengono testati valori del parametro q da -5 a 5 (escludendo lo 0) verificando la pendenza h(q) 

 

 

Se il segnale fosse monofrattale la curva sarebbe una linea orizzontale ma non e' il caso. Quindi il test mostra la presenza di trend multifrattale con q>0 (destra) influenzato dalle grandi fluttuazioni (picchi del grafico dei valori di Gamma o contrasti gelogici) o q<0 (sinistra) influenzati da regioni piu' omogenee o da rumore di fondo

La linea rossa (q=2) indica  la persistenza media dell'intera serie

nel grafico sottostante viene descritta l'intensita' della multifrattalita' (se non fossero dati multi frattali si ridurrebbe ad un punto. Piu' ampia e' l'apertura della curva maggiore e' l'intensita' della frattalita'

I dati quindi non possono essere descritti con un singolo esponente di Hurst 

 

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from scipy.stats import linregress

# ======================================================

# 0 — FUNZIONE MF-DFA MANUALE

# ======================================================

def mfdfa(series, qs, scales, order=1):

"""

Calcola l'esponente di Hurst generalizzato h(q) utilizzando l'algoritmo MF-DFA.

"""

N = len(series)

# 1. Integrazione del segnale (necessaria per il detrending)

Y = np.cumsum(series - np.mean(series))

h_q = []

for q in qs:

F_q_2 = [] # Momenti di fluttuazione di q-esimo ordine

for s in scales:

s = int(s)

# Suddivisione in segmenti di lunghezza s

segments = N // s

# Calcolo della varianza per tutti i segmenti

F_q_s_list = []

for i in range(segments):

idx = slice(i * s, (i + 1) * s)

Y_segment = Y[idx]

# 2. Detrending locale (regressione polinomiale di ordine 'order')

x_local = np.arange(s)

p = np.polyfit(x_local, Y_segment, order)

Y_fit = np.polyval(p, x_local)

# 3. Calcolo della varianza locale (residui al quadrato)

variance = np.mean((Y_segment - Y_fit)**2)

F_q_s_list.append(variance)

# 4. Calcolo della funzione di fluttuazione di q-esimo ordine

F_q_s_list = np.array(F_q_s_list)

if q != 2.0:

# Esecuzione del momento q-esimo (q/2 sulle varianze)

F_q = np.mean(F_q_s_list**(q / 2.0))**(1.0 / q)

else:

# Caso speciale q=2 (momento quadratico)

F_q = np.sqrt(np.mean(F_q_s_list))

F_q_2.append(F_q)

# 5. Regressione log-log: F_q(s) ~ s^h(q)

log_s = np.log10(scales)

log_F_q = np.log10(F_q_2)

# Filtra i valori infiniti o NaN che possono emergere da log(0)

valid_indices = np.isfinite(log_s) & np.isfinite(log_F_q)

if np.sum(valid_indices) > 1:

slope, _, _, _, _ = linregress(log_s[valid_indices], log_F_q[valid_indices])

h_q.append(slope)

else:

h_q.append(np.nan)

return np.array(h_q)

# ======================================================

# 1 — LETTURA FILE CSV E PREPARAZIONE DATI

# ======================================================

filename = "well_log.csv"

RIGHE_LIMITE = 12700

COLONNA_INDICE = 1

try:

df = pd.read_csv(filename, nrows=RIGHE_LIMITE)

except FileNotFoundError:

print(f"Errore: File '{filename}' non trovato.")

exit()

x = df.iloc[:, COLONNA_INDICE].values

x = pd.to_numeric(x, errors='coerce')

x = x[~np.isnan(x)]

# Riscalamento (opzionale)

xmin = np.min(x)

if xmin <= 0:

shift = abs(xmin) + 1e-6

x = x + shift

print(f"Riscalamento effettuato: shift = {shift:.6f}")

N = len(x)

print(f"Analisi eseguita su {N} campioni (Prime {RIGHE_LIMITE} righe del file).")

# ======================================================

# 2 — PARAMETRI MF-DFA

# ======================================================

# Scelta dell'intervallo q

Q_MIN = -5

Q_MAX = 5

Q_STEP = 0.5

qs = np.arange(Q_MIN, Q_MAX + Q_STEP, Q_STEP)

qs = qs[qs != 0] # Esclude q=0

# Scegliere un intervallo di scale logaritmico per una migliore regressione

min_scale = int(N / 50) # Scale minime per un detrending efficace (p.es. 50 campioni)

max_scale = int(N / 4) # Scala massima (1/4 della lunghezza totale)

scales = np.unique(np.logspace(np.log10(min_scale), np.log10(max_scale), num=50).astype(int))

ORDER = 1 # Ordine del polinomio di detrending locale

# ======================================================

# 3 — ESECUZIONE MF-DFA E CALCOLO h(q)

# ======================================================

print("Calcolo degli esponenti h(q)...")

hq = mfdfa(x, qs, scales, order=ORDER)

# Filtra i valori NaN (dovuti a una regressione fallita)

valid_hq = hq[np.isfinite(hq)]

valid_qs = qs[np.isfinite(hq)]

# ======================================================

# 4 — ANALISI E INTERPRETAZIONE DEI RISULTATI

# ======================================================

# 4a. Analisi monofattale (q=2)

if 2.0 in valid_qs:

q_mono_index = np.where(valid_qs == 2.0)[0][0]

hurst_monofractal = valid_hq[q_mono_index]

else:

# Interpolazione per q=2 se non è stato calcolato

hurst_monofractal = np.interp(2.0, valid_qs, valid_hq)

# 4b. Analisi della Multifrattalità

h_max = np.max(valid_hq)

h_min = np.min(valid_hq)

Delta_h = h_max - h_min

# ======================================================

# 5 — PLOT DEI RISULTATI

# ======================================================

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=False)

# --- A) Grafico h(q) vs q (Test di Multifattalità)

ax1.plot(valid_qs, valid_hq, 'o-', color='darkblue', markersize=4)

ax1.axhline(hurst_monofractal, color='red', linestyle='--',

label=f'H monofattale (q=2): {hurst_monofractal:.3f}')

ax1.set_title('Esponente di Hurst Generalizzato $h(q)$')

ax1.set_xlabel('Momento di Fluttuazione $q$')

ax1.set_ylabel('Esponente di Scaling $h(q)$')

ax1.grid(True)

ax1.legend()

# --- B) Lo Spettro Multifrattale f(alpha) vs alpha (Singolarità)

# Calcolo dello Spettro tramite Trasformata di Legendre

tau_q = valid_qs * valid_hq - 1

alpha = np.gradient(tau_q, valid_qs)

f_alpha = valid_qs * alpha - tau_q

ax2.plot(alpha, f_alpha, 's-', color='green', markersize=4)

ax2.set_title('Spettro Multifrattale $f(\\alpha)$')

ax2.set_xlabel('Esponente di Singolarità $\\alpha$')

ax2.set_ylabel('Dimensione Frattale $f(\\alpha)$')

ax2.grid(True)

plt.tight_layout()

plt.show()

print("\n" + "=" * 50)

print(f"ANALISI MULTIFRATTALE (MF-DFA) - RISULTATI")

print(f"Δh = h_max - h_min = {Delta_h:.4f}")

if Delta_h > 0.1:

print("CONCLUSIONE: Il segnale è **fortemente multifrattale** (Δh > 0.1).")

print("La sua complessità non è uniforme: sono presenti strutture a diverse scale e diversi tipi di singolarità (anomalie e omogeneità).")

else:

print("CONCLUSIONE: Il segnale è prevalentemente **monofattale** (Δh <= 0.1).")

print(f"Esponente di Hurst Monofattale (q=2): H ≈ {hurst_monofractal:.3f}")

print("=" * 50)

 

Frattalita' in misure geofisiche di pozzo

Premessa: dopo questo post mi sono lanciato in questa avventura nei frattali nella geologia non avendo nessuna nei guai in cui mi sarei infilato, di quando fosse malevola la natura che non solo e' frattale ma multifrattale e come sarei stato inadeguato dal punto di vista matematico. Detto questo l'argomento e' cosi complesso che viene diviso in due post,il presente ed il successivo 

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

Sempre leggendo Fractals and Chaos in Geology and Geophysics 2nd edition -- Donald Lawson Turcotte -- 2nd ed_, Cambridge, U_K, New York, England, 1997 si ha che le misure geofisiche seguono una legge frattale

Per fare una prova mi sono scaricato un log di pozzo da  KGS

I dati sono in formato LAS che non indicano il formato per laserscanner ma un formato testuale per log di pozzo per possono essere letti con la libreria lasio di Python

import lasio

# Load LAS file

las = lasio.read("1055385870.las")

# Convert to pandas DataFrame

df = las.df()

# Save to CSV

df.to_csv("well_log.csv", index=True)

io ho scaricato del tutto a caso il file 1055385870.las queste le prime righe

 ~Version Information
 VERS.                2.0:  CWLS LOG ASCII STANDARD -VERSION 2.0
 WRAP.                 NO:  ONE LINE PER DEPTH STEP
~WELL INFORMATION BLOCK
#MNEM.UNIT       DATA                    DESCRIPTION OF MNEMONIC
#---------    -------------            ------------------------
 STRT.FT  0.300                          :START DEPTH
 STOP.FT  1279.300                       :STOP DEPTH
 STEP.FT  0.100                          :STEP    UP_HOLE  
 NULL.    -999.25                        :NULL VALUE
 COMP.    LEIS, VICTOR J.                :COMPANY
 WELL.    MENTZER #N-1                                                 :WELL
 FLD .     NEOSHO FALLS - LEROY          :FIELD\LOCATION
 LOC .    26    23S   17E                :LOCATION
 CNTY.    ALLEN                          :COUNTY
 STAT.    KANSAS                         :STATE
 CTRY.    USA                            :COUNTRY
 SRVC.                                   :SERVICE COMPANY
 DATE.    05/06/2011                     :LOG DATE
 API .   15- 001-30189-00-00                 :API NUMBER
 LIC .                                   :LICENSE NUMBER
~Curve Information Block
#MNEM.UNIT      API CODE     Curve Description
#---------    -------------  ------------------------
 DEPT     .FT        00 001 00 00          :  1  DEPTH   
 COND(DI) .MMHO/M    00 000 00 00          :  2  DEEP INDUCTION COND 
 RES(GRD) .OHM-M     00 220 00 00          :  3  GUARD RESISTIVITY   
 SP       .MV        00 010 00 00          :  4  SPON POTENTIAL      
 RES(DI)  .OHM-M     00 120 00 00          :  5  DEEP INDUCTION RES  
 RES(MI)  .OHM-M     00 120 00 00          :  6  MED INDUCTION RES   
 GAMMA    .API-GR    00 310 00 00          :  7  GAMMA RAY           
 COMP     .G/CC      00 356 00 00          :  8  DEN COMPENSATION    
 CALIPER  .INCH      00 280 00 00          :  9  SHORT ARM CALIPER   
 DEN(CDL) .G/CC      00 350 00 00          : 10  COMPENSATED DENSITY 
 TEMP     .DEG F     00 660 00 00          : 11  TEMPERATURE         
 POR(DEN) .PERCENT   00 890 00 00          : 12  DENSITY POROSITY    
 DEL TEMP .DEG F     00 000 00 00          : 13  SAMPLE DELTA TEMP   
~Parameter Information Block
#MNEM.UNIT   Information     Description
#---------- ---------------  ------------------------
 FILE.      PROCESSED          :File Type
 FIID.      9541               :File Type Identifier
 VERS.      3.60B              :System Version
 SER .      1                  :System Serial Number
 TRUK.      .3879              :Truck Calibration Number

 

Ho scelto di lavorare sul parametro Gamma. 

Ci sono metodi per evidenziare la frattalita'

1) e' basato sull'analisi della frequenza delle classi di una misura

2) e' basata  S(f)∝f^-β

Le ipotesi di base per la verifica della frattalita' sono 

a) non e' necessariamente stazionario, puo' avere un trend e la varianza puo' crescere nel tempo ma lo sono le differenze tra X(t) e X(t+1). Cio' implica l'autosomiglianza ovvero lo stesso comportamento a diverse scale

b) deve sottostare ad un legge di potenza. Questa puo ' essere intepretata come 

    -  frequenza come potenza N(>x)∝x^(-alpha). la somma di elementi superiori ad una soglia e' direttamente proporzionale al valore della soglia alla potenza di -alpha

    - power spectral power : diretta proporzionalita' tra la potenza di spettro ad una determinata frequenza e la frequenza elevata alla potenze di -Beta     PSD(f)∝f^-β

Proviamo con i dati reali  partendo prima dalla frequenza

Nel grafico soprastante ci sono le classi con il numero di misure inferiori ad una determinata soglia

  

Se si plotta in grafico log-log si ha un chiaro trend lineare nella parte finale del grafico 

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression

# ======================================================

# 1 — LETTURA FILE CSV E COLONNA DI INTERESSE

# ======================================================

filename = "well_log.csv"

df = pd.read_csv(filename)

# Usa la colonna RES(GRD) → indice 2

col_index = 6

x = df.iloc[:, col_index].values

print(x)

# RIMUOVI i NaN (prime 100 righe vuote)

x = x[~np.isnan(x)]

if len(x) < 10:

raise ValueError("Troppi valori NaN in RES(GRD). Servono più dati validi.")

# ======================================================

# 2 — RISCALAMENTO SE CI SONO VALORI NEGATIVI

# ======================================================

xmin_original = np.min(x)

if xmin_original <= 0:

shift = abs(xmin_original) + 1e-6

x = x + shift

print(f"Riscalamento effettuato: aggiunto shift = {shift:.6f}")

# ======================================================

# 3 — CREAZIONE DI 25 CLASSI TRA MINIMO E MASSIMO

# ======================================================

xmin = np.min(x)

xmax = np.max(x)

thresholds = np.linspace(xmin, xmax, 25)

# ======================================================

# 4 — PER OGNI SOGLIA, CONTA GLI ELEMENTI >= SOGLIA

# ======================================================

counts = np.array([np.sum(x >= thr) for thr in thresholds])

# FILTRO valori validi

valid = np.logical_and(counts > 0, thresholds > 0)

thresholds = thresholds[valid]

counts = counts[valid]

# controlla se sono rimasti valori

if len(thresholds) < 3:

raise ValueError("Troppi valori invalidi dopo il filtraggio. Controlla la colonna RES(GRD).")

# ======================================================

# 5 — LOGARITMI

# ======================================================

logT = np.log10(thresholds)

logN = np.log10(counts)

# ======================================================

# 6 — REGRESSIONE SOLO PER log(threshold) ∈ [2.25, 2.5]

# ======================================================

mask = (logT >= 1.9) & (logT <= 4)

if np.sum(mask) >= 2:

logT_fit = logT[mask].reshape(-1, 1)

logN_fit = logN[mask]

model = LinearRegression()

model.fit(logT_fit, logN_fit)

slope = model.coef_[0]

intercept = model.intercept_

print("\n=== RISULTATI REGRESSIONE ===")

print("Intervallo usato: log10(soglia) ∈ [1.9, 4.0]")

print(f"Coefficiente angolare (slope) = {slope}")

print(f"Intercetta = {intercept}")

logN_pred = model.predict(logT_fit)

else:

print("\n⚠ Nessun dato nel range log10(t) ∈ [1.9, 4.0].")

print("La regressione NON è stata eseguita.")

logT_fit, logN_fit, logN_pred = None, None, None

# ======================================================

# 7 — ISTOGRAMMA DELLE CLASSI

# ======================================================

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

plt.bar(range(len(counts)), counts, width=0.8)

plt.xlabel("Classe (0–24) Parametro Gamma")

plt.ylabel("Numero elementi ≥ soglia della classe")

plt.title("Istogramma delle 25 classi Parametro Gamma")

plt.grid(axis='y')

plt.show()

# ======================================================

# 8 — GRAFICO LOG–LOG CON REGRESSIONE

# ======================================================

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

plt.scatter(logT, logN, s=20, label="Dati (log-log)")

if logT_fit is not None:

plt.plot(logT_fit.flatten(), logN_pred, linewidth=2,

label="Regressione (2.25–2.5)")

plt.xlabel("log10(soglia)")

if logT_fit is not None:

# Plot regressione

plt.plot(logT_fit.flatten(), logN_pred, linewidth=2,

label="Regressione (2.25–2.5)")

# Testo sul grafico con slope (3 cifre significative)

txt = f"slope = {slope:.3g}"

# posizione: leggermente a destra del primo punto valido

x_pos = logT_fit.min()

y_pos = logN_pred.max()

plt.text(x_pos, y_pos, txt, fontsize=12, color="red")

plt.ylabel("log10(N ≥ soglia)")

plt.title("Analisi log–log con regressione su intervallo selezionato")

plt.legend()

plt.grid(True)

plt.show()

usando invece l'analisi di densita' spettrale

 

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from numpy.fft import rfft, rfftfreq

from scipy.signal import welch

# ======================================================

# CONFIGURAZIONE

# ======================================================

CSV_FILE = "well_log.csv" # <-- cambia questo

WINDOW = 400 # sliding window

STEP = 1

COLUMN_DEPTH = 0

COLUMN_RXORT = 6

# ======================================================

# 1. LETTURA CSV

# ======================================================

df = pd.read_csv(CSV_FILE)

depth = df.iloc[:, COLUMN_DEPTH].values

rxort = df.iloc[:, COLUMN_RXORT].values

# ======================================================

# 2. FILTRO VALORI INVALIDI

# ======================================================

bad_values = [-9999, -999.25]

valid_mask = (

~np.isnan(depth) &

~np.isnan(rxort) &

np.isfinite(depth) &

np.isfinite(rxort) &

~np.isin(rxort, bad_values)

)

depth = depth[valid_mask]

rxort = rxort[valid_mask]

print("Campioni dopo il filtro:", len(rxort))

# ======================================================

# 3. FREQUENZA DI CAMPIONAMENTO

# ======================================================

dz = np.mean(np.diff(depth))

fs = 1.0 / dz # samples per meter

# ======================================================

# 4. SPETTROGRAMMA FFT A FINESTRA MOBILE

# ======================================================

num_samples = len(rxort)

freqs = rfftfreq(WINDOW, d=1/fs)

spectra = []

centers = []

for i in range(0, num_samples - WINDOW, STEP):

segment = rxort[i:i+WINDOW]

segment = segment - np.mean(segment)

F = rfft(segment)

P = np.abs(F)**2

spectra.append(P)

centers.append(depth[i + WINDOW//2])

spectra = np.array(spectra)

centers = np.array(centers)

# ======================================================

# 5. PLOT SPETTROGRAMMA + LOG ORIGINALE

# ======================================================

fig = plt.figure(figsize=(12, 8))

gs = fig.add_gridspec(2, 1, height_ratios=[2, 1], hspace=0.25)

# ---- (A) Spectrogram ----

ax1 = fig.add_subplot(gs[0, 0])

im = ax1.imshow(

spectra.T,

aspect='auto',

origin='lower',

extent=[centers.min(), centers.max(), freqs.min(), freqs.max()]

)

ax1.set_ylabel("Frequency (cycles/m)")

ax1.set_title("Spectrogram (window = 400 samples)")

# plt.colorbar(im, ax=ax1).set_label("Power")

# ---- (B) Log ----

ax2 = fig.add_subplot(gs[1, 0], sharex=ax1)

ax2.plot(depth, rxort, linewidth=0.8)

ax2.set_xlabel("Depth (m)")

ax2.set_ylabel("")

ax2.grid(True)

plt.setp(ax1.get_xticklabels(), visible=False)

plt.show()

# ======================================================

# 6. PSD WELCH + FRACTAL FIT

# ======================================================

nper = min(1024, len(rxort))

freqs_psd, psd = welch(rxort, fs=fs, nperseg=nper)

print("PSD samples:", len(freqs_psd))

print("Freq range:", freqs_psd.min(), " - ", freqs_psd.max())

# ---- Intervallo utile per fit ----

mask = (freqs_psd > 0.001) & (freqs_psd < 0.2)

print("Punti nel fit:", np.sum(mask))

if np.sum(mask) < 2:

raise ValueError("Intervallo di frequenza troppo stretto o dati insufficienti!")

# ---- Fit lineare in scala log-log ----

p = np.polyfit(np.log(freqs_psd[mask]), np.log(psd[mask]), 1)

beta = -p[0]

H = (beta - 1) / 2

D = 2 - H

print("\n========== RISULTATI FRACTAL ==========")

print("Spectral β =", beta)

print("Hurst H =", H)

print("Fractal D =", D)

# ======================================================

# 7. PLOT PSD + FIT

# ======================================================

plt.loglog(freqs_psd, psd, label="PSD")

plt.loglog(freqs_psd[mask], np.exp(np.polyval(p, np.log(freqs_psd[mask]))),

label="Fit", linewidth=2)

plt.xlabel("Frequency (cycles/m)")

plt.ylabel("PSD")

plt.title("Spectral fractal analysis")

plt.legend()

plt.show()

 

Spectral β = -0.594 (pendenza della retta di approssimazione)

Per questo valore negativo indica che le basse frequenze hanno un maggior potenza rispetto alle alte. Dal punto di vista assoluto 0.6 e' molto vicino a zero ad  indicare che i dati sono domininati dal rumore bianco. 

L'esponente di Hurst viene calcolato come H=(β +1)/2 quindi 0.797 e la dimensione frattale D=2-β pari a 1.2 

La varazione dei parametri frattali e' compresa  

0<H<1

1<D<2

-1<β<3

Un ulteriore approccio e' quello di utilizzare l'analisi wavelet

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import pywt

from scipy.stats import linregress

# ======================================================

# 1 — LETTURA FILE CSV E COLONNA DI INTERESSE

# ======================================================

filename = "well_log.csv"

df = pd.read_csv(filename,nrows=2500)

# Utilizzo la colonna GAMMA (indice 6) come nei codici precedenti

x = df.iloc[:, 6].values # Dati di log (GAMMA)

# ======================================================

# 2 — PULIZIA DATI E RISCALAMENTO

# ======================================================

# Converti in numerico e gestisci NaN

x = pd.to_numeric(x, errors='coerce')

# Rimuovi i valori NaN.

x = x[~np.isnan(x)]

# Riscalamento se ci sono negativi

xmin = np.min(x)

if xmin <= 0:

shift = abs(xmin) + 1e-6

x = x + shift

print(f"Riscalamento effettuato: shift = {shift:.6f}")

# ======================================================

# 3 — PARAMETRI WAVELET

# ======================================================

wavelet = "cmor1.5-1.0"

scales = np.arange(1, 100)

# ======================================================

# 4 — CALCOLO CWT e POTENZA

# ======================================================

coeffs, frequencies = pywt.cwt(x, scales, wavelet, sampling_period=1)

power = np.abs(coeffs)**2

# ======================================================

# 5 — CALCOLO DELL'ESPONENTE DI HURST (H)

# ======================================================

# Calcola la Potenza Media su ciascuna scala (Varianza Wavelet)

scale_mean_power = np.mean(power, axis=1)

# Trasformazione log-log

log_scales = np.log10(scales)

log_mean_power = np.log10(scale_mean_power)

# Selezione della regione di scaling (dall'10% all'80% delle scale)

min_scale_index = int(len(scales) * 0.1)

max_scale_index = int(len(scales) * 0.8)

H_log_scales = log_scales[min_scale_index:max_scale_index]

H_log_mean_power = log_mean_power[min_scale_index:max_scale_index]

# Regressione Lineare

slope, intercept, r_value, p_value, std_err = linregress(H_log_scales, H_log_mean_power)

# Calcolo di H

hurst_exponent = (slope + 1) / 2

fractal_dimension = 2 - hurst_exponent

print("-" * 50)

print("Risultati dell'Analisi di Frattalità Wavelet:")

print(f"Esponente di Hurst (H): {hurst_exponent:.4f}")

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

print(f"Pendenza della regressione (2H - 1): {slope:.4f}")

print("-" * 50)

# ======================================================

# 6 — PLOT DEI RISULTATI (ALLINEAMENTO E RIMOZIONE COLORBAR)

# ======================================================

# NOTA: Usiamo sharex=True per allineare A e B. C (la regressione) deve essere separato.

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8),

sharex=True, gridspec_kw={'height_ratios': [3, 1]})

# --- A) SCALOGRAMMA WAVELET (SENZA SCALA DI POTENZA)

extent = [0, len(x), 0, len(scales)]

im = ax1.imshow(power, extent=extent, cmap="jet", aspect="auto", origin="lower")

ax1.set_title(f"Scalogramma Wavelet (Potenza per Scala)")

ax1.set_ylabel("Scala Wavelet")

# Rimozione della colorbar (scala di potenza)

# plt.colorbar(im, ax=ax1, orientation='vertical', label='Potenza') <-- RIMOSSA

# --- B) PLOT XY SOTTO (SEGNALE ORIGINALE)

ax2.plot(x, color="black")

ax2.set_xlabel("Indice campione (profondità)")

ax2.set_ylabel("Valore (GAMMA)")

ax2.grid(True)

plt.tight_layout()

plt.show()

# --- C) REGRESSIONE PER H (PLOT SEPARATO E NON ALLINEATO)

# Questo grafico non ha l'asse X della profondità, quindi deve rimanere separato

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

plt.loglog(scales, scale_mean_power, 'o', markersize=3, label='Potenza Media')

plt.loglog(scales[min_scale_index:max_scale_index],

intercept + slope * H_log_scales,

'r--', label=f'Regressione: pendenza={slope:.2f}')

plt.xlabel('Scala Wavelet (log)')

plt.ylabel('Potenza Media (log)')

plt.title(f'Calcolo dell\'Esponente di Hurst (H={hurst_exponent:.4f})')

plt.legend()

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

plt.show()