Animazione di grafico con Matplotlib

 

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import pandas as pd
df = pd.read_csv("tempo.csv", dtype=str,sep=";")
print(df)
my_data = np.genfromtxt('./scatterplot.csv', delimiter=';')
x = my_data[:,1]
y = my_data[:,0]
fig = plt.figure()
ax = fig.add_subplot()
plt.ylim(-5, 22)
plt.ylabel ('Temperatura (C)')
plt.xlim(-59, -53)
plt.xlabel ('Estensimetro (mm)')
tempo = ax.text(0.8, 0.9, '', horizontalalignment='center',verticalalignment='center', transform=ax.transAxes)
graph, = plt.plot([], [], 'o', markersize=1)
def animate(i):
    graph.set_data(x[:i+1], y[:i+1])
    tempo.set_text(df.Data[i])
    return graph,
ani = FuncAnimation(fig, animate, frames=2000, interval=50)
#FFwriter = animation.FFMpegWriter(fps=10)
#ani.save('animation.mp4')
plt.show()

Calcolo di angolo tra vettori con NumPy

Calcolo di angolo tra due vettori con numpy secondo la formuale sottostante (prodotto scalare tra i due vettori diviso il prodotto delle norme dei vettori)

https://www.wikihow.it/Calcolare-l%27Angolo-tra-Due-Vettori

import numpy as np
A = np.array([3,3])
B = np.array([2,1])
dot_product = np.dot(A,B)
normA = (np.linalg.norm(A))
normB = (np.linalg.norm(B))
cos_angolo = dot_product/(normA*normB)
angolo_rad = np.arccos(cos_angolo)
print(np.rad2deg(angolo_rad))

Coregistrazione immagini ottiche

 Questo post illustra come e' stato risolto il problema della GIF sottostante

Le immagini ottiche sono riprese da una camera PTZ che nel corso del giorno si sposta su diverse inquadrature...cio' comporta che non ritorna mai nelle impostazione perfette del giorno precedente rendendo impossibile l'analisi al computer ma anche la semplice interpretazione umana di eventuali modifiche

Per risolvere il problema sono state applicate tecniche che riguardano sia la correzione geometrica che dello spazio colore

Per prima cosa i vari fotogrammi sono stati coregistrati in modo software con il seguente script Python

///////////////////////////////////////////////////////

import cv2

import numpy as np

import sys, getopt

file2 = sys.argv[1:][0]

#print(file2)

#print("luca")

#print("/home/luca/ptz/150.217.73.10/share/transi"+file2)

# Read the images to be aligned

im1 = cv2.imread("/home/luca/ptz/150/share/transi/1.jpg");

im2 = cv2.imread("/home/luca/ptz/150/share/transi/"+file2);

# Convert images to grayscale

im1_gray = cv2.cvtColor(im1,cv2.COLOR_BGR2GRAY)

im2_gray = cv2.cvtColor(im2,cv2.COLOR_BGR2GRAY)

# Find size of image1

sz = im1.shape

# Define the motion model

warp_mode = cv2.MOTION_TRANSLATION

# Define 2x3 or 3x3 matrices and initialize the matrix to identity

if warp_mode == cv2.MOTION_HOMOGRAPHY :

warp_matrix = np.eye(3, 3, dtype=np.float32)

else :

warp_matrix = np.eye(2, 3, dtype=np.float32)

# Specify the number of iterations.

number_of_iterations = 5000;

# Specify the threshold of the increment

# in the correlation coefficient between two iterations

termination_eps = 1e-10;

# Define termination criteria

criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations, termination_eps)

# Run the ECC algorithm. The results are stored in warp_matrix.

(cc, warp_matrix) = cv2.findTransformECC (im1_gray,im2_gray,warp_matrix, warp_mode, criteria)

if warp_mode == cv2.MOTION_HOMOGRAPHY :

# Use warpPerspective for Homography

im2_aligned = cv2.warpPerspective (im2, warp_matrix, (sz[1],sz[0]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)

else :

# Use warpAffine for Translation, Euclidean and Affine

im2_aligned = cv2.warpAffine(im2, warp_matrix, (sz[1],sz[0]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP);

# Show final results

#cv2.imshow("Image 1", im1)

#cv2.imshow("Image 2", im2)

#cv2.imshow("Aligned Image 2", im2_aligned)

cv2.imwrite("/home/luca/ptz/150/share/transi/all_"+file2,im2_aligned)

print("allinea")

print(file2)

///////////////////////////////////////////////////////

per la registrazione dell'istogramma dei colori e' stato usato

///////////////////////////////////////////////////////

from skimage import exposure

import matplotlib.pyplot as plt

import argparse

import cv2

# construct the argument parser and parse the arguments

ap = argparse.ArgumentParser()

ap.add_argument("-s", "--source", required=True,help="numero")

args = vars(ap.parse_args())

src = cv2.imread("/home/luca/ptz/150/share/transi/all_"+str(args["source"])+".jpg")

ref = cv2.imread("/home/luca/ptz/150/share/transi/all_1.jpg")

# determine if we are performing multichannel histogram matching

# and then perform histogram matching itself

print("[INFO] performing histogram matching...")

multi = True if src.shape[-1] > 1 else False

matched = exposure.match_histograms(src, ref, multichannel=multi)

cv2.imwrite("/home/luca/ptz/150/share/transi/orfeo"+str(args["source"])+".jpg",matched)

///////////////////////////////////////////////////////

per finire le immagini sono passate all'algoritmo di Change Detection di Orfeo ToolBox

otbcli_MultivariateAlterationDetector -in1 ~/ptz/150/share/transi/orfeo$j.jpg -in2 /home/luca/ptz/150.217.73.10/share/transi/orfeo$prec.jpg -out /home/luca/ptz/150/share/transi/change$j.jpg

ed il risultato viene elaborato da uno script in Python per esaltare il contrasto della mappa di CD

import cv2

import numpy as np

import sys, getopt

file2 = sys.argv[1:][0]

img = cv2.imread("/home/luca/ptz/150/transi/change"+file2+".jpg", cv2.IMREAD_COLOR)

# normalize float versions

#norm_img1 = cv2.normalize(img, None, alpha=0, beta=1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)

norm_img2 = cv2.normalize(img, None, alpha=0, beta=1.2, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)

# scale to uint8

#norm_img1 = (255*norm_img1).astype(np.uint8)

norm_img2 = np.clip(norm_img2, 0, 1)

norm_img2 = (255*norm_img2).astype(np.uint8)

# write normalized output images

cv2.imwrite("/home/luca/ptz/150/transi/change_stretch"+file2+".jpg",norm_img2)

#cv2.imwrite("zelda1_bm20_cm20_normalize2.jpg",norm_img2)

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:])

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)

Mandelbrot M1 Metal Python

 Dopo il test di Cuda e' arrivato il momento di Metal, il linguaggio per la GPU di Apple

Gli ersempi in C++ sono un po' complicati per una semplice prova ed ho optato per la versione in Python in cui il kernel e' comunque scritto in C++ MetalCompute (negli esempi della libreria c'e' gia' una implementazione di Mandelbrot, l'ho scritta da zero per imparare...sicuramente la mia versione e' meno efficiente)

La logica seguita e' quella di Cuda, ogni thread si occupa di una riga

Si puo' solo passare un solo parametro come array di float al kernel e si ha in uscita un solo array di float. L'array di input e' generato con NumPy con le coordinate X,Y dei punti dell'insieme 

import metalcompute as mc

from PIL import Image as im

mc.init()

code = """

#include <metal_stdlib>

using namespace metal;

kernel void mandelbrot(

const device float* arr [[ buffer(0) ]],

device float* out [[ buffer(1) ]],

uint id [[ thread_position_in_grid ]]) {

float a = arr[id];

float b;

float x_new,y_new,x,y;

int iterazioni = 255;

for (int s=0; s<640;s++)

{

x=0.0;

y=0.0;

b = arr[640+s];

for (int k=0;k<iterazioni;k++)

{

x_new = (x*x)-(y*y)+a;

y_new = (2.0*x*y)+b;

if (((x_new*x_new)+(y_new*y_new))>4)

{

out[id*640+s] = (float)(k%2)*255;

k=iterazioni;

}

x = x_new;

y = y_new;

}

}

}

"""

mc.compile(code, "mandelbrot")

import numpy as np

from time import time as now

dimensioni = 640

x_start = -2.0

x_stop = 0.75

x_res = (x_stop-x_start)/dimensioni

x = np.arange(start=x_start, stop=x_stop,step = x_res,dtype='f')

y_start = -1.5

y_stop = 1.5

y_res = (y_stop-y_start)/dimensioni

y = np.arange(start=y_start, stop=y_stop,step = y_res,dtype='f')

arr = np.concatenate((x,y))

print(arr.shape)

out_buf = np.empty(dimensioni*dimensioni,dtype='f')

start = now()

mc.run(arr, out_buf, dimensioni)

end = now()

immagine = out_buf.reshape(dimensioni,dimensioni)

print(end-start)

immagine = np.int_(immagine)

np.savetxt("dati.txt",immagine.astype(int),fmt='%i')

img = im.fromarray(np.uint8(immagine) , 'L')

img.save('m1.png')

img.show()

Preparare i dati per Tensorflow

Questo spero si tratti l'inizio di alcuni post sulla possibilita' di usare Tensorflow per stimare la concentrazione di alcune specie chimiche partendo da spettri IR

Per questo scopo c'e' bisogno di un database in cui sono presenti spettri IR ed analisi composizionali del materiale. L'unico database pubblico che conosco e' usgs_splib07 da questo link

Per semplicita' sono stati selezionati gli spettri a 480 bande

		SiO2	Al2O3	MgO
1	Actinolite HS116	57,8	0,22	22,4
2	Actinolite HS22	54,5	1,64	21,69
3	Actinolite HS315	57,7	1,38	24,58
4	Albite NMNHC5390	68,18	20,07	0,02
5	Almandine HS114	40,4	22,97	0,53
6	Almandine WS475	37,98	22,14	7,35
7	Almandine WS476	39,84	22,04	10,69
8	Almandine WS478	39,09	22,05	5,74
9	Almandine WS479	37,34	21,93	8,83
10	Andradite NMNH113829	35,5	0,29	0,14
11	Andradite WS474	35,85	0,085	0,67
12	Annite WS661	33,48	14,51	0,006
13	Antigorite NMNH96917a	41,3	1,59	36,3
14	Augite NMNH120049	48,93	8,44	14,37
15	Biotite WS660	39,52	11,05	14,45
16	Chlorite SMR-13	31	17,3	30,2
17	Clinochlore GDS158	27,5	19,8	24,7
18	Clinochlore GDS159	31,3	19,4	34,7
19	Clinochlore NMNH83369	32	16	33,9
20	Diopside NMNHR18685	54,77	0,49	18,26
21	Elbaite NMNH94217	36,83	40,56	0,03
22	Epidote GDS26	37,14	23,05	0,06
23	Grossular NMNH155371	38,84	17,95	0,15
24	Grossular WS484	39,71	22,58	0,019
25	Halloysite CM13	46,3	35,2	0,31
26	Hectorite SHCa-1	34,7	0,69	15,3
27	Hydrogrossular NMNH120555	31,44	21,4	0,17
28	Hypersthene NMNHC2368	51,32	4,89	26,09
29	Illite GDS4 (Marblehead)	51,62	23,96	3,83
30	Illite IMt-1	52,1	21,9	2,39
31	Kaolinite CM5	44,4	37,9	0,11
32	Kaolinite CM7	44,6	36,4	0,14
33	Kaolinite CM9	47,1	37,4	0,16
34	Kaolinite KGa-1 (wxl)	45	38	0,02
35	Kaolinite KGa-2 (pxl)	44,2	37,2	0,04
36	Lizardite NMNHR4687	40,97	0	40,43
37	Marialite NMNH126018-2	59,88	20,59	0
38	Meionite WS701	49,6	25,45	0,06
39	Mizzonite NMNH113775	49,6	25,43	0,03
40	Montmorillonite CM20	52	18,36	3,32
41	Montmorillonite CM26	63,4	18,2	2,26
42	Montmorillonite CM27	57,6	20,9	2,66
43	Montmorillonite+Illite CM37	46,5	22,7	2,74
44	Montmorillonite+Illite CM42	54,8	21	3,4
45	Montmorillonite SAz-1	60,4	17,6	6,46
46	Montmorillonite SCa-2	52,4	15	6,68
47	Montmorillonite STx-1	69,6	16,3	3,56
48	Montmorillonite SWy-1	62,9	19,3	2,8
49	Muscovite GDS108	44	35	0,57
50	Nontronite NG-1	80	2	0,42
51	Nontronite SWa-1	74,6	7,69	3,32
52	Olivine GDS70 Fo89	41,09	0	49,29
53	Olivine GDS71 Fo91	40,06	0	50,7
54	Olivine KI3005 Fo11	30,11	0	4,42
55	Olivine KI3054 Fo66	36,3	0	32,62
56	Olivine KI3188 Fo51	34,34	0	23,8
57	Olivine KI3189 Fo60	35,47	0	29,49
58	Palygorskite CM46a	56	10,6	8,56
59	Pectolite NMNH94865	53,37	0,03	0,03
60	Phlogopite GDS20	40,3	14,3	26,4
61	Phlogopite WS496	43,015	11,6	29,62
62	Quartz HS117 Aventurine	100	0	0
63	Rhodonite NMNHC6148	46,3	0,05	0,41
64	Richterite NMNH150800	54,87	1,97	22,92
65	Saponite SapCa-1	51,9	3,66	0,88
66	Sepiolite SepNev-1.AcB	50,9	0,45	22,3
67	Talc HS21	60,47	0,12	30,29
68	Talc WS659	58,38	0,18	31,9
69	Thuringite SMR-15	24,7	21,1	16,8
70	Tremolite NMNH117611	54,6	1,88	23,13
71	Vermiculite VTx-1.a	19,5	1,53	27,7

I dati sono stati estratti dalla directory ASCIIData/ASCIIdata_splib07a/ChapterM_Minerals selezionando gli spettri denominati BECKb_AREF.
Un esempio di file e'

---------------------------------------
 splib07a Record=120: Actinolite HS116.3B           BECKb AREF
-1.2300000e+034
 3.3700544e-002
 3.5717472e-002
 3.5481069e-002
 3.5344660e-002
 3.4786128e-002
 3.3701397e-002
 3.3403333e-002
 3.3196956e-002
 3.3779550e-002
 3.2682978e-002
 3.3049587e-002
 3.4390874e-002
 3.4718834e-002
 3.6168687e-002
...............................................
---------------------------------------

Si iniziano quindi a preparare i dati per inserirli in un array compatibile con NumPy,formato idoneo al successivo utilizzo con Tensorflow

per prima cosa si toglie la prima riga di intestazione con sed

for f in *.txt ; do sed -i '1d' $f; done

viene quindi eseguito uno script python che legge i dati riga per riga e restituisce una stringa lineare con i dati separati da un punto e virgola

converti.py
=================================
with open("1.txt",'r') as f:
while True:
buf = f.readline()
if not buf:
break
print (buf.strip()),
print ';',
=================================

for f in *.txt ; do python converti.py > a$f; done

i nuovi file avranno il carattere a anteposto al precedente nome file

for f in *.txt ; do python converti.py > a$f; done

si tolgono quindi gli spazi bianchi

for f in a*.txt ; do cat $f  | tr - d "[:space:]" > b$f ; done

i nuovi file avranno il carattere b anteposto al precedente nome file

siamo quasi arrivati alla fine. Ogni file dello spettro e' costituito da un riga, si deve quindi creare un unico file in cui ogni riga rappresenta uno spettro

for f in ba*.txt; do (cat "${f}"; echo) >> completo.txt; done
from numpy import genfromtxt
data = genfromtxt('completo.txt', delimiter=',')

Da notare bene che tensorflow non funziona su macchine con processori non recenti perche' vengono usate delle istruzioni cpu recenti (AVX and AVX2 instruction sets)

Per installare Tensorflow su Debian con virtualenv si usa la sequenza (ripresa da questo sito)

si installa python3

apt-get install python3-venv python3 

mkdir spectra

cd spectra

python3 -m venv venv

source venv/bin/activate

pip install --upgrade pip

pip install --upgrade tensorflow

per verificare l'installazione si usa 

python -c 'import tensorflow as tf; print(tf.__version__)'

per terminare virtualenv si digita

deactivate

A questo punto ho diviso il file di 71 spettri in una serie di train di 60 spettri ed una serie di test di 11 spettri (x_train,x_test) . Per quanto riguarda le concentrazioni il file delle concentrazioni (y_train, y_test) i valori di concentrazioni sono stati divisi in 10 classi con la seguente regola

0-9.9% : classe 0
10-19.9% : classe  1
20-29.9% : classe  2
30-39.9% : classe  3
40-49.9% : classe  4
50-59.9% : classe  5
60-69.9% : classe  6
70-79.9% : classe  7
80-89.9% : classe  8
90-100% : classe  9

per prova ho usato le concentrazioni di SiO2 (uno dei parametri piu' difficili da determinare da uno spettro IR nel rilevamento iperspettrale)

Prima di procedere i valori No data contenuti negli spettri e indicati dal valore -1.2300000e+034 sono stati convertiti in 0

provando con un rete neurale con delle impostazioni  non ottimizzate

==============================================================
import numpy as np
from numpy import genfromtxt

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers

x_train = genfromtxt("x_train.txt",delimiter=';',usecols=np.arange(0,479))
x_test = genfromtxt("x_test.txt",delimiter=';',usecols=np.arange(0,479))

y_train = genfromtxt("y_train.txt",delimiter=';',usecols=np.arange(0,1))
y_test = genfromtxt("y_test.txt",delimiter=';',usecols=np.arange(0,1))

y_train = y_train.transpose()

#print (x_train)


y_train = y_train.astype(int)
y_test = y_test.astype(int)
#print (y_train)

print ("====================")
print ("X trains" + str(x_train.shape))
print ("Y trains" + str(y_train.shape))
print ("X test" + str(x_test.shape))
print ("Y test" + str(y_test.shape))
print ("====================")

model = tf.keras.Sequential()
model.add(layers.Dense(128,activation='relu'))
model.add(layers.Dropout(0.2))
model.add(layers.Dense(10,activation='softmax'))

model.compile(optimizer='adam',loss='sparse_categorical_crossentropy',metrics=['accuracy'])

model.fit(x_train,y_train,epochs=16)
model.evaluate(x_test,y_test,verbose=2)
==============================================================

il risultato e' il seguente

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

X trains(60, 479)

Y trains(60,)

X test(11, 479)

Y test(11,)

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

2019-11-15 13:35:58.075622: I tensorflow/core/platform/cpu_feature_guard.cc:142] Your CPU supports instructions that this TensorFlow binary was not compiled to use: AVX2 FMA

2019-11-15 13:35:58.099810: I tensorflow/core/platform/profile_utils/cpu_utils.cc:94] CPU Frequency: 2294770000 Hz

2019-11-15 13:35:58.100500: I tensorflow/compiler/xla/service/service.cc:168] XLA service 0x38bfc90 executing computations on platform Host. Devices:

2019-11-15 13:35:58.100551: I tensorflow/compiler/xla/service/service.cc:175]   StreamExecutor device (0): Host, Default Version

Train on 60 samples

Epoch 1/16

60/60 [==============================] - 0s 8ms/sample - loss: 2.2551 - accuracy: 0.1000

Epoch 2/16

60/60 [==============================] - 0s 85us/sample - loss: 1.7444 - accuracy: 0.3833

Epoch 3/16

60/60 [==============================] - 0s 90us/sample - loss: 1.6373 - accuracy: 0.3333

Epoch 4/16

60/60 [==============================] - 0s 83us/sample - loss: 1.6465 - accuracy: 0.3167

Epoch 5/16

60/60 [==============================] - 0s 96us/sample - loss: 1.5436 - accuracy: 0.3833

Epoch 6/16

60/60 [==============================] - 0s 127us/sample - loss: 1.6167 - accuracy: 0.3667

Epoch 7/16

60/60 [==============================] - 0s 92us/sample - loss: 1.4574 - accuracy: 0.3833

Epoch 8/16

60/60 [==============================] - 0s 100us/sample - loss: 1.5219 - accuracy: 0.3500

Epoch 9/16

60/60 [==============================] - 0s 81us/sample - loss: 1.4713 - accuracy: 0.4000

Epoch 10/16

60/60 [==============================] - 0s 81us/sample - loss: 1.4284 - accuracy: 0.3500

Epoch 11/16

60/60 [==============================] - 0s 105us/sample - loss: 1.4529 - accuracy: 0.4167

Epoch 12/16

60/60 [==============================] - 0s 107us/sample - loss: 1.3913 - accuracy: 0.4167

Epoch 13/16

60/60 [==============================] - 0s 92us/sample - loss: 1.4574 - accuracy: 0.4000

Epoch 14/16

60/60 [==============================] - 0s 82us/sample - loss: 1.4370 - accuracy: 0.3500

Epoch 15/16

60/60 [==============================] - 0s 80us/sample - loss: 1.4666 - accuracy: 0.3667

Epoch 16/16

60/60 [==============================] - 0s 88us/sample - loss: 1.3726 - accuracy: 0.4500

11/1 - 0s - loss: 3.1148 - accuracy: 0.0000e+00

Impiegando i dati di concentrazione di Mg0 le cose migliorano

60/60 [==============================] - 1s 9ms/sample - loss: 1.8941 - accuracy: 0.3833
Epoch 2/16
60/60 [==============================] - 0s 194us/sample - loss: 1.4062 - accuracy: 0.6500
Epoch 3/16
60/60 [==============================] - 0s 167us/sample - loss: 1.1572 - accuracy: 0.6500
Epoch 4/16
60/60 [==============================] - 0s 187us/sample - loss: 1.1490 - accuracy: 0.6333
Epoch 5/16
60/60 [==============================] - 0s 173us/sample - loss: 1.0612 - accuracy: 0.6333
Epoch 6/16
60/60 [==============================] - 0s 195us/sample - loss: 1.0615 - accuracy: 0.6500
Epoch 7/16
60/60 [==============================] - 0s 193us/sample - loss: 1.0023 - accuracy: 0.6333
Epoch 8/16
60/60 [==============================] - 0s 164us/sample - loss: 1.0067 - accuracy: 0.6667
Epoch 9/16
60/60 [==============================] - 0s 221us/sample - loss: 0.9609 - accuracy: 0.6667
Epoch 10/16
60/60 [==============================] - 0s 185us/sample - loss: 0.9922 - accuracy: 0.6333
Epoch 11/16
60/60 [==============================] - 0s 163us/sample - loss: 0.9188 - accuracy: 0.6333
Epoch 12/16
60/60 [==============================] - 0s 234us/sample - loss: 0.9704 - accuracy: 0.6333
Epoch 13/16
60/60 [==============================] - 0s 176us/sample - loss: 0.9234 - accuracy: 0.6833
Epoch 14/16
60/60 [==============================] - 0s 202us/sample - loss: 0.9051 - accuracy: 0.6833
Epoch 15/16
60/60 [==============================] - 0s 181us/sample - loss: 0.8635 - accuracy: 0.7167
Epoch 16/16
60/60 [==============================] - 0s 171us/sample - loss: 0.8449 - accuracy: 0.7000
11/1 - 0s - loss: 2.3372 - accuracy: 0.2727