'''
Created on 28 mar 2017

@author: parisi
'''
import sys
import lsqnonneg as lsq
import numpy as np
#from _ast import If
import math
import os

__version__ = 0.1
__date__ = '2017-06-13'
__updated__ = '2017-06-13'

def getCol (filename, key):
    with open(filename) as f:
        f=[x.strip() for x in f if x.strip()]
        stationsNameList = f[0].split()
        data=[tuple(map(float,x.split())) for x in f[1:]]

    print ('Nomi Stazioni:')
    print (stationsNameList)
    print ('-------------')
    print ('numero stazioni: %s' % len(stationsNameList)) 
    print ('numero ore: %s' % len(data)) 
    print ('-------------')

    staz = []
    for nStaz in range(0,len(stationsNameList)):
        staz.append([x[nStaz] for x in data])
            
    '''
    for val_staz in staz:
        print val_staz
    print '-------------'
    '''
    '''
    col_1=[x[0] for x in data]
    col_2=[x[1] for x in data]
    col_3=[x[2] for x in data]
    col_4=[x[3] for x in data]
    col_5=[x[4] for x in data]
    col_6=[x[5] for x in data]
    col_7=[x[6] for x in data]
    col_8=[x[7] for x in data]
    col_9=[x[8] for x in data]
    col_10=[x[9] for x in data]
    print('col_1',col_1)
    print('col_2',col_2)
    print('col_3',col_3)
    print('col_4',col_4) 
    print('col_5',col_5)
    print('col_6',col_6)
    print('col_7',col_7)
    print('col_8',col_8)
    print('col_9',col_9)
    print('col_10',col_10)
    '''
    print (stationsNameList[key])
    '''
    for stationsName in stationsNameList:
        print stationsName
    '''

    for q, a in zip(staz[3], staz[key]):
        
        print ('ora')
        print ('Stazione %s?  Stazione4 %s.' % (q, a))        

def getStation (filename, key):
    with open(filename) as f:
        f=[x.strip() for x in f if x.strip()]
        stationsNameList = f[0].split()

    print ('Nomi Stazioni:')
    print (stationsNameList)
    print ('-------------')
    print ('numero stazioni: %s' % len(stationsNameList)) 
    print ('-------------')

    if key<=len(stationsNameList):
        print (stationsNameList[key])
        return stationsNameList[key]
    else:
        print ('errore key')

def getStationData (filename, key):
    with open(filename) as f:
        f=[x.strip() for x in f if x.strip()]
        stationsNameList = f[0].split()
        data=[tuple(map(float,x.split())) for x in f[1:]]

    print ('Nomi Stazioni:')
    print (stationsNameList)
    print ('-------------')
    print ('numero stazioni: %s' % len(stationsNameList)) 
    print ('numero ore: %s' % len(data) )
    print ('-------------')

    staz = []
    for nStaz in range(0,len(stationsNameList)):
        staz.append([x[nStaz] for x in data])
            
    '''
    print stationsNameList[key]
    print staz[key]
    '''
    
    for q, a in zip(staz[3], staz[7]):
        print ('Stazione %s?  Stazione4 %s.' % (q, a)) 

    if key==31:
        print ('OK')
        return zip(staz[0], staz[3], staz[7], staz[8])


#def main(filename , key ):
def main( ):
    program_name = os.path.basename(sys.argv[0])
    program_version = "v%s" % __version__
    program_build_date = "%s" % __updated__

    program_version_string = '%%prog %s (%s)' % (program_version, program_build_date)

    isWindows = 0
    if os.name == "posix":
        # Unix/Linux/MacOS/BSD/etc
        print ('Sistema Operativo Linux')
    elif os.name in ("nt", "dos", "ce"):
        # DOS/Windows
        print ('Sistema Operativo Windows')
        isWindows = 1
    else:
        # Fallback for other operating systems.
        print ('Sistema Operativo Altro')

    #with open('config_in.txt') as f:
    #print (os.path.join(os.getcwd(), 'workdir','config_in.txt'))
    with open(os.path.join(os.getcwd(), 'workdir','config_in.txt')) as f:
        f=[x.strip() for x in f if x.strip()]
        workingDir = f[0].split()[0]
        inputPath = f[1].split()[0]
        inputFile = f[2].split()[0]
        outputPath = f[3].split()[0]
        outputFile = f[4].split()[0]
        ctrlFile = f[5].split()[0]
        stationId = f[6].split()[0]
        elaborationId = f[7].split()[0]
        
    #print ('os.getcwd():')
    #print (os.getcwd())
    #f = open('D:/Sviluppo/Terraria/projects2017/dev_tipa/workdir/workdir/Output/getcwd.txt', 'w')
    #f.write("%s " % os.getcwd())
    #f.close()  
    
    #print (os.path.join(os.getcwd(), inputPath , inputFile))
    #with open(inputFile) as f:
    with open(os.path.join(os.getcwd(), workingDir, inputPath , inputFile)) as f:
        f=[x.strip() for x in f if x.strip()]
        stationsNameList = f[0].split()
        data=[tuple(map(float,x.split())) for x in f[1:]]

    #with open(filename) as f:
    #    f=[x.strip() for x in f if x.strip()]
    #    stationsNameList = f[0].split()
    #    data=[tuple(map(float,x.split())) for x in f[1:]]

    print ('Nomi Stazioni:')
    print (stationsNameList)
    print ('-------------')
    print ('numero stazioni: %s' % len(stationsNameList)) 
    print ('numero ore: %s' % len(data)) 
    print ('-------------')

    staz = []
    for nStaz in range(0,len(stationsNameList)):
        staz.append([x[nStaz] for x in data])
    
    Nt = len(data)
    #print ('Nt: %s' % Nt)

    kconf=1.2816;       #Fasce di confidenza pari al 80%
    idrog_unitario=5;   #na=n di ordinate dell'idrogramma unitario
    ore_ritardo=6;      #Ni=n di ore di ritardo Sesto Calende-Ponte Coperto
    #k=6;               #k=n di ore di ritardo Ponte Valenza-Ponte Coperto
    #kt=6;              #kt=n di ore di ritardo Montecastello-Ponte Coperto
    ore_previsione = ore_ritardo;   #Np=n di ore di previsione
    finestra_temporale = 30;    #Nfin=30; %n di ore ampiezza finestra temporale
    
    '''
    Operazione di media mobile sui dati per regolarizzarne 
    l'andamento(operazione preliminare). Si deve fissare un 
    valore di soglia oltre il quale non si procede alla media
    in modo tale da non tagliare troppo i picchi
    '''
    numero_media = 7;   #Numero di dati impiegati per la media
    soglia_pc=56.5;  #Soglia per i dati del Ponte Coperto di Pavia
    soglia_sc=3;    #Soglia per i dati del Sesto Calende
    soglia_va=3;    #Soglia per i dati di Valenza
    soglia_mc=5;    #Soglia per i dati di Montecastello
    
    #indice_valle = 0
    #soglia_valle = 0
    #soglia_monte_1 = 0
    #soglia_monte_2 = 0
    #soglia_monte_3 = 0
    #indice_monte_1 = 0
    #indice_monte_2 = 0
    #indice_monte_3 = 0
    # Assegnazione indice stazione valle in base alla elaborazione
    print (' stationId: %s' % stationId)
    if stationId=='1000006':
        indice_valle = 0
    elif stationId=='3000066':
        indice_valle = 1
    elif stationId=='3000079':
        indice_valle = 3
    else: print ('conf ERR')

    # Assegnazione soglie in base alla elaborazione e indici stazioni di monte
    if elaborationId=='31':
        print (' conf 31')
        #return zip(staz[0], staz[3], staz[7], staz[8])
        soglia_valle = soglia_pc
        soglia_monte_1 = soglia_sc
        soglia_monte_2 = soglia_va
        soglia_monte_3 = soglia_mc
        indice_monte_1 = 4
        indice_monte_2 = 7
        indice_monte_3 = 8
    elif elaborationId=='21':
        print (' conf 21')
        soglia_valle = soglia_pc
        soglia_monte_1 = soglia_sc
        soglia_monte_2 = soglia_mc
        indice_monte_1 = 4
        indice_monte_2 = 5
    elif elaborationId=='22':
        print (' conf 22')
        soglia_valle = soglia_pc
        soglia_monte_1 = soglia_sc
        soglia_monte_2 = soglia_mc # da sistemre la soglia
        indice_monte_1 = 4
        indice_monte_2 = 6
    else: print (' conf ERR')

    valle_medmob = []
    monte1_medmob = []
    monte2_medmob = []
    monte3_medmob = []
    
    indice_media=int(numero_media/2);
    
    #valle_medmob = staz[0];
    print (' indice_valle: %s ' % indice_valle )
    print (' valori_valle: %s '  % list(staz[indice_valle]))
    print (' indice_monte_1: %s ' % indice_monte_1 )
    print (' valori_monte_1: %s '  % list(staz[indice_monte_1]))
    print (' indice_monte_2: %s ' % indice_monte_2 )
    print (' valori_monte_2: %s '  % list(staz[indice_monte_2]))
    if elaborationId=='31':
        print (' indice_monte_3: %s ' % indice_monte_3 )
        print (' valori_monte_3: %s '  % list(staz[indice_monte_3]))
    # per duplicare una lista --> list(a) or a[:]
    valle_medmob = list(staz[indice_valle]);
    for i in range(indice_media+1,Nt-indice_media):
        if valle_medmob[i] <= soglia_valle:
            valle_medmob[i] = 0;
        
        for j in range(0,numero_media):
            valle_medmob[i] = valle_medmob[i]+staz[indice_valle][i+j-indice_media]

        valle_medmob[i] = valle_medmob[i]/numero_media
    
    #monte1_medmob = staz[3];
    monte1_medmob = list(staz[indice_monte_1]);
    for i in range(indice_media+1,Nt-indice_media):
        if monte1_medmob[i] <= soglia_monte_1:
            monte1_medmob[i] = 0;
        
        #print '------------ start'
        for j in range(0,numero_media):
            monte1_medmob[i] = monte1_medmob[i]+staz[indice_monte_1][i+j-indice_media]
        #    print 'i: %s' % monte1_medmob[i]
        
        #print '------------ fine'
        #print 'i: %s' % i
        #print 'monte1_medmob[i]: %s' % monte1_medmob[i]
        #print 'monte1_medmob[i]: %s' % numero_media
        monte1_medmob[i] = monte1_medmob[i]/numero_media
        #print 'monte1_medmob[i]/numero_media: %s' % monte1_medmob[i]
        #print 'staz[3][i]: %s' % staz[3][i]

    #monte2_medmob = staz[7];
    monte2_medmob = list(staz[indice_monte_2]);
    for i in range(indice_media+1,Nt-indice_media):
        if monte2_medmob[i] <= soglia_monte_2:
            monte2_medmob[i] = 0;
            
        for j in range(0,numero_media):
            monte2_medmob[i] = monte2_medmob[i]+staz[indice_monte_2][i+j-indice_media]

        monte2_medmob[i] = monte2_medmob[i]/numero_media
    
    if elaborationId=='31':
        #monte3_medmob = staz[8];
        monte3_medmob = list(staz[indice_monte_3]);
        for i in range(indice_media+1,Nt-indice_media):
            if monte3_medmob[i] <= soglia_monte_3:
                monte3_medmob[i] = 0;
                
            for j in range(0,numero_media):
                monte3_medmob[i] = monte3_medmob[i]+staz[indice_monte_3][i+j-indice_media]
    
            monte3_medmob[i] = monte3_medmob[i]/numero_media
       
    #print (' valle_medmob: %s' % valle_medmob)
    #print (' monte1_medmob: %s' % monte1_medmob)
    #print (' monte2_medmob: %s' % monte2_medmob)
    #print (' monte3_medmob: %s' % monte3_medmob)
    
    '''
    Calcolo differenze di altezze a Pavia, a Sesto a Valenza e a Montecastello 
    '''
    valle_diffmedmob = []
    monte1_diffmedmob = []
    monte2_diffmedmob = []
    if elaborationId=='31':
        monte3_diffmedmob = []
    
    for i in range(0,Nt-1):
        valle_diffmedmob.append( valle_medmob[i+1] - valle_medmob[i]);
        monte1_diffmedmob.append(monte1_medmob[i+1] - monte1_medmob[i]);
        monte2_diffmedmob.append(monte2_medmob[i+1] - monte2_medmob[i]);
        if elaborationId=='31':
            monte3_diffmedmob.append(monte3_medmob[i+1] - monte3_medmob[i]);
    #print 'monte1            : %s' % staz[3]
    #print 'monte1_medmob     : %s' % monte1_medmob
    #print 'monte1_diffmedmob : %s' % monte1_diffmedmob
    '''
    Prima parte: ricostruzione dell'andamento delle differenze di altezza 
                 a Pavia nella prima finestra temporale e previsione per le prime Np ore
    '''
    #monte1_fi = []
    monte1_fi = [[0 for col in range(0,idrog_unitario)] for row in range(0,finestra_temporale-(idrog_unitario+ore_previsione)+1)]
    #monte2_fi =  [[]]
    monte2_fi = [[0 for col in range(0,idrog_unitario)] for row in range(0,finestra_temporale-(idrog_unitario+ore_previsione)+1)]
    #monte3_fi = [[]]
    if elaborationId=='31':
        monte3_fi = [[0 for col in range(0,idrog_unitario)] for row in range(0,finestra_temporale-(idrog_unitario+ore_previsione)+1)]
    else:
        monte3_fi = 0
    
    ftot = [[0 for col in range(0,idrog_unitario*3)] for row in range(0,finestra_temporale-(idrog_unitario+ore_previsione)+1)]

    for t in range(idrog_unitario+ore_previsione,finestra_temporale+1): #11-30
        for i in range(0,idrog_unitario): #1-5
            #print ('---')
            #print (t-idrog_unitario-ore_previsione)
            #print (t - ore_previsione - i + 1)
            
            monte1_fi[t-idrog_unitario-ore_previsione][i] = monte1_diffmedmob[t - ore_previsione - i + 1]
            monte2_fi[t-idrog_unitario-ore_previsione][i] = monte2_diffmedmob[t - ore_previsione - i + 1]
            if elaborationId=='31':
                monte3_fi[t-idrog_unitario-ore_previsione][i] = monte3_diffmedmob[t - ore_previsione - i + 1]
            #monte2_fi[t,i] = monte2_diffmedmob[finestra_temporale - ore_ritardo + t - i + 1]
            #monte3_fi[t,i] = monte3_diffmedmob[finestra_temporale - ore_ritardo + t - i + 1]
    #ftot = list(set(monte1_fi + monte2_fi)) #monte1_fi + monte2_fi #+ monte3_fi
    #ftot = monte1_fi + monte2_fi + monte3_fi;
    
    #print (' monte1_fi: %s' % len(monte1_fi) )
    #print (' monte1_fi: %s' % len(monte1_fi[0])) 
    ftot = getUnion(monte1_fi, monte2_fi, monte3_fi)
    print (' ftot')
    #print (ftot)


    #y=np.array([np.array(ftoti) for ftoti in ftot])
    #np.array(ftot)
    #y = np.asarray(ftot)
    
    valle_diffmedmob_substring = np.asarray([valle_diffmedmob[i] for i in range(idrog_unitario+ore_previsione,finestra_temporale+1)])
    #valle_diffmedmob_substring = valle_diffmedmob_substring.T
    teta = lsq.lsqnonneg(ftot, valle_diffmedmob_substring);
    print (' teta ')
    #print ('teta: %s' % teta )
    
    qd1 = np.dot(ftot, teta)

    monte1_fi1 = [[0 for col in range(0,idrog_unitario)] for row in range(0,ore_previsione)]
    monte2_fi1 = [[0 for col in range(0,idrog_unitario)] for row in range(0,ore_previsione)]
    if elaborationId=='31':
        monte3_fi1 = [[0 for col in range(0,idrog_unitario)] for row in range(0,ore_previsione)]
    else:
        monte3_fi1 = 0
    ftot1 = [[0 for col in range(0,idrog_unitario*3)] for row in range(0,ore_previsione)]
    
    for t in range(0,ore_previsione): #1-6
        for i in range(0,idrog_unitario): #1-5
            monte1_fi1[t][i] = monte1_diffmedmob[finestra_temporale - ore_ritardo + t - i + 1]
            monte2_fi1[t][i] = monte2_diffmedmob[finestra_temporale - ore_ritardo + t - i + 1]
            if elaborationId=='31':
                monte3_fi1[t][i] = monte3_diffmedmob[finestra_temporale - ore_ritardo + t - i + 1]

    ftot1 = getUnion(monte1_fi1, monte2_fi1, monte3_fi1)
    
    #??? Cancello qr1???
    #qr1 = valle_diffmedmob_substring + qd1
    
    qd2 = np.dot(ftot1, teta)
    
    qr2 = [0 for row in range(0,ore_previsione+1)] # 1-7
    qr2[0] = valle_diffmedmob[finestra_temporale+1]
    for m in range(0, ore_previsione):
        qr2[m+1] = qr2[m] + qd2[m]
    
     
    '''
    #Seconda parte:ciclo con spostamento della finestra temporale di 1 ora in avanti, 
    #stima dei parametri e previsione per le successive Np ore
    '''
    
    valle_diffmedmob_substring_plus = np.asarray([valle_diffmedmob[i] for i in range(idrog_unitario+ore_previsione+1,finestra_temporale+1+1)])
    teta1 = lsq.lsqnonneg(ftot, valle_diffmedmob_substring_plus);
            
    monte1_fin1 = [[0 for col in range(0,idrog_unitario)] for row in range(0,ore_previsione)]
    monte2_fin1 = [[0 for col in range(0,idrog_unitario)] for row in range(0,ore_previsione)]
    if elaborationId=='31':
        monte3_fin1 = [[0 for col in range(0,idrog_unitario)] for row in range(0,ore_previsione)]
    else:
        monte3_fin1 = 0
    ftotn1 = [[0 for col in range(0,idrog_unitario*3)] for row in range(0,ore_previsione)]
    for t in range(0,ore_previsione): #1-6
        for i in range(0,idrog_unitario): #1-5
            monte1_fin1[t][i] = monte1_diffmedmob[finestra_temporale - ore_ritardo + t - i + 1]
            monte2_fin1[t][i] = monte2_diffmedmob[finestra_temporale - ore_ritardo + t - i + 1]
            if elaborationId=='31':
                monte3_fin1[t][i] = monte3_diffmedmob[finestra_temporale - ore_ritardo + t - i + 1]

    ftotn1 = getUnion(monte1_fin1, monte2_fin1, monte3_fin1)
    
    qdn3 = np.dot(ftotn1, teta1)
    
    qr3 = [0 for row in range(0,ore_previsione+1)] # 1-7
    
    #qr3[0] = valle_diffmedmob[finestra_temporale+1]
    qr3[0] = staz[indice_valle][finestra_temporale+16] # dato staz valle con 16 is j 11-16)
    for m in range(0, ore_previsione):
        qr3[m+1] = qr3[m] + qdn3[m]
        
    print (' qr3 ')
    #print ('qr3: %s' % qr3 )
    
    for j in range(Nt - finestra_temporale - ore_previsione, Nt - finestra_temporale - 1):
        # fase taratura
        monte1_fin2 = [[0 for col in range(0,idrog_unitario)] for row in range(0,finestra_temporale)]
        monte2_fin2 = [[0 for col in range(0,idrog_unitario)] for row in range(0,finestra_temporale)]
        if elaborationId=='31':
            monte3_fin2 = [[0 for col in range(0,idrog_unitario)] for row in range(0,finestra_temporale)]
        else:
            monte3_fin2 = 0
        ftotn2 = [[0 for col in range(0,idrog_unitario*3)] for row in range(0,ore_previsione)]
        for t in range(j,finestra_temporale+j-1): #11-16
            for i in range(0,idrog_unitario): #1-5
                monte1_fin2[t-j][i] = monte1_diffmedmob[t - ore_ritardo - i + 1]
                monte2_fin2[t-j][i] = monte2_diffmedmob[t - ore_ritardo - i + 1]
                if elaborationId=='31':
                    monte3_fin2[t-j][i] = monte3_diffmedmob[t - ore_ritardo - i + 1]

        ftotn2 = getUnion(monte1_fin2, monte2_fin2, monte3_fin2)
        
        valle_diffmedmob_substring2 = np.asarray([valle_diffmedmob[i] for i in range(j,finestra_temporale+j)])
        teta1 = lsq.lsqnonneg(ftotn2, valle_diffmedmob_substring2);
    
        qdn = np.dot(ftotn2, teta1)
        qrn = [staz[indice_valle][row] for row in range(j,finestra_temporale+j)] + qdn
        yeffn = [staz[indice_valle][row] for row in range(j+1,finestra_temporale+j+1)]
        
        v = qrn - yeffn
        dpf = len(v)
        
        varv = 0
        for jj in range(1,dpf):
            varv = varv + np.dot(v[jj],v[jj])
        
        vv =  varv/dpf
        
        # Covarianze #
        cov = [0 for row in range(0,dpf)]
        for im in range(0,dpf):
            covv = 0
            for il in range (0, dpf-im):
                covv = covv + np.dot(v[il],v[il + im])
            
            cov[im] = covv/(dpf-im)
        
        #print (' Covarianze '  )
        #print (' Covarianze: %s' % cov)
        
        # matrice varianza errore v #
        mcov = [[0 for col in range(0,dpf)] for row in range(0,dpf)]
        for mx in range(0,dpf):
            for inn in range(0,dpf):
                if inn == mx :
                    mcov[mx][inn] = vv
                elif inn>mx :
                    mcov[mx][inn] = cov[inn-mx]
                elif inn<mx :
                    mcov[mx][inn] = cov[mx-inn]
        
        # calcolo della matrice varianza dell'errore di stima dei parametri teta
        #vteta1 = np.linalg.inv(np.dot(ftotn2,np.linalg.inv(mcov),ftotn2))
        #print (' ftotn2: %s' % ftotn2 )
        #print (' ftotn2.T: %s' % ftotn2.T)
        #print (' matrice varianza errore di stima dei parametri teta' )
        #print (' matrice varianza errore di stima dei parametri teta: %s' % mcov)
        print (' Determinante mcov: %s' % np.linalg.det(mcov))
        if(np.linalg.det(mcov) == 0):
             print (' LA MATRICE HA DETERMINANTE ZERO: la matrice mcov non e invertibile')
             continue
        A0 = np.linalg.inv(mcov)
        #print (' A0 ---------------------------------------------')
        #print (' A0: %s' % A0)
        A1 = np.dot(A0,ftotn2)
        #print (' A1 ---------------------------------------------')
        #print (' A1: %s' % A1)
        A2 = np.dot(ftotn2.T,A1)
        #print (' A2 ---------------------------------------------')
        #print (' A2: %s' % A2)
        
        print (' Determinante A2: %s' % np.linalg.det(A2))
        if(np.linalg.det(A2) == 0):
             print (' LA MATRICE HA DETERMINANTE ZERO: la matrice A2 non e invertibile')
             continue
        A3 = np.linalg.inv(A2)
        #print (' A3 ---------------------------------------------')
        #print (' A3: %s' % A3)
        #vteta1 = np.linalg.inv(np.dot(ftotn2.T,np.dot(np.linalg.inv(mcov),ftotn2)))
        vteta1 = A3
        
        # fase previsione
        for t in range(0,ore_previsione): #1-6
            for i in range(0,idrog_unitario): #1-5
                monte1_fin1[t][i] = monte1_diffmedmob[finestra_temporale + j - ore_ritardo + t - i + 1]
                monte2_fin1[t][i] = monte2_diffmedmob[finestra_temporale + j - ore_ritardo + t - i + 1]
                if elaborationId=='31':
                    monte3_fin1[t][i] = monte3_diffmedmob[finestra_temporale + j - ore_ritardo + t - i + 1]

        ftotn4 = getUnion(monte1_fin1, monte2_fin1, monte3_fin1)
            
        qdn3 = np.dot(ftotn1, teta1)
        qr3 = [0 for row in range(0,ore_previsione+1)] # 1-7
    
        #qr3[0] = valle_diffmedmob[finestra_temporale+1]
        qr3[0] = staz[indice_valle][finestra_temporale+16] # dato staz valle con 16 is j 11-16)
        for m in range(0, ore_previsione):
            qr3[m+1] = qr3[m] + qdn3[m]
        
        #print ('qr3: %s' % qr3 )
        
        # calcolo della matrice varianza delle previsioni delle differenze di altezza idrometrica
        A4 = ftotn4.T
        A5 = np.dot(vteta1,A4)
        A6 = np.dot(ftotn4,A5)
        #vprev = np.dot(ftotn4,vteta1,ftotn4')
        vprev = A6
        
        zup = [0 for row in range(0,ore_previsione+1)]
        zdn = [0 for row in range(0,ore_previsione+1)]
        var = [0 for row in range(0,ore_previsione)]
        sc = [0 for row in range(0,ore_previsione)]
        
        zup[0] = qr3[0]
        zdn[0] = qr3[0]
        for inn in range(0, ore_previsione):
            var1=0
            for jn in range(0, inn):
                var1 = var1 + abs(vprev[jn][jn])
            var[inn] = var1
            sc[inn] = math.sqrt( var[inn] )
            zup [inn+1] = qr3 [inn + 1] + np.dot(kconf,sc[inn])
            zdn [inn+1] = qr3 [inn + 1] - np.dot(kconf,sc[inn])
        
        #print ('qr3: %s' % qr3 )
        #print ('zup: %s' % zup )
        #print ('zdn: %s' % zdn )
        #print ('FINE ciclo J ----------------------- %s' % j )
    
    #print ('FINE FOR J-----------------------' )
    print ('qr3: %s' % qr3 )
    try: 
        # se zup non e definito vuol dire che una matrice non e invertibile
        zup
    except NameError:
        # gestisco il caso di una matrice non e invertibile terminando l'algoritmo (scrivendo il ctrlFile)
        flag_errore = 1
        print (' ************************************************************')
        print (' * NON E STATO POSSIBILE COMPLETARE LA SEGUENTE SIMULAZIONE *')
        print (' * statione: %s                                        *' % stationId)
        print (' * elaborazione: %s                                         *' % elaborationId)
        print (' ************************************************************')
        print('output: %s ' % os.path.join(os.getcwd(), workingDir, outputPath , ctrlFile))
        f = open(os.path.join(os.getcwd(), workingDir, outputPath , ctrlFile),'w')
        f.write("%s " % flag_errore)
        f.close()
        return 0
    
    print ('zup: %s' % zup )
    print ('zdn: %s' % zdn )
     
    #print ('Arrivato qui') 
    # prova per eseguire da java il py
    #f = open('livelli.txt', 'w')
    #f = open('D:/Sviluppo/Terraria/projects2017/devTipaPy/src/livelli.txt', 'w')
    #f = open('D:/Sviluppo/Terraria/projects2017/dev_tipa/workdir/workdir/Output/livelli_prev.txt', 'w')
    print('output: %s ' % os.path.join(os.getcwd(), workingDir, outputPath , outputFile))
    f = open(os.path.join(os.getcwd(), workingDir, outputPath , outputFile),'w')
    for p in range(0, ore_previsione+1):
        f.write("%s " % zdn[p])
        f.write("%s\n" % zup[p])
    f.close()
    
    # controllo
    flag_errore = 1
    if( all(x for x in zdn) and all(x for x in zup) ):
        flag_errore = 0 
    else: 
        flag_errore = 1

    print('output: %s ' % os.path.join(os.getcwd(), workingDir, outputPath , ctrlFile))
    f = open(os.path.join(os.getcwd(), workingDir, outputPath , ctrlFile),'w')
    f.write("%s " % flag_errore)
    f.close()

    #return flag_errore


def getUnion(a, b, c=0):
    # find all unique values   -len(monte1_fi[0]
    maxRow = len(a);
    maxCol = len(a[0]);
    results=[[0 for col in range(0,maxCol)] for row in range(0,maxRow)]
    for colonna in range(0,len(a)): #1-20
        riga_results=[]
        for riga  in range(0,len(a[0])): #1-5
            #if result in unique_results:
            #    continue
            riga_results.append(a[colonna][riga])
        for riga  in range(0,len(a[0])): #1-5
            riga_results.append(b[colonna][riga])
        
        # gestisce il caso con tre matrici start
        if c!=0:
            for riga  in range(0,len(a[0])): #1-5
                riga_results.append(c[colonna][riga])
        # gestisce il caso con tre matrici end

        #results.append(riga_results)
        results[colonna]=np.asarray(riga_results)

    return np.asarray(results)

       
if __name__ == '__main__':
    
    #getCol('altezze.txt',3)
    #getStation('altezze.txt',3)
    #print getStationData('altezze.txt',31)
    
    #livGuardia = 58.5
    
    #sys.exit(main('altezze.txt',31))
    #stations = main('altezze.txt',31)
    #stations = main('D:/Sviluppo/Terraria/projects2017/devTipaPy/src/altezze.txt',31)
    #stations = main('D:/Sviluppo/Terraria/projects2017/dev_tipa/workdir/workdir/Input/altezze.txt',31)
    sys.exit(main())
    #print stations