#exec(open("Bolsas.py").read()) # Para anotar junto de um ponto, ver: # https://www.geeksforgeeks.org/python/line-chart-in-matplotlib-python/ import gudhi import gudhi.representations import numpy as np import matplotlib.pyplot as plt from scipy.spatial.distance import pdist, squareform import random import os import sys import warnings warnings.filterwarnings('ignore') if sys.platform == "win32": exec(open("smeUtils.py").read()) exec(open("smeGudhi.py").read()) else: from smeUtils import * from smeGudhi import * #exec(open("smeUtils.py").read()) #exec(open("smeGudhi.py").read()) exec(open("BolsasAux.py").read()) exec(open("auxMatlabEx.py").read()) # Exemplo: # # rc = BolsasFile("NY_20080102_20081231__319_dim250.txt", edgeMax=0.2, maxDim=4, sparse=None, skip=0, cols=[0,1,2,3], persLim=0.03) # rc = BolsasFile("/db/Economia/DataOrg2026/NY_20080102_20081231__319_dim250.txt", edgeMax=0.12, maxDim=3, sparse=None, skip=0, cols=[1,2,3], persLim=0.03) # rc = Bolsas(1998, 0, edgeMax=0.2, maxDim=4, sparse=None, skip=0, cols=[1,2,3], persLim=0.03) # rc = Bolsas(1998, 0, edgeMax=0.15, maxDim=4, sparse=None, skip=0, cols=[0,1,2,3], persLim=0.03) # fnames = { 1989:["NY_19890421_19891229__318_dim175.txt","NY_19890421_19890630__318_dim49.txt" ,"NY_19890703_19891229__318_dim125.txt"], 1990:["NY_19900102_19901231__318_dim252.txt","NY_19900102_19900629__318_dim125.txt","NY_19900702_19901231__318_dim126.txt"], 1991:["NY_19910102_19911231__318_dim252.txt","NY_19910102_19910628__318_dim124.txt","NY_19910701_19911231__318_dim127.txt"], 1992:["NY_19920102_19921231__318_dim253.txt","NY_19920102_19920630__318_dim125.txt","NY_19920701_19921231__318_dim127.txt"], 1993:["NY_19930104_19931231__318_dim252.txt","NY_19930104_19930630__318_dim124.txt","NY_19930701_19931231__318_dim127.txt"], 1994:["NY_19940103_19941230__318_dim251.txt","NY_19940103_19940630__318_dim124.txt","NY_19940701_19941230__318_dim126.txt"], 1995:["NY_19950103_19951229__318_dim251.txt","NY_19950103_19950630__318_dim125.txt","NY_19950703_19951229__318_dim125.txt"], 1996:["NY_19960102_19961231__318_dim253.txt","NY_19960102_19960628__318_dim125.txt","NY_19960701_19961231__318_dim127.txt"], 1997:["NY_19970102_19971231__318_dim252.txt","NY_19970102_19970630__318_dim124.txt","NY_19970701_19971231__318_dim127.txt"], 1998:["NY_19980102_19981231__318_dim251.txt","NY_19980102_19980630__318_dim123.txt","NY_19980701_19981231__318_dim127.txt"], 1999:["NY_19990104_19991231__318_dim251.txt","NY_19990104_19990630__318_dim123.txt","NY_19990701_19991231__318_dim127.txt"], 2000:["NY_20000103_20001229__318_dim251.txt","NY_20000103_20000630__318_dim125.txt","NY_20000703_20001229__318_dim125.txt"], 2001:["NY_20010102_20011231__318_dim247.txt","NY_20010102_20010629__318_dim124.txt","NY_20010702_20011231__318_dim122.txt"], 2002:["NY_20020102_20021231__318_dim251.txt","NY_20020102_20020628__318_dim123.txt","NY_20020701_20021231__318_dim127.txt"], 2003:["NY_20030102_20031231__318_dim251.txt","NY_20030102_20030630__318_dim123.txt","NY_20030701_20031231__318_dim127.txt"], 2004:["NY_20040102_20041231__318_dim251.txt","NY_20040102_20040630__318_dim123.txt","NY_20040701_20041231__318_dim127.txt"], 2005:["NY_20050103_20051230__318_dim251.txt","NY_20050103_20050630__318_dim124.txt","NY_20050701_20051230__318_dim126.txt"], 2006:["NY_20060103_20061229__318_dim250.txt","NY_20060103_20060630__318_dim124.txt","NY_20060703_20061229__318_dim125.txt"], 2007:["NY_20070103_20071231__318_dim250.txt","NY_20070103_20070629__318_dim123.txt","NY_20070702_20071231__318_dim126.txt"], 2008:["NY_20080102_20081231__318_dim250.txt","NY_20080102_20080630__318_dim124.txt","NY_20080701_20081231__318_dim125.txt"], 2009:["NY_20090102_20090420__318_dim73.txt"] } Sectores = {'10':'Energy','15':'Materials','20':'Industrials','25':'Cons. Discret.','30':'Cons. Staples', '35':'Health Care', '40':'Financials','45':'Inform Technol.','50':'Telecom. Services','55':'Utilities', '60':'Real Estate', '95':'Não Classif.'} def graph_NascMorte (rc, x, y): if not x in ['D','P','C','M']: print("Valor inválido para 'x'") return if len(y) == 0: print("Valor inválido para 'y'") return vMais = [] vMenos = [] vPers = [] vPCent = [] vDim = [] for v1 in rc.holes: vDim.append(v1['dimensão']) vMenos.append(v1['morte'] - v1['nascimento']) vMais.append((v1['morte'] + v1['nascimento']) / 2.) vPers.append(v1['persistência']) vPCent.append(v1['ponto_central']) print(f"x; {x}") vx1 = [] if x == 'D': vx1 = vDim elif x == 'P': vx1 = vPers elif x == 'C': vx1 = vMais elif x == 'M': vx1 = vMenos else: print("***** Valor inválido para 'x'") return if 'D' in y: plt.plot(vx1, vDim) if 'P' in y: plt.plot(vx1, vPers) if 'C' in y: plt.plot(vx1, vMais) if 'M' in y: plt.plot(vx1, vMenos) plt.show() def graph_NascMorte2 (rc): vMais = [] vMenos = [] vPers = [] vPCent = [] for x in rc.holes: vMais.append(x['nascimento'] + x['morte']) vMenos.append(x['nascimento'] - x['morte']) vPers.append(x['persistência']) vPCent.append(x['ponto_central']) plt.plot(vMenos, vMais) plt.plot(vMenos, vPers) plt.plot(vMenos, vPCent) plt.grid(visible=True) plt.show() return {'vMais':vMais, 'vMenos':vMenos} # Exemplo: # rc = BolsasFile("/db/Economia/DataOrg2026/NY_20080102_20081231__319_dim250.txt", edgeMax=0.12, maxDim=3, sparse=None, skip=0, cols=[1,2,3], persLim=0.03) def BolsasFile (fname, edgeMax=0.2, maxDim=2, sparse=None, skip=0, cols=[0,1], persLim=0.03, graphShow=False, graphTitle='', verbose=1, printFName=False, edgeMaxVol=True, persLimEdgeMax=0.1): # !!!!! #pts = smeFileReadNumRand ("/db/Gudhi/TDA-tutorial-master/datasets/tore3D_1307.off", skip=2, nRand=200) #pts = smeFileReadNumRand("./Coord30.txt", cols=[1,2], nRand=100) #pts = gudhi.read_points_from_off_file(off_file="vDim02_01_to_02_test.txt") #pts = [[0.1, 0.2, 0.3, 0.4, 0.5], [0.13, 0.12, 0.23, 0.14, 0.45], [0.23, 0.32, 0.13, 0.04, 0.25], [0.13, 0.32, 0.45, 0.21, 0.34], [0.03, 0.34, 0.25, 0.41, 0.19]] #pts = gera2circulos(centros=[[-3,0],[3,0]]) #pts = gera2esferas(qtPts=500, centros=[[-3,0],[3,0]], ruido=0.0) #pts = readDataFileNY(cols=cols, skip=skip) #fname = "NY_19890420__20090420__319_dim250.txt" #pts = smeFileReadNum(fname, cols=cols, skip=skip) vData = bolsaFileRead(fname, cols=cols, skip=skip) if edgeMaxVol == True: edgeMax = pow(vData[1][0]['vol'], 1/3) # !!!!! if persLimEdgeMax > 0: persLim = persLimEdgeMax * edgeMax # !!!!! pts = vData[0] cols = vData[2] if len(pts) == 0 or len(cols) == 0: print("\n***** Não existem colunas v álidas.\n") return -1 #print(f"fname: {fname}") #print(f"edgeMax: {edgeMax} ; maxDim{maxDim} ; persLim: {persLim}") #print(f"Len(pts): {len(pts)} ; len(cols): {len(cols)}") #pts = gera2esferas(qtPts=2000, centros=[[-2,0,0],[2,0,0]], raio=1, ruido=0.01) pts = np.vstack(pts) if graphShow: smeGhudiPlotRipsComplexFromPts(pts, max_edge_length=edgeMax, title=f"Bolsas {graphTitle}") for r in range(0, 20, 15): x = float(r) / 10 xx = f"Raio: {x:.2f}" #smeGhudiPlotRipsComplexFromPts(pts, max_edge_length=x, title=xx) #plt.show() # Cria objecto para guardar RipsComplex rc = criaRipsComplex(pts, fname, edgeMax, sparse=sparse, maxDim=maxDim, cols=cols, skip=skip, persLim=persLim) if rc.error < 0: print(f"Erro a executar 'criaRipsComplex'. Erro nº: {rc.error} ; File:\n fname: {fname}") return rc rc.infos = vData[1] rc.infos[0]['fname']=fname rc.vol = rc.infos[0]['vol'] if graphShow: rc.visualizar_pontos("Pontos") #print(f"persistence: {rc.persistence}") #5/0 if printFName: print(f"fname: {fname}") rc.barcode = rc.simplex_tree.persistence_intervals_in_dimension(0) for d in range(1, rc.maxDim + 1): dim_barcode = rc.simplex_tree.persistence_intervals_in_dimension(d) rc.barcode = np.vstack([rc.barcode, dim_barcode]) if len(rc.barcode) > 0 else dim_barcode # barcodes diagram if graphShow: rc.visualizar_diag_persistencia() # persistence diagrams if graphShow: gudhi.plot_persistence_diagram(rc.persistence) #plt.show() # Números de Betti: rc.numBetti = rc.simplex_tree.betti_numbers() if verbose > 0: print(f"Números de Betti: {rc.numBetti}") # Cálculo dos buracos: rc.holes = rc.findHoles(limiar_persistencia=persLim) rc.numHoles = len(rc.holes) if verbose > 0: print(f"Buracos significativos encontrados: {rc.numHoles}") plt.show() # Dados if verbose > 0: print(f"Volume: {rc.vol} ; distMédia: {rc.meanDist:.4f} +- {rc.meanDistErr:.4f}") return rc def Bolsas (f1, f2, edgeMax=0.2, maxDim=2, sparse=None, skip=0, cols=[0,1], persLim=0.03, graphShow=True): fname = "/db/Economia/DataOrg2026/" + fnames[f1][f2] print(f"fname: {fname}") rc = BolsasFile (fname, edgeMax, maxDim, sparse, skip, cols, persLim, graphShow) return rc # rc = BolsasGetMat("20-Jan-2005", "10-Feb-2005", "/tmp/Economia", edgeMax=0.2, maxDim=4, sparse=None, skip=0, cols=[0,1,2,3], persLim=0.03) def BolsasGetMat (date1, date2, dirOut, edgeMax=0.2, maxDim=2, sparse=None, skip=0, cols=[0,1], persLim=0.03, graphShow=True, delFiles=False): fname = runMatlabFromPython(date1, date2, dirOut, delFiles) print(f"fname: {fname}") rc = BolsasFile (fname, edgeMax, maxDim, sparse, skip, cols, persLim, graphShow) return rc def loopBolsasFromTo_getFiles (tpLoop, f1, f2): if not f1 in fnames: print("\n ***** Erro: O primeiro ano que existe é 1989\n") return anos = range(f1, f2+1) vFiles = [] if tpLoop == 'Ano': filesAno = [1, 2] else: filesAno = [0] n1 = 0 for ano in anos: if not ano in fnames: break lenFAno = len(fnames[ano]) if tpLoop == 'Ano': dic1 = {'N':n1, 'Ano':ano, 'qt':1, 'pos':0, 'file':fnames[ano][0], 'rc':0} vFiles.append(dic1) n1 += 1 else: if lenFAno > 1: vFiles.append({'N':n1, 'Ano':ano, 'qt':2, 'pos':1, 'file':fnames[ano][1], 'rc':0}) if lenFAno > 2: vFiles.append({'N':n1+1, 'Ano':ano, 'qt':2, 'pos':2, 'file':fnames[ano][2], 'rc':0}) n1 += 2 return vFiles # # tpLoop: 'Ano', 'Sem' # f1: Primeiro ano # f2: Último ano # # Retorna um vector de estruturas com todas as informações # # vrc = BolsasFromTo('Ano', 1993, 2004, edgeMax=0.2, maxDim=4, sparse=None, skip=0, cols=[0,1,2,3], persLim=0.03) def BolsasFromTo (tpLoop, f1, f2, edgeMax=0.2, maxDim=2, sparse=None, skip=0, cols=[0,1], persLim=0.03): n1 = 0 vData = loopBolsasFromTo_getFiles (tpLoop, f1, f2) dname = "/db/Economia/DataOrg2026/" n1 = 0 for x1 in vData: print(f"Ano: {x1['Ano']} ; Qt: {x1['qt']} ; File: {x1['file']}") fname = dname + x1['file'] print(f"fname: {fname}") vData[n1]['rc'] = BolsasFile (fname, edgeMax, maxDim, sparse, skip, cols, persLim, graphShow=False) if vData[n1]['rc'] == -1: print(f"\n***** Erro(BolsasFromTo) vData[{n1}] - BolsasFile retorna '-1'") return -1 n1 += 1 return vData # BolsasExecSeq("/db/Economia/DataOrg2026_Testes/DataFiles_Len06_01-Jul-1989_31-Jul-1989.txt", edgeMax=0.2, maxDim=2, sparse=None, skip=0, cols=[0,1], persLim=0.03, graphShow=True): def BolsasExecSeq (fname, edgeMax=0.2, maxDim=2, sparse=None, skip=0, cols=[0,1], persLim=0.03, graphShow=False): f1 = open(fname, "rt") if not f1: print("***** Erro: Não conseguiu abrir o ficheiro") return -1 linha = f1.readline() while linha: linha = linha.strip() spLine = linha.split("|") lenSpLine = len(spLine) if lenSpLine != 3: print(f"***** Erro a ler o ficheiro '{fname}'") return -1 date1 = spLine[0] date2 = spLine[1] fname2 = spLine[2] title1 = date1 + " - " + date2 print(f"File: {fname2}") rc = BolsasFile(fname2, edgeMax=0.2, maxDim=3, sparse=None, skip=0, cols=[1,2,3], persLim=0.03, graphShow=graphShow, graphTitle=title1) if rc == -1: print("\n***** Erro(BolsasExecSeq): BolsasFile retorna '-1'") linha = f1.readline() # rc = BolsasFile("/db/Economia/DataOrg2026/NY_20010101_20011231__148_dim147.txt", edgeMax=0.12, maxDim=3, sparse=None, skip=0, cols=[1,2,3], persLim=0.1) # L = BolsasLandscape(rc.persistence, dim=1) # L = BolsasLandscape (rc.persistence, dim=2) def BolsasLandscape2 (persistence, dim=-1, resolution=100): from sklearn.preprocessing import MinMaxScaler from gudhi.representations import (DiagramSelector,Clamping,Landscape,Silhouette,BettiCurve,ComplexPolynomial,TopologicalVector,DiagramScaler,BirthPersistenceTransform,PersistenceImage,PersistenceWeightedGaussianKernel,Entropy,PersistenceScaleSpaceKernel,SlicedWassersteinDistance,SlicedWassersteinKernel,PersistenceFisherKernel,WassersteinDistance,) D1 = [] if dim == -1: for p1 in persistence: D1.append(p1[1]) else: for p1 in persistence: if p1[0] == dim: D1.append(p1[1]) D1 = np.array(D1) proc1 = DiagramSelector(use=True, point_type="finite") proc2 = DiagramScaler(use=True, scalers=[([0, 1], MinMaxScaler())]) proc3 = DiagramScaler(use=True, scalers=[([1], Clamping(maximum=0.9))]) D1 = proc3(proc2(proc1(D1))) LS = Landscape(resolution=resolution) L = LS(D1) plt.plot(L) # plt.plot(L[:100]) # plt.plot(L[100:200]) # plt.plot(L[200:300]) plt.title("Landscape") plt.grid(visible=True) plt.show() return L def BolsasLandscape (persistence): from sklearn.preprocessing import MinMaxScaler from gudhi.representations import (DiagramSelector,Clamping,Landscape,Silhouette,BettiCurve,ComplexPolynomial,TopologicalVector,DiagramScaler,BirthPersistenceTransform,PersistenceImage,PersistenceWeightedGaussianKernel,Entropy,PersistenceScaleSpaceKernel,SlicedWassersteinDistance,SlicedWassersteinKernel,PersistenceFisherKernel,WassersteinDistance,) D1=[] # guarda nasc e morte de cada ponto do diag.persistencia em dim>0 for p in rc.persistence: dim=p[0] bd=p[1] de=bd[1] if dim > 0 and not np.isinf(de): D1.append(p[1]) X=np.array(D1) proc1 = DiagramSelector(use=True, point_type="finite") proc2 = DiagramScaler(use=True, scalers=[([0, 1], MinMaxScaler())]) proc3 = DiagramScaler(use=True, scalers=[([1], Clamping(maximum=0.9))]) D = proc3(proc2(proc1(X))) LS = Landscape(resolution=1000) L = LS(D) plt.plot(L[:1000]) plt.plot(L[1000:2000]) plt.plot(L[2000:3000]) V=rc.infos[0]['vol'] Y=[ V, [rc.fname[44:61]]] plt.title(Y) plt.grid(visible=True) plt.show()