High Order Nostationary Fuzzy Time Series - HONSFTS

pyFTS/nonstationary/common.py

@@ -124,18 +124,18 @@ class FuzzySet(FS.FuzzySet):
     def get_midpoint(self, t):
 
         self.perturbate_parameters(t)
+        param = self.perturbated_parameters[t]
 
         if self.mf == Membership.gaussmf:
-            return self.perturbated_parameters[t][0]
+            return param[0]
         elif self.mf == Membership.sigmf:
-            return self.perturbated_parameters[t][1]
+            return param[1]
         elif self.mf == Membership.trimf:
-            return self.perturbated_parameters[t][1]
+            return param[1]
         elif self.mf == Membership.trapmf:
-            param = self.perturbated_parameters[t]
             return (param[2] - param[1]) / 2
         else:
-            return self.perturbated_parameters[t]
+            return param
 
     def get_lower(self, t):
 

pyFTS/nonstationary/flrg.py

@@ -13,7 +13,7 @@ class NonStationaryFLRG(flrg.FLRG):
     def get_membership(self, data, t, window_size=1):
         ret = 0.0
         if isinstance(self.LHS, (list, set)):
-            assert len(self.LHS) == len(data)
+            #assert len(self.LHS) == len(data)
 
             ret = min([self.LHS[ct].membership(dat, common.window_index(t - (self.order - ct), window_size))
                        for ct, dat in enumerate(data)])
@@ -31,20 +31,20 @@ class NonStationaryFLRG(flrg.FLRG):
         else:
             return self.LHS[-1].get_midpoint(common.window_index(t, window_size))
 
-
     def get_lower(self, t, window_size=1):
-        if self.lower is None:
+        if len(self.RHS) > 0:
-            if len(self.RHS) > 0:
+            if isinstance(self.RHS, (list, set)):
-                self.lower = min([r.get_lower(common.window_index(t, window_size)) for r in self.RHS])
+                return min([r.get_lower(common.window_index(t, window_size)) for r in self.RHS])
-            else:
+            elif isinstance(self.RHS, dict):
-                self.lower = self.LHS[-1].get_lower(common.window_index(t, window_size))
+                return min([self.RHS[r].get_lower(common.window_index(t, window_size)) for r in self.RHS.keys()])
-
+        else:
-        return self.lower
+            return self.LHS[-1].get_lower(common.window_index(t, window_size))
 
     def get_upper(self, t, window_size=1):
-        if self.upper is None:
+        if len(self.RHS) > 0:
-            if len(self.RHS) > 0:
+            if isinstance(self.RHS, (list, set)):
-                self.upper = min([r.get_upper(common.window_index(t, window_size)) for r in self.RHS])
+                return max([r.get_upper(common.window_index(t, window_size)) for r in self.RHS])
-            else:
+            elif isinstance(self.RHS, dict):
-                self.upper = self.LHS[-1].get_upper(common.window_index(t, window_size))
+                return max([self.RHS[r].get_upper(common.window_index(t, window_size)) for r in self.RHS.keys()])
-        return self.upper
+        else:
+            return self.LHS[-1].get_upper(common.window_index(t, window_size))

pyFTS/nonstationary/honsfts.py

@@ -91,6 +91,8 @@ class HighOrderNonStationaryFTS(hofts.HighOrderFTS):
                 for st in rhs:
                     flrgs[flrg.strLHS()].appendRHS(st)
 
+        # flrgs = sorted(flrgs, key=lambda flrg: flrg.get_midpoint(0, window_size=1))
+
         return flrgs
 
     def train(self, data, sets=None, order=2, parameters=None):
@@ -108,6 +110,65 @@ class HighOrderNonStationaryFTS(hofts.HighOrderFTS):
         window_size = parameters if parameters is not None else 1
         self.flrgs = self.generate_flrg(ndata, window_size=window_size)
 
+    def _affected_flrgs(self, sample, k, time_displacement, window_size):
+        # print("input: " + str(ndata[k]))
+
+        affected_flrgs = []
+        affected_flrgs_memberships = []
+
+        lags = {}
+
+        for ct, dat in enumerate(sample):
+            tdisp = common.window_index((k + time_displacement) - (self.order - ct), window_size)
+            sel = [ct for ct, set in enumerate(self.sets) if set.membership(dat, tdisp) > 0.0]
+
+            if len(sel) == 0:
+                sel.append(common.check_bounds_index(dat, self.sets, tdisp))
+
+            lags[ct] = sel
+
+        # Build the tree with all possible paths
+
+        root = tree.FLRGTreeNode(None)
+
+        self.build_tree(root, lags, 0)
+
+        # Trace the possible paths and build the PFLRG's
+
+        for p in root.paths():
+            path = list(reversed(list(filter(None.__ne__, p))))
+            flrg = HighOrderNonStationaryFLRG(self.order)
+
+            for kk in path:
+                flrg.appendLHS(self.sets[kk])
+
+            affected_flrgs.append(flrg)
+            # affected_flrgs_memberships.append(flrg.get_membership(sample, disp))
+
+            #                print(flrg.strLHS())
+
+            # the FLRG is here because of the bounds verification
+            mv = []
+            for ct, dat in enumerate(sample):
+                td = common.window_index((k + time_displacement) - (self.order - ct), window_size)
+                tmp = flrg.LHS[ct].membership(dat, td)
+                # print('td',td)
+                # print('dat',dat)
+                # print(flrg.LHS[ct].name, flrg.LHS[ct].perturbated_parameters[td])
+                # print(tmp)
+
+                if (tmp == 0.0 and flrg.LHS[ct].name == self.sets[0].name and dat < self.sets[0].get_lower(td)) \
+                        or (tmp == 0.0 and flrg.LHS[ct].name == self.sets[-1].name and dat > self.sets[-1].get_upper(
+                            td)):
+                    mv.append(1.0)
+                else:
+                    mv.append(tmp)
+            # print(mv)
+
+            affected_flrgs_memberships.append(np.prod(mv))
+
+        return [affected_flrgs, affected_flrgs_memberships]
+
     def forecast(self, data, **kwargs):
 
         time_displacement = kwargs.get("time_displacement",0)
@@ -122,54 +183,32 @@ class HighOrderNonStationaryFTS(hofts.HighOrderFTS):
 
         for k in np.arange(self.order, l+1):
 
-            #print("input: " + str(ndata[k]))
-
-            disp = common.window_index(k + time_displacement, window_size)
-
-            affected_flrgs = []
-            affected_flrgs_memberships = []
-
-            lags = {}
-
             sample = ndata[k - self.order: k]
 
-            for ct, dat in enumerate(sample):
+            affected_flrgs, affected_flrgs_memberships = self._affected_flrgs(sample, k,
-                tdisp = common.window_index((k + time_displacement) - (self.order - ct), window_size)
+                                                                              time_displacement, window_size)
-                sel = [ct for ct, set in enumerate(self.sets) if set.membership(dat, tdisp) > 0.0]
 
-                if len(sel) == 0:
+            #print([str(k) for k in affected_flrgs])
-                    sel.append(common.check_bounds_index(dat, self.sets, tdisp))
+            #print(affected_flrgs_memberships)
-
-                lags[ct] = sel
-
-            # Build the tree with all possible paths
-
-            root = tree.FLRGTreeNode(None)
-
-            self.build_tree(root, lags, 0)
-
-            # Trace the possible paths and build the PFLRG's
-
-            for p in root.paths():
-                path = list(reversed(list(filter(None.__ne__, p))))
-                flrg = HighOrderNonStationaryFLRG(self.order)
-
-                for kk in path:
-                    flrg.appendLHS(self.sets[kk])
-
-                affected_flrgs.append(flrg)
-                affected_flrgs_memberships.append(flrg.get_membership(ndata[k - self.order: k], disp))
-
-            #print(affected_sets)
 
             tmp = []
-            for ct, aset in enumerate(affected_flrgs):
+            tdisp = common.window_index(k + time_displacement, window_size)
-                if aset.strLHS() in self.flrgs:
+            if len(affected_flrgs) == 0:
-                    tmp.append(self.flrgs[aset.strLHS()].get_midpoint(tdisp) *
+                tmp.append(common.check_bounds(sample[-1], self.sets, tdisp))
-                               affected_flrgs_memberships[ct])
+            elif len(affected_flrgs) == 1:
+                if affected_flrgs[0].strLHS() in self.flrgs:
+                    flrg = affected_flrgs[0]
+                    tmp.append(self.flrgs[flrg.strLHS()].get_midpoint(tdisp))
                 else:
-                    tmp.append(aset.LHS[-1].get_midpoint(tdisp))
+                    tmp.append(flrg.LHS[-1].get_midpoint(tdisp))
-
+            else:
+                for ct, aset in enumerate(affected_flrgs):
+                    if aset.strLHS() in self.flrgs:
+                        tmp.append(self.flrgs[aset.strLHS()].get_midpoint(tdisp) *
+                                   affected_flrgs_memberships[ct])
+                    else:
+                        tmp.append(aset.LHS[-1].get_midpoint(tdisp)*
+                                   affected_flrgs_memberships[ct])
             pto = sum(tmp)
 
             #print(pto)
@@ -182,7 +221,9 @@ class HighOrderNonStationaryFTS(hofts.HighOrderFTS):
 
     def forecastInterval(self, data, **kwargs):
 
-        time_displacement = kwargs.get("time_displacement",0)
+        time_displacement = kwargs.get("time_displacement", 0)
+
+        window_size = kwargs.get("window_size", 1)
 
         ndata = np.array(self.doTransformations(data))
 
@@ -190,21 +231,48 @@ class HighOrderNonStationaryFTS(hofts.HighOrderFTS):
 
         ret = []
 
-        for k in np.arange(0, l):
+        for k in np.arange(self.order, l + 1):
 
-            tdisp = k + time_displacement
+            sample = ndata[k - self.order: k]
 
-            affected_sets = [ [set.name, set.membership(ndata[k], tdisp)]
+            affected_flrgs, affected_flrgs_memberships = self._affected_flrgs(sample, k,
-                              for set in self.sets if set.membership(ndata[k], tdisp) > 0.0]
+                                                                              time_displacement, window_size)
+
+            # print([str(k) for k in affected_flrgs])
+            # print(affected_flrgs_memberships)
 
             upper = []
             lower = []
-            for aset in affected_sets:
+
-                lower.append(self.flrgs[aset[0]].get_lower(tdisp) * aset[1])
+            tdisp = common.window_index(k + time_displacement, window_size)
-                upper.append(self.flrgs[aset[0]].get_upper(tdisp) * aset[1])
+            if len(affected_flrgs) == 0:
+                aset = common.check_bounds(sample[-1], self.sets, tdisp)
+                lower.append(aset.get_lower(tdisp))
+                upper.append(aset.get_upper(tdisp))
+            elif len(affected_flrgs) == 1:
+                if affected_flrgs[0].strLHS() in self.flrgs:
+                    flrg = affected_flrgs[0]
+                    lower.append(self.flrgs[flrg.strLHS()].get_lower(tdisp))
+                    upper.append(self.flrgs[flrg.strLHS()].get_upper(tdisp))
+                else:
+                    lower.append(flrg.LHS[-1].get_lower(tdisp))
+                    upper.append(flrg.LHS[-1].get_upper(tdisp))
+            else:
+                for ct, aset in enumerate(affected_flrgs):
+                    if aset.strLHS() in self.flrgs:
+                        lower.append(self.flrgs[aset.strLHS()].get_lower(tdisp) *
+                                   affected_flrgs_memberships[ct])
+                        upper.append(self.flrgs[aset.strLHS()].get_upper(tdisp) *
+                                   affected_flrgs_memberships[ct])
+                    else:
+                        lower.append(aset.LHS[-1].get_lower(tdisp) *
+                                   affected_flrgs_memberships[ct])
+                        upper.append(aset.LHS[-1].get_upper(tdisp) *
+                                   affected_flrgs_memberships[ct])
 
             ret.append([sum(lower), sum(upper)])
 
+
         ret = self.doInverseTransformations(ret, params=[data[self.order - 1:]])
 
         return ret

pyFTS/nonstationary/nsfts.py

@@ -96,12 +96,18 @@ class NonStationaryFTS(fts.FTS):
             tmp = []
 
             if len(affected_sets) == 1 and self.method == 'fuzzy':
-                tmp.append(affected_sets[0][0].get_midpoint(tdisp))
+                aset = affected_sets[0][0]
+                if aset.name in self.flrgs:
+                    tmp.append(self.flrgs[aset.name].get_midpoint(tdisp))
+                else:
+                    tmp.append(aset.get_midpoint(tdisp))
             else:
                 for aset in affected_sets:
                     if self.method == 'fuzzy':
                         if aset[0].name in self.flrgs:
                             tmp.append(self.flrgs[aset[0].name].get_midpoint(tdisp) * aset[1])
+                        else:
+                            tmp.append(aset[0].get_midpoint(tdisp) * aset[1])
                     elif self.method == 'maximum':
                         if aset.name in self.flrgs:
                             tmp.append(self.flrgs[aset.name].get_midpoint(tdisp))
@@ -120,7 +126,7 @@ class NonStationaryFTS(fts.FTS):
 
     def forecastInterval(self, data, **kwargs):
 
-        time_displacement = kwargs.get("time_displacement",0)
+        time_displacement = kwargs.get("time_displacement", 0)
 
         window_size = kwargs.get("window_size", 1)
 
@@ -132,16 +138,46 @@ class NonStationaryFTS(fts.FTS):
 
         for k in np.arange(0, l):
 
+            # print("input: " + str(ndata[k]))
+
             tdisp = common.window_index(k + time_displacement, window_size)
 
-            affected_sets = [ [set.name, set.membership(ndata[k], tdisp)]
+            if self.method == 'fuzzy':
-                              for set in self.sets if set.membership(ndata[k], tdisp) > 0.0]
+                affected_sets = [[set, set.membership(ndata[k], tdisp)]
+                                 for set in self.sets if set.membership(ndata[k], tdisp) > 0.0]
+            elif self.method == 'maximum':
+                mv = [set.membership(ndata[k], tdisp) for set in self.sets]
+                ix = np.ravel(np.argwhere(mv == max(mv)))
+                affected_sets = [self.sets[x] for x in ix]
+
+            if len(affected_sets) == 0:
+                if self.method == 'fuzzy':
+                    affected_sets.append([common.check_bounds(ndata[k], self.sets, tdisp), 1.0])
+                else:
+                    affected_sets.append(common.check_bounds(ndata[k], self.sets, tdisp))
 
             upper = []
             lower = []
-            for aset in affected_sets:
+
-                lower.append(self.flrgs[aset[0]].get_lower(tdisp) * aset[1])
+            if len(affected_sets) == 1:
-                upper.append(self.flrgs[aset[0]].get_upper(tdisp) * aset[1])
+                #print(2)
+                aset = affected_sets[0][0]
+                if aset.name in self.flrgs:
+                    lower.append(self.flrgs[aset.name].get_lower(tdisp))
+                    upper.append(self.flrgs[aset.name].get_upper(tdisp))
+                else:
+                    lower.append(aset.get_lower(tdisp))
+                    upper.append(aset.get_upper(tdisp))
+            else:
+                for aset in affected_sets:
+                    #print(aset)
+                    if aset[0].name in self.flrgs:
+                        lower.append(self.flrgs[aset[0].name].get_lower(tdisp) * aset[1])
+                        upper.append(self.flrgs[aset[0].name].get_upper(tdisp) * aset[1])
+                    else:
+                        lower.append(aset[0].get_lower(tdisp) * aset[1])
+                        upper.append(aset[0].get_upper(tdisp) * aset[1])
+
 
             ret.append([sum(lower), sum(upper)])
 

pyFTS/pwfts.py

@@ -60,7 +60,7 @@ class ProbabilisticWeightedFLRG(hofts.HighOrderFLRG):
         for count, set in enumerate(self.LHS):
             mv.append(set.membership(x[count]))
 
-        min_mv = np.prod(mv)
+        min_mv = np.min(mv)
         return min_mv
 
     def partition_function(self, uod, nbins=100):
@@ -73,6 +73,7 @@ class ProbabilisticWeightedFLRG(hofts.HighOrderFLRG):
         return self.Z
 
     def get_midpoint(self):
+        '''Return the expectation of the PWFLRG, the weighted sum'''
         return sum(np.array([self.get_RHSprobability(s) * self.RHS[s].centroid
                              for s in self.RHS.keys()]))
 
@@ -495,6 +496,9 @@ class ProbabilisticWeightedFTS(ifts.IntervalFTS):
 
     def forecastDistribution(self, data, **kwargs):
 
+        if not isinstance(data, (list, set, np.ndarray)):
+            data = [data]
+
         smooth = kwargs.get("smooth", "none")
         nbins = kwargs.get("num_bins", 100)
 

pyFTS/tests/nonstationary.py

@@ -56,24 +56,28 @@ ws=12
 trainp = passengers[:ts]
 testp = passengers[ts:]
 
-tmp_fsp = Grid.GridPartitioner(trainp[:ws], 15)
+tmp_fsp = Grid.GridPartitioner(trainp[:50], 10)
+
 
 fsp = common.PolynomialNonStationaryPartitioner(trainp, tmp_fsp, window_size=ws, degree=1)
 
 
-#nsftsp = honsfts.HighOrderNonStationaryFTS("", partitioner=fsp)
+nsftsp = honsfts.HighOrderNonStationaryFTS("", partitioner=fsp)
-nsftsp = nsfts.NonStationaryFTS("", partitioner=fsp, method='fuzzy')
+#nsftsp = nsfts.NonStationaryFTS("", partitioner=fsp, method='fuzzy')
 
-#nsftsp.train(trainp, order=1, parameters=ws)
+nsftsp.train(trainp, order=2, parameters=ws)
 
-print(fsp)
+#print(fsp)
 
 #print(nsftsp)
 
-#tmpp = nsftsp.forecast(passengers[55:65], time_displacement=55, window_size=ws)
+tmpp = nsftsp.forecast(passengers[101:104], time_displacement=101, window_size=ws)
+tmpi = nsftsp.forecastInterval(passengers[101:104], time_displacement=101, window_size=ws)
 
-#print(passengers[100:120])
+#print(passengers[101:104])
-#print(tmpp)
+print([k[0] for k in tmpi])
+print(tmpp)
+print([k[1] for k in tmpi])
 
 #util.plot_sets(fsp.sets,tam=[10, 5], start=0, end=100, step=2, data=passengers[:100],
 #               window_size=ws, only_lines=False)

pyFTS/tests/pwfts.py

@@ -54,20 +54,28 @@ pfts1.shortname = "1st Order"
 
 #print(pfts1_enrollments)
 
-tmp = pfts1.forecast(data[3000:3020])
+#tmp = pfts1.forecast(data[3000:3020])
 
-tmp = pfts1.forecastInterval(data[3000:3020])
+#tmp = pfts1.forecastInterval(data[3000:3020])
 
-tmp = pfts1.forecastAheadInterval(data[3000:3020],20)
+tmp = pfts1.forecastDistribution(data[3500])
 
-tmp = pfts1.forecastAheadDistribution(data[3000:3020],20, method=3, h=0.45, kernel="gaussian")
+p = 0
-print(tmp[0])
+for b in tmp[0].bins:
+    p += tmp[0].density(b)
+
+print(p)
+
+#tmp = pfts1.forecastAheadInterval(data[3000:3020],20)
+
+#tmp = pfts1.forecastAheadDistribution(data[3000:3020],20, method=3, h=0.45, kernel="gaussian")
+#print(tmp[0])
 
 #print(tmp[0].quantile([0.05, 0.95]))
 
 #pfts1_enrollments.AprioriPDF
 #norm = pfts1_enrollments.global_frequency_count
-#uod = pfts1_enrollments.get_UoD()
+#uod = pfts1.get_UoD()
 
 #for k in sorted(pfts1_enrollments.flrgs.keys())
 #    flrg = pfts1_enrollments.flrgs[k]