AutoKeras LSTM per estensimetro

Questo e' un tentativo di applicazione di AutoKeras (AutoML ovvero machine learning automatizzato basato su Keras) ai dati del post precedente

La differenza sostanziale e' che la dimensione del train set deve essere un multiplo intero della batch size  (in questo caso 42 e 4200)

# -*- coding: utf-8 -*-

"""Autokeras timeseries_forecaster

Automatically generated by Colaboratory.

Original file is located at

https://colab.research.google.com/drive/1mQHWTqwyRKhShtTrwNXBvPHgOkU4ghsX

"""

!pip install autokeras

import pandas as pd

import tensorflow as tf

import numpy as np

import matplotlib.pyplot as plt

import autokeras as ak

!rm detrend2.csv

!wget http://c1p81.altervista.org/detrend2.zip

!unzip detrend2.zip

dataset=pd.read_csv(r'detrend2.csv', sep=',')

# cancella la prima colonna del tempo

#dataset_c.drop('Data',axis=1,inplace=True)

dataset= dataset_c[1:6001]

dataset.head()

dataset.shape

val_split = int(len(dataset) * 0.7)

data_train = dataset[:val_split]

validation_data = dataset[val_split:]

data_x = data_train[

[

"Est.[mm]"

]

].astype("float64")

data_x_val = validation_data[

[

"Est.[mm]"

]

].astype("float64")

# Data with train data and the unseen data from subsequent time steps.

data_x_test = dataset[

[

"Est.[mm]"

]

].astype("float64")

data_y = data_train["Est.[mm]"].astype("float64")

data_y_val = validation_data["Est.[mm]"].astype("float64")

##print(data_x.shape)

#print(data_y.shape)

#print(data_y)

#la batch size deve essere un divisore

#della serie di train altrimenti genera errore

predict_from = 1

#lunghezza dei dati nel futuro

predict_until = 42

lookback = 30

clf = ak.TimeseriesForecaster(

lookback=lookback,

predict_from=predict_from,

predict_until=predict_until,

max_trials=4,

objective="val_loss",

)

clf.fit(

x=data_x,

y=data_y,

validation_data=(data_x_val, data_y_val),

batch_size=42,

epochs=20

)

# Predizione

predictions = clf.predict(data_x_test)

print(predictions.shape)

print(data_x_val.shape)

print(data_y_val.shape)

loss,acc = clf.evaluate(data_x_val, data_y_val, batch_size=42, verbose=1)

print('Accuracy: %.3f' % acc)

plt.plot(predictions)

plt.title("Estensimetro")

plt.xlabel('')

plt.ylabel('')

plt.show()

model = clf.export_model()

# summarize the loaded model

model.summary()

# save the best performing model to file

model.save('model_sonar.h5')

LSTM Time Series con Tensorflow per la geologia

Questo post e' un tentativo di applicazione delle reti neurali LSTM per l'analisi di dati (sulla base di questo articolo) derivanti da un sensore ad applicazione geologica, nello specifico un estensimetro. I dati sono dati reali che sono stati anonimizzati ma sono relativi ad un movimento di versante attivo. Sono originali sono acquisiti con passo di 20 minuti e coprono circa 10 mesi e prevedono misure anche meteo

Il primo tentativo e' stato quello di esaminare la curva completa ma presto mi sono reso conto che la rete neurale non riusciva mai a convergere. Ho limitato i dati quindi dall'inizio fino alla primo movimento avvenuto l'8 maggio

Come si osserva dal grafico ridotto sottostante oltre ad un trend generale in crescita dei valori dell'estensimetro vi sono variazioni cicliche a scala corta

in generale le variazioni cicliche a breve periodo di un estensimetro sono relative a dilatazioni termiche

Plottando i dati si ha comunque una scarsa correlazione

Un approccio di analisi base prevede di trovare una curva che approssima i dati (in questo caso la migliore approssimazione e' stata una polinomio di secondo grado)

facendo uno zoom si ha vede che le variazioni a breve periodo non sono sinusoidali (in pratica si osserva solo dei picchi di discesa ) e che approssimando con la curva di interpolazione si ha una fascia di incertezza dell'ordine di 1 mm

Vediamo se analizzando la serie tempo con una rete LSTM in Tensorflow si riesce ad ottenere un modello che fitta meglio con i dati

Per prima cosa devo dire che i primi tentativi sono stati piuttosto deludenti..mi sono trovato spesso in situazioni in cui la loss dei dati di training era ottima mentre era pessima  (e spesso in crescita con le epochs) quella del set di validazione. 

In altri casi mi sono trovato ad avere un modello decente che aveva un offset costante rispetto ai dati di controllo (vedi grafico sottostante)

Una soluzione e' stata quella di inserire dei layer di dropout nella rete e quella di effettuare un detrending (sottrarre ai dati grezzi il valore della funzione che approssima la tendenza generale)

Un altro aspetto determinante per un modello ottimale consiste nell'individuazione della corretta batch size (nel caso specifico batch size inferiore a 50 non riusciva a seguire le variazioni dei dati)

un dettaglio dei dati dopo il detrend

il miglior risultato che sono riuscito ad ottenere e' il seguente

questi sono i dati di train contro dati reali. Si nota che il modello segue l'andamento dei dati reali ma mostra un ritardo nei minimi relativi del trend generale

 

questi sono invece i dati di validazione. Si osserva una buona sincronia sulle variazioni a breve periodo

========================================

# -*- coding: utf-8 -*-

"""Quincineto_stima_errore_detrend2.ipynb

Automatically generated by Colaboratory.

Original file is located at

https://colab.research.google.com/drive/1fFEfuUj5iwNucKDOfVpEzvIYK_1KebH8

### Elaborazione Quincineto

"""

import pandas as pd

import numpy as np

from sklearn.preprocessing import MinMaxScaler

import matplotlib.pyplot as plt

import tensorflow as tf

import os

# Commented out IPython magic to ensure Python compatibility.

#prepara tensorboard

# %load_ext tensorboard

"""

Download dati"""

!wget http://c1p81.altervista.org/detrend2.zip

!unzip detrend2.zip

"""Carica i dati in matrice"""

df=pd.read_csv(r'detrend2.csv', sep=',', header=0,index_col=['Data'])

df.head()

plt.plot(df['Est.[mm]'])

plt.title("Estensimetro")

plt.xlabel('Datetime')

plt.ylabel('Est')

plt.show()

#finestra mobile dei dati

n_past = 300

n_future = 100

n_features = 1

# divide il dataset in 75% train, 25 % test

righe = df.shape[0]

t = np.round(righe*0.75,0)

print(t)

train_df,test_df = df[1:6000], df[6000:]

# il dataset di test non entra nel calcolo della rete neurale

# e' lo stato futuro che la rete deve prevedere

# qui si trova la rottura della tendenza

plt.plot(test_df['Est.[mm]'])

plt.title("Train dataset estensimetro Dettaglio")

plt.xlabel('Datetime')

plt.ylabel('Est')

plt.show()

#riscala i dati per il calcolo della rete

train = train_df

scalers={}

for i in train_df.columns:

scaler = MinMaxScaler(feature_range=(-1,1))

s_s = scaler.fit_transform(train[i].values.reshape(-1,1))

s_s=np.reshape(s_s,len(s_s))

scalers['scaler_'+ i] = scaler

train[i]=s_s

test = test_df

for i in train_df.columns:

scaler = scalers['scaler_'+i]

s_s = scaler.transform(test[i].values.reshape(-1,1))

s_s=np.reshape(s_s,len(s_s))

scalers['scaler_'+i] = scaler

test[i]=s_s

"""**Converting the series to samples for supervised learning**"""

def split_series(series, n_past, n_future):

#

# n_past ==> no of past observations

#

# n_future ==> no of future observations

#

X, y = list(), list()

for window_start in range(len(series)):

past_end = window_start + n_past

future_end = past_end + n_future

if future_end > len(series):

break

# slicing the past and future parts of the window

past, future = series[window_start:past_end, :], series[past_end:future_end, :]

X.append(past)

y.append(future)

return np.array(X), np.array(y)

X_train, y_train = split_series(train.values,n_past, n_future)

X_train = X_train.reshape((X_train.shape[0], X_train.shape[1],n_features))

y_train = y_train.reshape((y_train.shape[0], y_train.shape[1], n_features))

X_test, y_test = split_series(test.values,n_past, n_future)

X_test = X_test.reshape((X_test.shape[0], X_test.shape[1],n_features))

y_test = y_test.reshape((y_test.shape[0], y_test.shape[1], n_features))

X_test.shape

print(train)

# E1D1

# n_features ==> no of features at each timestep in the data.

#

encoder_inputs = tf.keras.layers.Input(shape=(n_past, n_features))

encoder_l1 = tf.keras.layers.LSTM(300, return_state=True)

encoder_outputs1 = encoder_l1(encoder_inputs)

encoder_states1 = encoder_outputs1[1:]

#

decoder_inputs = tf.keras.layers.RepeatVector(n_future)(encoder_outputs1[0])

#

decoder_l1 = tf.keras.layers.LSTM(300, return_sequences=True,dropout=0.5,recurrent_dropout=0.5)(decoder_inputs,initial_state = encoder_states1)

decoder_outputs1 = tf.keras.layers.TimeDistributed(tf.keras.layers.Dense(n_features))(decoder_l1)

#

model_e1d1 = tf.keras.models.Model(encoder_inputs,decoder_outputs1)

model_e1d1.summary()

reduce_lr = tf.keras.callbacks.LearningRateScheduler(lambda x: 1e-3 * 0.90 ** x)

epoche = 15

#tensorboard start

import datetime

log_dir = "Logs/fit/" + datetime.datetime.now().strftime("%Y%m%d-%H%M%S")

tensorboard_callback = tf.keras.callbacks.TensorBoard(log_dir=log_dir, histogram_freq=1)

#tensorboard start

model_e1d1.compile(optimizer=tf.keras.optimizers.Adam(), loss=tf.keras.losses.Huber())

history_e1d1=model_e1d1.fit(X_train,y_train,epochs=epoche,validation_data=(X_test,y_test),batch_size=164,verbose=1,callbacks=[reduce_lr,tensorboard_callback])

plt.plot(history_e1d1.history['loss'])

plt.plot(history_e1d1.history['val_loss'])

plt.title("Model Loss")

plt.xlabel('Epochs')

plt.ylabel('Loss')

plt.legend(['Train', 'Valid'])

plt.show()

plt.plot(history_e1d1.history['lr'])

plt.title("Model Lr")

plt.xlabel('Epochs')

plt.ylabel('Lr')

plt.show()

"""Predizione dei dati"""

pred1_e1d1=model_e1d1.predict(X_test)

pred_e1d1=model_e1d1.predict(X_train)

"""Ritorna dai valori scalati ai valori reali"""

for index,i in enumerate(train_df.columns):

scaler = scalers['scaler_'+i]

pred1_e1d1[:,:,index]=scaler.inverse_transform(pred1_e1d1[:,:,index])

pred_e1d1[:,:,index]=scaler.inverse_transform(pred_e1d1[:,:,index])

y_train[:,:,index]=scaler.inverse_transform(y_train[:,:,index])

y_test[:,:,index]=scaler.inverse_transform(y_test[:,:,index])

#print(y_train)

np.savetxt('array_ytrain_X.csv', y_train[:,:,0], delimiter=';', fmt='%f')

np.savetxt('array_ytrain_Y.csv', y_train[:,:,1], delimiter=';', fmt='%f')

np.savetxt('array_ytrain_Est.csv', y_train[:,:,2], delimiter=';', fmt='%f')

np.savetxt('array_ytrain_T.csv', y_train[:,:,3], delimiter=';', fmt='%f')

np.savetxt('array_ytrain_bat.csv', y_train[:,:,4], delimiter=';', fmt='%f')

np.savetxt('array_ytest_X.csv', y_test[:,:,0], delimiter=';', fmt='%f')

np.savetxt('array_ytest_Y.csv', y_test[:,:,1], delimiter=';', fmt='%f')

np.savetxt('array_ytest_Est.csv', y_test[:,:,2], delimiter=';', fmt='%f')

np.savetxt('array_ytest_T.csv', y_test[:,:,3], delimiter=';', fmt='%f')

np.savetxt('array_ytest_bat.csv', y_test[:,:,4], delimiter=';', fmt='%f')

print(pred1_e1d1[:,:,0])

"""**Checking Error** """

#errore medio tra modello e dati reali

from sklearn.metrics import mean_absolute_error

print(y_test.shape)

#print(pred1_e1d1[:,:,2])

#print(y_test[:,:,2])

#print(y_test.shape[0])

diff = pred1_e1d1[:,:,0]-y_test[:,:,0]

print("Diff")

#print(diff)

print(diff.shape)

#print(np.square(diff))

sommaq=np.sum(np.square(diff[0]))

print(sommaq.shape)

s_err1=np.sqrt(sommaq/y_test.shape[0]*y_test.shape[1])

print()

print(sommaq)

#err1 = np.square(pred1_e1d1[:,:,2]-y_test[:,:,2])

#s_err1 = np.sqrt(np.sum(err1, axis=0))

#print("Errore quadratico medio (mm)")

print(s_err1)

#print(y_test)

plt.plot(y_test[:,99,0])

plt.plot(pred1_e1d1[:,99,0])

plt.title("Modello vs monitoraggio")

plt.xlabel('Tempo')

plt.ylabel('Est')

plt.legend(['Dati reali', 'Modello'])

plt.show()

# Modello contro monitoraggio nei dati di addestramento

plt.plot(y_train[:,99,2])

plt.plot(pred_e1d1[:,99,2])

plt.title("")

plt.xlabel('Tempo')

plt.ylabel('Est (mm)')

plt.legend(['Monitoraggio', 'Modello'])

plt.show()

# dettaglio del grafico Modello contro dati reali

# si nota come il modello copia in modo idoneo anche

# le variazioni a corto periodo

plt.plot(y_train[1:1000,99,2])

plt.plot(pred_e1d1[1:1000,99,2])

plt.title("")

plt.xlabel('Tempo')

plt.ylabel('Est (mm)')

plt.legend(['Monitoraggio', 'Modello'])

plt.show()

# Commented out IPython magic to ensure Python compatibility.

# %tensorboard --logdir Logs/fit

Utilizzando l'esempio del Time Series Forecasting direttamente dal sito di Tensorflow si hanno risultati similari

Soliflusso

Un esempio di soliflusso nei prati della Consuma. 

Circuito elettrico per freatimetro

 Mi e' capitato di trovare in ufficio un freatimetro abbandonato in un armadio e, visto che non funzionava, ho provato ad aprirlo

Avevo sempre immaginato che il circuito di un freatimetro fosse una semplice batteria ed un buzzer.

Con sorpresa ho trovato all'interno un integrato LM393N, un comparatore di tensione 

Questa immagine e' stata ribaltata sull'asse verticale in modo da poter avere la corrispondenza con l'immagine soprastante

Controllando le piste ci sono un paio di componenti (CIC e L1) assenti nonostante le saldature non siano presenti...non sara' cosi' facile rimetterlo in funzione.. oltre all'ovvio cavo negativo della batteria che si e' dissaldato

Iphone 12 Pro per la geologia

Avevo provato tempo fa ad usare il tablet Google Tango per scansionare un affioramento roccioso

Adesso che anche Iphone 12 pro (e Ipad Pro 2020) ha un lidar volevo vedere la possibilita' di usare il telefono di Apple per applicazione geologiche

 

Mettendo a confronto Tango con Iphone si nota che entrambi hanno un raggio di scansione massimo tra i 3 ed i 4 metri ma sicuramente la IMU di Iphone e' migliore in quanto spostandosi non decorrela in modo significativo anche su lunghezze superiori ai 10 m

Scansione con Iphone

Scansione con Google Tango Tablet

La cosa che mi impressiona sempre e' come sia possibile ottenere delle misure di lunghezza senza la necessita' di ricalibrare i dati ma questo aspetto e' comune sia a Tango che ad IPhone

Misura di strato : a sinistra di misura reale, a destra misura dal mesh

Non esiste una applicazione ufficiale per esportare i cloudpoint o la mesh da IPhone ma sono disponibili 3DScanner App, PolyCam e Heges

I dati comunque sono leggibili ed elaboraibili da CloudCompare 

PCL e Tango tablet

Ho ritirato fuori il mio tablet Google Yellowstone Tango per provare il trattamento dati delle nuvole dei punti con PCL

I dati reali sono stati presi presso la cava di Maiano (Fiesole) che e' la palestra di geologia per numerose generazioni di geologi fiorentini 

Ho ripreso con lo scanner del tablet questo dettaglio 

in particolare volevo vedere se riuscivo a misurare l'angolo tra i due piani indicati nella foto sottostante

(quello in blu e' il piano del fronte di scavo, quello in rosso e' relativo ad una frattura)

Usando la app Clino Fieldmove e Innstereo 

a posteriori la scelta delle superfici non e' stata felicissima perche' sono tutte ad alto angolo e si distinguono poco sullo stereoplot

Con il tablet la distanza massima a cui era possibile avere risposta dallo scanner era di circa 3 m. Per prova la superficie e' stata bagnata per vedere se il laser infrarosso era influenzato dall'umidita' della parete ma non si verifiche significative tra superficie asciutta e bagnata

I dati sono stati salvati in PLY e sono stati elaborati con PCL in QtCreator

Il programma dopo aver caricato il file PLY fa un sottocampionamento dei dati con VoxelGrid (passando da oltre 200.000 superfici a poco oltre 50 superfici) e calcola le  normali. I dati sono visualizzati con il visualizzatore interno a PLC

leggermente modificato da https://github.com/jeffdelmerico/pointcloud_tutorial

========================================

cmake_minimum_required(VERSION 3.5)

project(maiano LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 11)

set(CMAKE_CXX_STANDARD_REQUIRED ON)

find_package(PCL 1.3 REQUIRED COMPONENTS common io features visualization)

include_directories(${PCL_INCLUDE_DIRS})

link_directories(${PCL_LIBRARY_DIRS})

add_definitions(${PCL_DEFINITIONS})

add_executable(maiano main.cpp)

target_link_libraries(maiano ${PCL_LIBRARIES})

========================================

#include <iostream>

#include <pcl/point_types.h>

#include <pcl/point_cloud.h>

#include <pcl/io/ply_io.h>

#include <pcl/io/pcd_io.h>

#include <pcl/filters/voxel_grid.h>

#include <pcl/kdtree/kdtree_flann.h>

#include <pcl/features/normal_3d.h>

#include <pcl/visualization/pcl_visualizer.h>

using namespace std;

void

downsample (pcl::PointCloud<pcl::PointXYZRGB>::Ptr &points, float leaf_size,

            pcl::PointCloud<pcl::PointXYZRGB>::Ptr &downsampled_out)

{

    cout << "PointCloud before filtering: " << points->width * points->height

           << " data points (" << pcl::getFieldsList (*points) << ")." << std::endl;

  pcl::VoxelGrid<pcl::PointXYZRGB> vox_grid;

  vox_grid.setLeafSize (leaf_size, leaf_size, leaf_size);

  vox_grid.setInputCloud (points);

  vox_grid.filter (*downsampled_out);

  cout << "PointCloud after filtering: " << downsampled_out->width * downsampled_out->height

         << " data points (" << pcl::getFieldsList (*downsampled_out) << ")." << std::endl;

}

void compute_surface_normals (pcl::PointCloud<pcl::PointXYZRGB>::Ptr &points, float normal_radius,

                                    pcl::PointCloud<pcl::Normal>::Ptr &normals_out)

{

  pcl::NormalEstimation<pcl::PointXYZRGB, pcl::Normal> norm_est;

  // Use a FLANN-based KdTree to perform neighborhood searches

  norm_est.setSearchMethod (pcl::search::KdTree<pcl::PointXYZRGB>::Ptr

                            (new pcl::search::KdTree<pcl::PointXYZRGB>));

  // Specify the size of the local neighborhood to use when computing the surface normals

  norm_est.setRadiusSearch (normal_radius);

  // Set the input points

  norm_est.setInputCloud (points);

  // Estimate the surface normals and store the result in "normals_out"

  norm_est.compute (*normals_out);

}

void visualize_normals (const pcl::PointCloud<pcl::PointXYZRGB>::Ptr points,

                        const pcl::PointCloud<pcl::PointXYZRGB>::Ptr normal_points,

                        const pcl::PointCloud<pcl::Normal>::Ptr normals)

{

  // Add the points and normals to the vizualizer

  pcl::visualization::PCLVisualizer viz;

  viz.addPointCloud (points, "points");

  viz.addPointCloud (normal_points, "normal_points");

  viz.addPointCloudNormals<pcl::PointXYZRGB, pcl::Normal> (normal_points, normals, 1, 0.01, "normals");

  // Give control over to the visualizer

  viz.spin ();

}

int main (int argc, char** argv)

{

    // Load data from pcd

    pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZRGB>);

    //if (pcl::io::loadPCDFile<pcl::PointXYZRGB> ("../data/robot1.pcd", *cloud) == -1) //* load the file

    if (pcl::io::loadPLYFile<pcl::PointXYZRGB> ("maiano.ply", *cloud) == -1) //* load the file

    {

        PCL_ERROR ("Couldn't read file robot1.pcd \n");

        return (-1);

    }

    // Point Clouds to hold output

    pcl::PointCloud<pcl::PointXYZRGB>::Ptr downsampled (new pcl::PointCloud<pcl::PointXYZRGB>);

    pcl::PointCloud<pcl::Normal>::Ptr normals (new pcl::PointCloud<pcl::Normal>);

    // Downsample the cloud

    const float voxel_grid_leaf_size = 1.0;

    downsample (cloud, voxel_grid_leaf_size, downsampled);

    // Compute surface normals

    const float normal_radius = 1.0;

    compute_surface_normals (downsampled, normal_radius, normals);

    visualize_normals(cloud, downsampled, normals);

    int i = 0;

    for (pcl::Normal n : *normals) {

        std::cout << n.normal_x << "," << n.normal_y << "," << n.normal_z <<"\n";

        i++;

    }

    return(0);

}

========================================

Plottando le normali si ha questo risultato

Al di la' dell'orientamento verso il Nord (che lo scanner non registra in quanto non calibra i dati con la bussola) il programma estra sia le suprfici ad alto angolo che alcun suborizzontali.. non so quanto i dati alla fine siano corretti....i dati delle normali vengono forniti come componenti x,y,z di un versore con origine 0,0,0....non so se ho fatto qualche sbaglio nella conversione in strike e dip