Slightly different metrics

discretesampling/base/algorithms/smc.py

@@ -7,6 +7,13 @@
 from discretesampling.base.algorithms.smc_components.normalisation import normalise
 from discretesampling.base.algorithms.smc_components.effective_sample_size import ess
 from discretesampling.base.algorithms.smc_components.resampling import systematic_resampling
+from discretesampling.base.algorithms.smc_components.knapsack_resampling import knapsack_resampling
+from discretesampling.base.algorithms.smc_components.minError_ImportanceResampling import min_error_continuous_state_resampling, min_error_importance_resampling
+from discretesampling.base.algorithms.smc_components.variational_resampling import  kl
+from discretesampling.base.algorithms.smc_components.importance_resampling_version3 import importance_resampling_v3
+from discretesampling.base.algorithms.smc_components.residual_resampling import residual_resampling
+
+
 
 
 class DiscreteVariableSMC():
@@ -37,26 +44,57 @@ def __init__(self, variableType, target, initialProposal, proposal=None,
         self.initialProposal = initialProposal
         self.target = target
 
-    def sample(self, Tsmc, N, seed=0, verbose=True):
+    def sample(self, Tsmc, N,a, resampling, seed=0, verbose=True):
+
         loc_n = int(N/self.exec.P)
+
         rank = self.exec.rank
         mvrs_rng = RNG(seed)
         rngs = [RNG(i + rank*loc_n + 1 + seed) for i in range(loc_n)]  # RNG for each particle
 
         initialParticles = [self.initialProposal.sample(rngs[i], self.target) for i in range(loc_n)]
         current_particles = initialParticles
-        logWeights = np.array([self.target.eval(p) - self.initialProposal.eval(p, self.target) for p in initialParticles])
+
+        logWeights = np.array([self.target.eval(p)[0] - self.initialProposal.eval(p, self.target) for p in initialParticles])
 
         display_progress_bar = verbose and rank == 0
         progress_bar = tqdm(total=Tsmc, desc="SMC sampling", disable=not display_progress_bar)
-
         for t in range(Tsmc):
+            tot_new_possibilities_for_predictions = []
             logWeights = normalise(logWeights, self.exec)
             neff = ess(logWeights, self.exec)
 
             if math.log(neff) < math.log(N) - math.log(2):
-                current_particles, logWeights = systematic_resampling(
-                    current_particles, logWeights, mvrs_rng, exec=self.exec)
+
+
+                if (resampling == "systematic"):
+                    current_particles, logWeights = systematic_resampling(
+                        current_particles, logWeights, mvrs_rng, exec=self.exec)
+
+                elif (resampling == "knapsack"):
+                    current_particles, logWeights, _ = knapsack_resampling(
+                        current_particles, np.exp(logWeights), mvrs_rng)
+
+                elif (resampling == "min_error"):
+                    current_particles, logWeights, _ = min_error_continuous_state_resampling(
+                        current_particles, np.exp(logWeights), mvrs_rng, N)
+
+                elif (resampling == "variational"):
+                    new_ancestors, logWeights = kl(logWeights)
+                    current_particles = np.array(current_particles)[new_ancestors].tolist()
+
+                elif (resampling == "min_error_imp"):
+                    current_particles, logWeights= min_error_importance_resampling(
+                        current_particles, np.exp(logWeights), mvrs_rng, N)
+
+                elif (resampling == "CIR"):
+                    current_particles, logWeights= importance_resampling_v3(
+                        current_particles, np.exp(logWeights), mvrs_rng, N)
+
+                elif (resampling == "residual"):
+                    current_particles, logWeights= residual_resampling(
+                        current_particles, np.exp(logWeights), mvrs_rng, N)
+
 
             new_particles = copy.copy(current_particles)
 
@@ -65,7 +103,7 @@ def sample(self, Tsmc, N, seed=0, verbose=True):
             # Sample new particles and calculate forward probabilities
             for i in range(loc_n):
                 forward_proposal = self.proposal
-                new_particles[i] = forward_proposal.sample(current_particles[i], rng=rngs[i])
+                new_particles[i] = forward_proposal.sample(current_particles[i], a, rng=rngs[i])
                 forward_logprob[i] = forward_proposal.eval(current_particles[i], new_particles[i])
 
             if self.use_optimal_L:
@@ -79,13 +117,15 @@ def sample(self, Tsmc, N, seed=0, verbose=True):
                     Lkernel = self.Lkernel
                     reverse_logprob = Lkernel.eval(new_particles[i], current_particles[i])
 
-                current_target_logprob = self.target.eval(current_particles[i])
-                new_target_logprob = self.target.eval(new_particles[i])
+                current_target_logprob, current_possibilities_for_predictions = self.target.eval(current_particles[i])
+                new_target_logprob, new_possibilities_for_predictions = self.target.eval(new_particles[i])
 
                 logWeights[i] += new_target_logprob - current_target_logprob + reverse_logprob - forward_logprob[i]
-
+
+                tot_new_possibilities_for_predictions.append(new_possibilities_for_predictions)
+            #if t<Tsmc:
             current_particles = new_particles
             progress_bar.update(1)
 
         progress_bar.close()
-        return current_particles
+        return current_particles,tot_new_possibilities_for_predictions, logWeights

discretesampling/base/algorithms/smc_components/importance_resampling_version3.py

@@ -0,0 +1,125 @@
+import numpy as np
+def importance_resampling_v3(x, w, mvrs_rng , N=None):
+
+    if N is None:
+        N = len(w)
+
+    #x_new = np.zeros_like(x)
+    x_new = []
+    w_new = np.zeros_like(w)
+
+    if np.sum(w) != 1:
+        w = w / np.sum(w)
+
+    S = [i for i in range(N) if w[i] < 1 / N]
+
+    NS = [i for i in range(N) if w[i] > 2 / N]
+
+    if not NS:
+        x_new, w_new, _ = systematic_resampling(x, w, mvrs_rng, N)
+        return x_new, w_new
+
+    w_tot = np.sum(w[S])
+    wwt = w_tot * N
+    noninteger = np.ceil(wwt) - wwt
+
+    WS = np.concatenate([w[S], [noninteger / N]])
+    Nsmall = int(np.sum(WS) * N)
+
+    WS = WS / np.sum(WS)
+
+    xS = []
+    for i in S:
+        xS.append(x[i])
+
+    xS.append(x[NS[0]])
+
+    x_small, w_small, _ = systematic_resampling(xS, WS, mvrs_rng, Nsmall)
+    w_small = np.exp(w_small)
+    w_small = np.ones(N) / N
+
+    for i in range(Nsmall):
+        #x_new[i] = x_small[i]
+        x_new.append(x_small[i])
+        w_new[i] = w_small[i]
+
+    NSindices = np.setdiff1d(np.arange(N), S)
+
+    wNS = w[NS[0]]
+    wNS = wNS * N - noninteger
+    ww = np.zeros_like(w)
+    for i in NSindices:
+        ww[i] = w[i]
+
+    if NS:
+        i = next(iter(NS))  # Get the first element from the set
+        ww[i] = wNS / N
+
+    ww[NSindices] = ww[NSindices] / np.sum(ww[NSindices])
+
+    xnotS=[]
+    for i in NSindices:
+        xnotS.append(x[i])
+
+    x_NS, w_NS = importance_resampling(xnotS, ww[NSindices], mvrs_rng, N - Nsmall)
+    w_NS = np.exp(w_NS)
+
+    w_NS = w_NS * (N - Nsmall) / N
+
+    x_new[Nsmall:] = x_NS
+    w_new[Nsmall:] = w_NS
+
+    log_w_new = np.log(w_new)
+
+    return x_new, log_w_new
+
+
+########################################################################
+def importance_resampling(x, w, mvrs_rng, N=None):
+
+    if N is None:
+        N = len(w)
+
+    _, _, nc = systematic_resampling(x, w, mvrs_rng, N)
+
+    quantisedweights = nc / N
+
+    x_new, _, ncopies = systematic_resampling(x, quantisedweights, mvrs_rng, N)
+
+    ww = w / quantisedweights
+    w_new = np.repeat(ww, ncopies, axis=0)
+
+    w_new /= np.sum(w_new)
+
+    log_w_new = np.log(w_new)
+
+    return x_new, log_w_new
+
+########################################################################
+
+def systematic_resampling(x, w,  mvrs_rng , N=None):
+    if N is None:
+        N = len(w)
+
+    N = int(N)
+
+    x_new = []
+    w_new = np.ones(N) / N
+    log_w_new = np.log(w_new)
+
+    cw = np.cumsum(w)
+
+    u = mvrs_rng.uniform()
+
+    ncopies = np.zeros(len(x), dtype=int)
+
+    for i in range(N):
+        j = 0
+
+        while cw[j] < (i + u) / N:
+            j += 1
+
+        x_new.append(x[j])
+        ncopies[j] += 1
+
+    return x_new, log_w_new, ncopies

discretesampling/base/algorithms/smc_components/knapsack_resampling.py

@@ -0,0 +1,83 @@
+import numpy as np
+from discretesampling.base.algorithms.smc_components.resampling import check_stability
+from discretesampling.base.executor import Executor
+
+
+def knapsack(W, wt, val, n):
+    K = [[0 for x in range(W + 1)] for x in range(n + 1)]
+
+    # Build table K[][] in bottom up manner
+    for i in range(n + 1):
+        for w in range(W + 1):
+            if i == 0 or w == 0:
+                K[i][w] = 0
+            elif wt[i - 1] <= w:
+                K[i][w] = max(val[i - 1] + K[i - 1][w - wt[i - 1]], K[i - 1][w])
+            else:
+                K[i][w] = K[i - 1][w]
+
+    res = K[n][W]
+
+    zero_one_copies = np.zeros(n).astype('i')
+
+    w = W
+    for i in range(n, 0, -1):
+        if res <= 0:
+            break
+        if res == K[i - 1][w]:
+            continue
+        else:
+            zero_one_copies[i-1] = 1
+
+            # Since this weight is included its value is deducted
+            res = res - val[i - 1]
+            w = w - wt[i - 1]
+
+    return zero_one_copies, K[n][W]
+
+
+def knapsack_resampling(x, w, mvrs_rng, exec=Executor()):
+    x_new = []
+    E = len(w)
+    T = mvrs_rng.uniform(low=0.01, high=1.0)
+    max_w = np.max(w)
+    if max_w > T:
+        T = max_w
+    C = int(E*T)
+    N = len(w)
+
+    w_knap = np.ceil(E * w).astype(int)
+    #w_knap[-1] = 53
+    value = (w_knap ** 2).astype(int)
+
+    zero_one_copies, res = knapsack(W=C, wt=w_knap, val=value, n=N)
+
+    wsum = np.sum(zero_one_copies * w)
+
+    ncopies = (N * (zero_one_copies * w)/wsum + 0.5).astype(int)
+    ncopies = check_stability(ncopies, exec)
+    #x_new = np.repeat(x, ncopies)
+    i = 0
+    for j in range(len(ncopies)):
+        for k in range(ncopies[j]):
+            i += 1
+            x_new.append(x[j])
+    log_new_w = np.repeat(np.log(1/N), N).astype('float32')
+
+    return x_new, log_new_w, ncopies
+
+# """
+# profit = [60, 100, 120]
+# weight = [10, 20, 30]
+# W = 50
+# n = len(profit)
+# zero_one_copies, res = knapsack(W, weight, profit, n)
+#
+# print(zero_one_copies)
+# print(res)
+# """
+#
+# w = np.array([1.0, 2, 3, 8, 16])/30
+# x = np.array([1.2, 2.3, 3.4, 4.5, 6.0])
+#
+# new_x, new_w, ncopies = knapsack_resampling(x, w)