From 65acdf03e3164c7b2594934ccd3ca43f8af69eab Mon Sep 17 00:00:00 2001 From: Shayan Bastani <32581929+shayanb2004@users.noreply.github.com> Date: Wed, 21 Sep 2022 23:39:22 +0430 Subject: [PATCH 1/2] updated for python 3.x --- financeExample.py | 4 +- ggs.py | 302 +++++++++++++++++++++++----------------------- 2 files changed, 154 insertions(+), 152 deletions(-) diff --git a/financeExample.py b/financeExample.py index 0a57d05..a0adacb 100644 --- a/financeExample.py +++ b/financeExample.py @@ -20,8 +20,8 @@ meancovs = GGSMeanCov(data, breakpoints = bp10, lamb = 1e-4, features = feats) -print "Breakpoints are at", bps -print "Objectives are", objectives +print ("Breakpoints are at", bps) +print ("Objectives are", objectives) # Plot objective vs. number of breakpoints plotVals = range(len(objectives)) diff --git a/ggs.py b/ggs.py index 5997ebb..8e1d765 100644 --- a/ggs.py +++ b/ggs.py @@ -14,101 +14,101 @@ from sys import platform as _platform - -#Find K breakpoints on the data at a specific lambda -#Returns: The K breakpoints, along with all intermediate breakpoints (for k < K) and their corresponding +# Find K breakpoints on the data at a specific lambda +# Returns: The K breakpoints, along with all intermediate breakpoints (for k < K) and their corresponding # covariance-regularized maximum likelihoods -def GGS(data, Kmax, lamb, features = [], verbose = False): +def GGS(data, Kmax, lamb, features=[], verbose=False): data = data.T - #Select the desired features + # Select the desired features if (features == []): features = range(data.shape[1]) - data = data[:,features] - m,n = data.shape + data = data[:, features] + m, n = data.shape - #Initialize breakpoints - breaks = [0,m+1] + # Initialize breakpoints + breaks = [0, m + 1] breakPoints = [breaks[:]] - plotPoints = [calculateLikelihood(data, breaks,lamb)] + plotPoints = [calculateLikelihood(data, breaks, lamb)] - #Start GGS Algorithm + # Start GGS Algorithm for z in range(Kmax): - numBreaks = z+1 + numBreaks = z + 1 newInd = -1 newVal = +1 - #For each segment, find breakpoint and increase in LL + # For each segment, find breakpoint and increase in LL for i in range(numBreaks): - tempData = data[breaks[i]:breaks[i+1], :] + tempData = data[breaks[i]:breaks[i + 1], :] ind, val = addBreak(tempData, lamb) - if(val < newVal): + if (val < newVal): newInd = ind + breaks[i] newVal = val - #Check if our algorithm is finished - if(newVal == 0): - print "We are done adding breakpoints!" - print breaks + # Check if our algorithm is finished + if (newVal == 0): + print("We are done adding breakpoints!") + print(breaks) return breaks, plotPoints - #Add new breakpoint + # Add new breakpoint breaks.append(newInd) breaks.sort() if (verbose == True): - print "Breakpoint occurs at sample number: ", newInd, ", LL = ", newVal - print len(breaks) - 2, breaks + print("Breakpoint occurs at sample number: ", newInd, ", LL = ", newVal) + print(len(breaks) - 2, breaks) - #Adjust current locations of the breakpoints - breaks = adjustBreaks(data,breaks,[newInd],lamb,verbose)[:] + # Adjust current locations of the breakpoints + breaks = adjustBreaks(data, breaks, [newInd], lamb, verbose)[:] - #Calculate likelihood - ll = calculateLikelihood(data,breaks,lamb) + # Calculate likelihood + ll = calculateLikelihood(data, breaks, lamb) breakPoints.append(breaks[:]) plotPoints.append(ll) return breakPoints, plotPoints -#Run cross-validation up to Kmax for a set of lambdas -#Return: train and test set likelihood for every K, lambda -def GGSCrossVal(data, Kmax=25, lambList = [0.1, 1, 10], features = [], verbose = False): + +# Run cross-validation up to Kmax for a set of lambdas +# Return: train and test set likelihood for every K, lambda +def GGSCrossVal(data, Kmax=25, lambList=[0.1, 1, 10], features=[], verbose=False): data = data.T if (features == []): features = range(data.shape[1]) - data = data[:,features] + data = data[:, features] origSize, n = data.shape np.random.seed(0) - ordering = range(origSize) + ordering = list(range(origSize)) random.shuffle(ordering) trainTestResults = [] - #For each lambda, run the 10 folds in parallel - numProcesses = min(multiprocessing.cpu_count(),10 ) - pool = multiprocessing.Pool(processes = numProcesses) + # For each lambda, run the 10 folds in parallel + numProcesses = min(multiprocessing.cpu_count(), 10) + pool = multiprocessing.Pool(processes=numProcesses) for lamb in lambList: mseList = [] trainList = [] - returnList = pool.map(multi_run_wrapper, [(0,data, Kmax, lamb, verbose, origSize, n, ordering), - (1,data, Kmax, lamb, verbose, origSize, n, ordering), - (2,data, Kmax, lamb, verbose, origSize, n, ordering), - (3,data, Kmax, lamb, verbose, origSize, n, ordering), - (4,data, Kmax, lamb, verbose, origSize, n, ordering), - (5,data, Kmax, lamb, verbose, origSize, n, ordering), - (6,data, Kmax, lamb, verbose, origSize, n, ordering), - (7,data, Kmax, lamb, verbose, origSize, n, ordering), - (8,data, Kmax, lamb, verbose, origSize, n, ordering), - (9,data, Kmax, lamb, verbose, origSize, n, ordering)]) - - #Accumulate results + returnList = pool.map(multi_run_wrapper, [(0, data, Kmax, lamb, verbose, origSize, n, ordering), + (1, data, Kmax, lamb, verbose, origSize, n, ordering), + (2, data, Kmax, lamb, verbose, origSize, n, ordering), + (3, data, Kmax, lamb, verbose, origSize, n, ordering), + (4, data, Kmax, lamb, verbose, origSize, n, ordering), + (5, data, Kmax, lamb, verbose, origSize, n, ordering), + (6, data, Kmax, lamb, verbose, origSize, n, ordering), + (7, data, Kmax, lamb, verbose, origSize, n, ordering), + (8, data, Kmax, lamb, verbose, origSize, n, ordering), + (9, data, Kmax, lamb, verbose, origSize, n, ordering)]) + + # Accumulate results for i in range(10): for j in returnList[i][0]: mseList.append(j) for j in returnList[i][1]: trainList.append(j) - #Get average of the 10 folds - plotVals = map(list, zip(*mseList)) - maxBreaks = max(plotVals[0])+1 + # Get average of the 10 folds + plotVals = list(map(list, zip(*mseList))) + maxBreaks = max(plotVals[0]) + 1 testAvg = [] for i in range(maxBreaks): num = 0 @@ -117,8 +117,8 @@ def GGSCrossVal(data, Kmax=25, lambList = [0.1, 1, 10], features = [], verbose = if (plotVals[0][j] == i): runsum = runsum + plotVals[1][j] num = num + 1 - testAvg.append(float(runsum)/num) - plotVals2 = map(list, zip(*trainList)) + testAvg.append(float(runsum) / num) + plotVals2 = list(map(list, zip(*trainList))) trainAvg = [] for i in range(maxBreaks): num = 0 @@ -127,101 +127,97 @@ def GGSCrossVal(data, Kmax=25, lambList = [0.1, 1, 10], features = [], verbose = if (plotVals[0][j] == i): runsum = runsum + plotVals2[1][j] num = num + 1 - trainAvg.append(float(runsum)/num) - - #Combine results for all lambdas into one list and return that + trainAvg.append(float(runsum) / num) + + # Combine results for all lambdas into one list and return that trainTestResults.append((lamb, (trainAvg, testAvg))) return trainTestResults - -#Find and return the means/regularized covariance of each segment for a given set of breakpoints -def GGSMeanCov(data, breakpoints, lamb, features = [], verbose = False): +# Find and return the means/regularized covariance of each segment for a given set of breakpoints +def GGSMeanCov(data, breakpoints, lamb, features=[], verbose=False): data = data.T - #Select the desired features + # Select the desired features if (features == []): features = range(data.shape[1]) - data = data[:,features] - m,n = data.shape + data = data[:, features] + m, n = data.shape numSegments = len(breakpoints) - 1 mean_covs = [] for i in range(numSegments): - - #Get mean and regularized covariance of current segment - tempData = data[breakpoints[i]:breakpoints[i+1],:] - m,n = tempData.shape + # Get mean and regularized covariance of current segment + tempData = data[breakpoints[i]:breakpoints[i + 1], :] + m, n = tempData.shape empMean = np.mean(tempData, axis=0) - empCov = np.cov(tempData.T,bias = True) - regularizedCov = empCov + float(lamb)*np.identity(n)/m + empCov = np.cov(tempData.T, bias=True) + regularizedCov = empCov + float(lamb) * np.identity(n) / m mean_covs.append((empMean, regularizedCov)) return mean_covs - - - - - - -#HELPER FUNCTIONS -def calculateLikelihood(data, breaks,lamb): +# HELPER FUNCTIONS +def calculateLikelihood(data, breaks, lamb): ll = 0 for i in range(len(breaks) - 1): - tempData = data[breaks[i]:breaks[i+1],:] - m,n = tempData.shape - empCov = np.cov(tempData.T,bias = True) - ll = ll - (m*np.linalg.slogdet(empCov + float(lamb)*np.identity(n)/m)[1] - float(lamb) * np.trace(np.linalg.inv(empCov + float(lamb)*np.identity(n)/m))) + tempData = data[breaks[i]:breaks[i + 1], :] + m, n = tempData.shape + empCov = np.cov(tempData.T, bias=True) + ll = ll - (m * np.linalg.slogdet(empCov + float(lamb) * np.identity(n) / m)[1] - float(lamb) * np.trace( + np.linalg.inv(empCov + float(lamb) * np.identity(n) / m))) return ll + def addBreak(data, lamb): - #Initialize parameters - m,n = data.shape + # Initialize parameters + m, n = data.shape origMean = np.mean(data, axis=0) - origCov = np.cov(data.T,bias = True) - origLL = m*np.linalg.slogdet(origCov + float(lamb)*np.identity(n)/m)[1] - float(lamb) * np.trace(np.linalg.inv(origCov + float(lamb)*np.identity(n)/m)) - totSum = m*(origCov+np.outer(origMean,origMean)) - muLeft = data[0,:]/n - muRight = (m * origMean - data[0,:])/(m-1) - runSum = np.outer(data[0,:],data[0,:]) - #Loop through all samples, find point where breaking the segment would have the largest LL increase - minLL = origLL + origCov = np.cov(data.T, bias=True) + origLL = m * np.linalg.slogdet(origCov + float(lamb) * np.identity(n) / m)[1] - float(lamb) * np.trace( + np.linalg.inv(origCov + float(lamb) * np.identity(n) / m)) + totSum = m * (origCov + np.outer(origMean, origMean)) + muLeft = data[0, :] / n + muRight = (m * origMean - data[0, :]) / (m - 1) + runSum = np.outer(data[0, :], data[0, :]) + # Loop through all samples, find point where breaking the segment would have the largest LL increase + minLL = origLL minInd = 0 - for i in range(2,m-1): - #Update parameters - runSum = runSum + np.outer(data[i-1,:],data[i-1,:]) - muLeft = ((i-1)*muLeft + data[i-1,:])/(i) - muRight = ((m-i+1) * muRight - data[i-1,:])/(m-i) - sigLeft = runSum/(i) - np.outer(muLeft, muLeft) - sigRight = (totSum - runSum)/(m-i) - np.outer(muRight,muRight) - - #Compute Cholesky, LogDet, and Trace - Lleft = np.linalg.cholesky(sigLeft + float(lamb)*np.identity(n)/i) - Lright = np.linalg.cholesky(sigRight + float(lamb)*np.identity(n)/(m-i)) - llLeft = 2*sum(map(math.log, np.diag(Lleft))) - llRight = 2*sum(map(math.log, np.diag(Lright))) - (trLeft, trRight) = (0,0) - if(lamb > 0): - trLeft = math.pow(np.linalg.norm(np.linalg.inv(Lleft)),2) - trRight = math.pow(np.linalg.norm(np.linalg.inv(Lright)),2) - LL = i*llLeft - float(lamb)*trLeft + (m-i)*llRight - float(lamb)*trRight - #Keep track of the best point so far - if(LL < minLL): + for i in range(2, m - 1): + # Update parameters + runSum = runSum + np.outer(data[i - 1, :], data[i - 1, :]) + muLeft = ((i - 1) * muLeft + data[i - 1, :]) / (i) + muRight = ((m - i + 1) * muRight - data[i - 1, :]) / (m - i) + sigLeft = runSum / (i) - np.outer(muLeft, muLeft) + sigRight = (totSum - runSum) / (m - i) - np.outer(muRight, muRight) + + # Compute Cholesky, LogDet, and Trace + Lleft = np.linalg.cholesky(sigLeft + float(lamb) * np.identity(n) / i) + Lright = np.linalg.cholesky(sigRight + float(lamb) * np.identity(n) / (m - i)) + llLeft = 2 * sum(map(math.log, np.diag(Lleft))) + llRight = 2 * sum(map(math.log, np.diag(Lright))) + (trLeft, trRight) = (0, 0) + if (lamb > 0): + trLeft = math.pow(np.linalg.norm(np.linalg.inv(Lleft)), 2) + trRight = math.pow(np.linalg.norm(np.linalg.inv(Lright)), 2) + LL = i * llLeft - float(lamb) * trLeft + (m - i) * llRight - float(lamb) * trRight + # Keep track of the best point so far + if (LL < minLL): minLL = LL minInd = i - #Return break, increase in LL - return (minInd,minLL-origLL) + # Return break, increase in LL + return (minInd, minLL - origLL) + -def adjustBreaks(data, breakpoints, newInd, lamb = 0, verbose = False, maxShuffles = 250): +def adjustBreaks(data, breakpoints, newInd, lamb=0, verbose=False, maxShuffles=250): bp = breakpoints[:] random.seed(0) - #Just one breakpoint, no need to adjust anything - if (len(bp) == 3): + # Just one breakpoint, no need to adjust anything + if (len(bp) == 3): return bp - #Keep track of what breakpoints have changed, so that we don't have to adjust ones which we know are constant + # Keep track of what breakpoints have changed, so that we don't have to adjust ones which we know are constant lastPass = dict() thisPass = dict() for b in bp: @@ -234,36 +230,40 @@ def adjustBreaks(data, breakpoints, newInd, lamb = 0, verbose = False, maxShuffl for b in bp: thisPass[b] = 0 switchAny = False - ordering = range(1,len(bp) - 1) + ordering = list(range(1, len(bp) - 1)) random.shuffle(ordering) for i in ordering: - #Check if we need to adjust it - if(lastPass[bp[i-1]] == 1 or lastPass[bp[i+1]] == 1 or thisPass[bp[i-1]] == 1 or thisPass[bp[i+1]] == 1): - tempData = data[bp[i-1]:bp[i+1], :] - ind, val = addBreak(tempData, lamb) - if (bp[i] != ind + bp[i-1] and val != 0): - lastPass[ind+bp[i-1]] = lastPass[bp[i]] + # Check if we need to adjust it + if (lastPass[bp[i - 1]] == 1 or lastPass[bp[i + 1]] == 1 or thisPass[bp[i - 1]] == 1 or thisPass[ + bp[i + 1]] == 1): + tempData = data[bp[i - 1]:bp[i + 1], :] + ind, val = addBreak(tempData, lamb) + if (bp[i] != ind + bp[i - 1] and val != 0): + lastPass[ind + bp[i - 1]] = lastPass[bp[i]] del lastPass[bp[i]] del thisPass[bp[i]] - thisPass[ind+bp[i-1]] = 1 + thisPass[ind + bp[i - 1]] = 1 if (verbose == True): - print "Moving", bp[i], "to", ind+bp[i-1], "length = ", tempData.shape[0], ind - bp[i] = ind + bp[i-1] + print("Moving", bp[i], "to", ind + bp[i - 1], "length = ", tempData.shape[0], ind) + bp[i] = ind + bp[i - 1] switchAny = True if (switchAny == False): return bp return bp + def multi_run_wrapper(args): return oneFold(*args) + + def oneFold(fold, data, breakpoints, lamb, verbose, origSize, n, ordering): # Remove 10% of data for test set mseList = [] trainList = [] - testSet = np.sort(ordering[(fold)*origSize/10:(fold+1)*origSize/10]) + testSet = np.sort(ordering[(fold) * int(origSize / 10):(fold + 1) * int(origSize / 10)]) mask = np.ones(origSize, dtype=bool) mask[testSet] = False - trainData = data[mask,:] + trainData = data[mask, :] # Solve for test and train error testSize = len(testSet) trainSize = origSize - testSize @@ -273,43 +273,45 @@ def oneFold(fold, data, breakpoints, lamb, verbose, origSize, n, ordering): (mse, currBreak) = (0, 1) temp = trainData[0:i[1]] empMean = np.mean(temp, axis=0) - empCov = np.cov(temp.T,bias = True) + float(lamb)*np.identity(n)/temp.shape[0] + empCov = np.cov(temp.T, bias=True) + float(lamb) * np.identity(n) / temp.shape[0] invCov = np.linalg.inv(empCov) - #Calculate test error + # Calculate test error for j in range(testSize): - #Find which break it's in + # Find which break it's in adj = testSet[j] - j - cb = max(sum(1 for k in i if k < adj),1) + cb = max(sum(1 for k in i if k < adj), 1) if (currBreak != cb): currBreak = cb - temp = trainData[i[currBreak-1]:i[currBreak]] + temp = trainData[i[currBreak - 1]:i[currBreak]] empMean = np.mean(temp, axis=0) - empCov = np.cov(temp.T,bias = True) + float(lamb)*np.identity(n)/temp.shape[0] - invCov = np.linalg.inv(empCov) - #Compute likelihood - ldet = 0.5*np.linalg.slogdet(invCov)[1] - ll = ldet - 0.5*(data[testSet[j]] - empMean).dot(invCov).dot((data[testSet[j]] - empMean)) - n*math.log(2*math.pi)/2 - mse = mse+ll - mseList.append((len(i)-2, mse/testSize)) - #Calculate training error + empCov = np.cov(temp.T, bias=True) + float(lamb) * np.identity(n) / temp.shape[0] + invCov = np.linalg.inv(empCov) + # Compute likelihood + ldet = 0.5 * np.linalg.slogdet(invCov)[1] + ll = ldet - 0.5 * (data[testSet[j]] - empMean).dot(invCov).dot((data[testSet[j]] - empMean)) - n * math.log( + 2 * math.pi) / 2 + mse = mse + ll + mseList.append((len(i) - 2, mse / testSize)) + # Calculate training error tErr = 0 currBreak = 1 temp = trainData[0:i[1]] empMean = np.mean(temp, axis=0) - empCov = np.cov(temp.T,bias = True) + float(lamb)*np.identity(n)/temp.shape[0] + empCov = np.cov(temp.T, bias=True) + float(lamb) * np.identity(n) / temp.shape[0] invCov = np.linalg.inv(empCov) - for j in range(1,trainSize): - if(j in i): + for j in range(1, trainSize): + if (j in i): currBreak = currBreak + 1 - temp = trainData[i[currBreak-1]:i[currBreak]] + temp = trainData[i[currBreak - 1]:i[currBreak]] empMean = np.mean(temp, axis=0) - empCov = np.cov(temp.T,bias = True) + float(lamb)*np.identity(n)/temp.shape[0] - invCov = np.linalg.inv(empCov) - #Compute likelihood - ldet = 0.5*np.linalg.slogdet(invCov)[1] - ll = ldet - 0.5*(trainData[j] - empMean).dot(invCov).dot((trainData[j] - empMean)) - n*math.log(2*math.pi)/2 - tErr = tErr+ll - trainList.append((len(i)-2, tErr/trainSize)) + empCov = np.cov(temp.T, bias=True) + float(lamb) * np.identity(n) / temp.shape[0] + invCov = np.linalg.inv(empCov) + # Compute likelihood + ldet = 0.5 * np.linalg.slogdet(invCov)[1] + ll = ldet - 0.5 * (trainData[j] - empMean).dot(invCov).dot((trainData[j] - empMean)) - n * math.log( + 2 * math.pi) / 2 + tErr = tErr + ll + trainList.append((len(i) - 2, tErr / trainSize)) return mseList, trainList From 672c600d8452cf28faceb5fce5cc50bae2b20e6c Mon Sep 17 00:00:00 2001 From: Shayan Bastani <32581929+shayanb2004@users.noreply.github.com> Date: Thu, 22 Sep 2022 15:29:30 +0330 Subject: [PATCH 2/2] fixed a bug for k_max selection fixed a bug where a K with value lower than K_max caused an error --- ggs.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/ggs.py b/ggs.py index 8e1d765..979d559 100644 --- a/ggs.py +++ b/ggs.py @@ -46,8 +46,13 @@ def GGS(data, Kmax, lamb, features=[], verbose=False): # Check if our algorithm is finished if (newVal == 0): print("We are done adding breakpoints!") - print(breaks) - return breaks, plotPoints + breaks = adjustBreaks(data, breaks, [newInd], lamb, verbose)[:] # new line + # Calculate likelihood #new line + ll = calculateLikelihood(data, breaks, lamb) # new line + breakPoints.append(breaks[:]) # new line + plotPoints.append(ll) # new line + print(breakPoints) + return breakPoints, plotPoints # Add new breakpoint breaks.append(newInd)