diff --git a/examples/QGA.png b/examples/QGA.png
new file mode 100644
index 0000000..a10d694
Binary files /dev/null and b/examples/QGA.png differ
diff --git a/examples/comparison-proba.py b/examples/comparison-proba.py
new file mode 100755
index 0000000..e118c42
--- /dev/null
+++ b/examples/comparison-proba.py
@@ -0,0 +1,79 @@
+#!/usr/bin/env python3
+
+import copy
+
+from pyrimidine import *
+from pyrimidine.benchmarks.optimization import *
+
+# generate a knapsack problem randomly
+n_bags = 50
+evaluate = Knapsack.random(n=n_bags)
+
+
+class YourIndividual(BinaryChromosome // n_bags):
+
+    def _fitness(self):
+        return evaluate(self.decode())
+
+
+class YourPopulation(HOFPopulation):
+    element_class = YourIndividual
+    default_size = 20
+
+
+from pyrimidine.deco import add_memory
+
+@add_memory({'measure_result': None, 'fitness': None})
+class MyIndividual(QuantumChromosome // n_bags):
+
+    def _fitness(self):
+        return evaluate(self.decode())
+
+    def backup(self, check=False):
+        f = super().fitness
+        if not check or (self.memory['fitness'] is None or f > self.memory['fitness']):
+            self._memory = {
+            'measure_result': self.measure_result,
+            'fitness': f
+            }
+
+
+class MyPopulation(HOFPopulation):
+
+    element_class = MyIndividual
+    default_size = 20
+
+    def init(self):
+        for i in self:
+            i.backup()
+        super().init()
+
+    def backup(self, check=True):
+        for i in self:
+            i.backup(check=check)
+
+    def transition(self, *args, **kwargs):
+        """
+        Update the `hall_of_fame` after each step of evolution
+        """
+        super().transition(*args, **kwargs)
+        self.backup()
+
+
+stat={'Mean Fitness': 'mean_fitness', 'Best Fitness': 'best_fitness'}
+mypop = MyPopulation.random()
+
+yourpop = YourPopulation([YourIndividual(i.decode()) for i in mypop])
+mydata = mypop.evolve(n_iter=100, stat=stat, history=True)
+yourdata = yourpop.evolve(n_iter=100, stat=stat, history=True)
+
+import matplotlib.pyplot as plt
+fig = plt.figure()
+ax = fig.add_subplot(111)
+yourdata[['Mean Fitness', 'Best Fitness']].plot(ax=ax)
+mydata[['Mean Fitness', 'Best Fitness']].plot(ax=ax)
+ax.legend(('Mean Fitness', 'Best Fitness', 'Mean Fitness(Quantum)', 'Best Fitness(Quantum)'))
+ax.set_xlabel('Generations')
+ax.set_ylabel('Fitness')
+ax.set_title(f'Demo of (Quantum)GA: {n_bags}-Knapsack Problem')
+plt.show()
diff --git a/examples/example1.py b/examples/example1.py
index 4c2ec66..366d63b 100755
--- a/examples/example1.py
+++ b/examples/example1.py
@@ -33,7 +33,7 @@ def _fitness(self):
         return _evaluate(self.decode())
 
 
-MyPopulation = HOFPopulation[MyIndividual] // 12
+MyPopulation = HOFPopulation[MyIndividual] // 16
 
 pop = MyPopulation.random()
 stat = {'Mean Fitness':'mean_fitness', 'Best Fitness':'best_fitness'}
@@ -47,5 +47,5 @@ def _fitness(self):
     data[['Mean Fitness', 'Best Fitness']].plot(ax=ax)
     ax.set_xlabel('Generations')
     ax.set_ylabel('Fitness')
-    ax.set_title('Demo of GA')
+    ax.set_title('Demo for GA')
     plt.show()
diff --git a/examples/example4.py b/examples/example4.py
index 3d43140..e5aced0 100755
--- a/examples/example4.py
+++ b/examples/example4.py
@@ -19,10 +19,14 @@ class _Chromosome(BinaryChromosome):
     def decode(self):
         return c(self)
 
+from pyrimidine.deco import fitness_cache
+
+@fitness_cache
 class ExampleIndividual(PolyIndividual):
     """
     You should implement the methods, cross, mute
     """
+
     element_class = _Chromosome
     default_size = 10
 
@@ -30,13 +34,13 @@ def _fitness(self):
         x = self.decode()
         return evaluate(x)
 
-
+@fitness_cache
 class MyIndividual(ExampleIndividual, SimulatedAnnealing):
 
     def get_neighbour(self):
         cpy = self.clone()
         r = randint(0, self.n_chromosomes-1)
-        cpy.chromosomes[r].mutate()
+        cpy[r].mutate()
         return cpy
 
 
@@ -57,6 +61,7 @@ def get_neighbour(self):
     LocalSearchPopulation.element_class = MyIndividual
 
     lga = LocalSearchPopulation.random(n_individuals=20, n_chromosomes=10, size=10)
+
     lga.mate_prob = 0.9
     d= lga.evolve(n_iter=10, stat=stat, history=True)
     d[['Mean Fitness', 'Best Fitness']].plot(ax=ax, style='.-')
diff --git a/examples/example5.py b/examples/example5.py
index a5eb380..d42371f 100755
--- a/examples/example5.py
+++ b/examples/example5.py
@@ -11,6 +11,9 @@
 evaluate = Function1DApproximation(function=lambda x:10*np.arctan(x), lb=-2, ub=2, basis=_my_basis)
 n_basis = len(evaluate.basis)
 
+from pyrimidine.deco import fitness_cache
+
+@fitness_cache
 class MyIndividual(makeIndividual(FloatChromosome, n_chromosomes=1, size=n_basis)):
     def _fitness(self):
         return evaluate(self.chromosome)
@@ -19,10 +22,10 @@ def _fitness(self):
 class MyPopulation(HOFPopulation):
     element_class = MyIndividual
 
-pop = MyPopulation.random(n_individuals=200)
+pop = MyPopulation.random(n_individuals=100)
 
 stat = {'Best Fitness': 'best_fitness', 'Mean Fitness': 'mean_fitness'}
-data = pop.evolve(n_iter=500, stat=stat, history=True)
+data = pop.evolve(n_iter=200, stat=stat, history=True)
 
 
 import matplotlib.pyplot as plt
diff --git a/pyrimidine/__init__.py b/pyrimidine/__init__.py
index 6ef402c..8d99087 100755
--- a/pyrimidine/__init__.py
+++ b/pyrimidine/__init__.py
@@ -16,7 +16,7 @@
 from .ba import *
 
 
-__version__ = "1.4.2"
+__version__ = "1.5.1"
 
 __template__ = """
 from pyrimidine import MonoBinaryIndividual
diff --git a/pyrimidine/base.py b/pyrimidine/base.py
index b9dfd53..0202a47 100755
--- a/pyrimidine/base.py
+++ b/pyrimidine/base.py
@@ -64,7 +64,7 @@ class MyPopulation(SGAPopulation):
 from .meta import *
 from .mixin import *
 
-from .deco import clear_cache
+from .deco import side_effect
 
 
 class BaseGene:
@@ -133,11 +133,21 @@ def decode(self):
 
     @classmethod
     def encode(cls, x):
+        # encode x to a chromosome
         raise NotImplementedError
 
     def equal(self, other):
         return np.array_equal(self, other)
 
+    def clone(self, *args, **kwargs):
+        raise NotImplementedError
+
+    def replicate(self):
+        # Replication operation of a chromosome
+        ind = self.clone()
+        ind.mutate()
+        return ind
+
 
 class BaseIndividual(FitnessMixin, metaclass=MetaContainer):
     """Base class of individuals
@@ -149,7 +159,8 @@ class BaseIndividual(FitnessMixin, metaclass=MetaContainer):
 
     element_class = BaseChromosome
     default_size = 2
-    alias = {"chromosomes": "elements"}
+    alias = {"chromosomes": "elements",
+    "n_chromosomes": "n_elements"}
 
     def __repr__(self):
         # seperate the chromosomes with $ 
@@ -205,19 +216,23 @@ def _fitness(self):
         else:
             raise NotImplementedError
 
-    def cross(self, other, k=None):
+    def clone(self, type_=None):
+        if type_ is None:
+            type_ = self.__class__
+        return type_([c.clone(type_=type_.element_class) for c in self])
+
+    def cross(self, other):
         # Cross operation of two individual
         return self.__class__([chromosome.cross(other_c) for chromosome, other_c in zip(self.chromosomes, other.chromosomes)])
 
-    @clear_cache
+    @side_effect
     def mutate(self, copy=False):
         # Mutating operation of an individual
-        self.clear_cache()
         for chromosome in self.chromosomes:
             chromosome.mutate()
         return self
 
-    def replicate(self, k=2):
+    def replicate(self):
         # Replication operation of an individual
         ind = self.clone()
         ind.mutate()
@@ -256,7 +271,6 @@ def __mul__(self, n):
             TYPE: BasePopulation
         """
         assert isinstance(n, np.int_) and n>0, 'n must be a positive integer'
-        from .population import StandardPopulation
         C = StandardPopulation[self.__class__]
         return C([self.clone() for _ in range(n)])
 
@@ -270,7 +284,7 @@ def __rmul__(self, other):
         return self.__class__([other * this for this in self.chromosomes])
 
 
-class BasePopulation(PopulationMixin, metaclass=MetaHighContainer):
+class BasePopulation(PopulationMixin, metaclass=MetaContainer):
     """The base class of population in GA
     
     Represents a state of a stachostic process (Markov process)
@@ -356,7 +370,6 @@ def select(self, n_sel=None, tournsize=None):
         else:
             raise Exception('No winners in the selection!')
 
-    @clear_cache
     def merge(self, other, n_sel=None):
         """Merge two populations.
 
@@ -366,16 +379,16 @@ def merge(self, other, n_sel=None):
         """
 
         if isinstance(other, BasePopulation):
-            self.individuals += other.individuals
+            self.extend(other.individuals)
         elif isinstance(other, typing.Iterable):
-            self.individuals.extend(other)
+            self.extend(other)
         else:
             raise TypeError("`other` should be a population or a list/tuple of individuals")
 
         if n_sel:
             self.select(n_sel)
 
-    @clear_cache
+    @side_effect
     def mutate(self, mutate_prob=None):
         """Mutate the whole population.
 
@@ -389,7 +402,6 @@ def mutate(self, mutate_prob=None):
             if random() < (mutate_prob or self.mutate_prob):
                 individual.mutate()
 
-    @clear_cache
     def mate(self, mate_prob=None):
         """Mate the whole population.
 
@@ -402,9 +414,10 @@ def mate(self, mate_prob=None):
         mate_prob = mate_prob or self.mate_prob
         offspring = [individual.cross(other) for individual, other in zip(self.individuals[::2], self.individuals[1::2])
         if random() < mate_prob]
-        self.individuals.extend(offspring)
+        self.extend(offspring)
         self.offspring = self.__class__(offspring)
 
+    @side_effect
     def local_search(self, *args, **kwargs):
         """Call local searching method
         
@@ -454,10 +467,14 @@ def rank(self, tied=False):
             for k, i in enumerate(sorted_list):
                 i.ranking = k / self.n_individuals
 
-    @clear_cache
     def cross(self, other):
         # Cross two populations as two individuals
         k = randint(1, self.n_individuals-2)
+        return self.__class__(self[k:] + other[:k])
+
+    def migrate(self, other):
+        # migrate between two populations
+        k = randint(1, self.n_individuals-2)
         self.individuals = self[k:] + other[:k]
         other.individuals = other[k:] + self[:k]
 
@@ -470,7 +487,7 @@ class BaseMultiPopulation(PopulationMixin, metaclass=MetaHighContainer):
     
     Attributes:
         default_size (int): the number of populations
-        element_class (TYPE): type of the populations
+        element_class (TYPE): the type of the populations
         elements (TYPE): populations as the elements
         fitness (TYPE): best fitness
     """
@@ -499,8 +516,8 @@ def random(cls, n_populations=None, *args, **kwargs):
     def migrate(self, migrate_prob=None):
         for population, other in zip(self[:-1], self[1:]):
             if random() < (migrate_prob or self.migrate_prob):
-                other.individuals.append(population.best_individual.clone())
-                population.individuals.append(other.best_individual.clone())
+                other.append(population.get_best_individual(copy=True))
+                population.append(other.get_best_individual(copy=True))
 
     def transition(self, *args, **kwargs):
         super().transition(*args, **kwargs)
@@ -509,10 +526,13 @@ def transition(self, *args, **kwargs):
     def best_fitness(self):
         return max(map(attrgetter('best_fitness'), self))
 
-    def get_best_individual(self):
+    def get_best_individual(self, copy=True):
         bests = map(methodcaller('get_best_individual'), self)
         k = np.argmax([b.fitness for b in bests])
-        return bests[k]
+        if copy:
+            return bests[k].clone()
+        else:
+            return bests[k]
 
     @property
     def individuals(self):
diff --git a/pyrimidine/benchmarks/cluster.py b/pyrimidine/benchmarks/cluster.py
index 1006cfb..d904b6a 100755
--- a/pyrimidine/benchmarks/cluster.py
+++ b/pyrimidine/benchmarks/cluster.py
@@ -10,6 +10,7 @@ class KMeans:
     ERM:
     min J(c,mu) = sum_c sum_{x:c} ||x-mu_c||
     """
+
     def __init__(self, X, n_components=2):
         self.X = X
         self.n_components = n_components
@@ -21,7 +22,6 @@ def random(N, p=2):
         X = np.vstack((X1, X2))
         return KMeans(X, n_components=2)
 
-
     def __call__(self, x):
         # xi = k  iff  Xi in k-class
         cs = set(x)
diff --git a/pyrimidine/benchmarks/matrix.py b/pyrimidine/benchmarks/matrix.py
index bd9fbdb..1178115 100755
--- a/pyrimidine/benchmarks/matrix.py
+++ b/pyrimidine/benchmarks/matrix.py
@@ -7,10 +7,9 @@
 
 class NMF:
     # M ~ A diag(C) B'
+
     def __init__(self, M):
         self.M = M
-        self.norm = LA.norm(M)
-
 
     @staticmethod
     def random(N=500, p=100):
@@ -20,14 +19,44 @@ def random(N=500, p=100):
             M[k] /= s[k]
         return NMF(M=M)
 
+    def __call__(self, A, B, C=None):
+        """A: N * K
+        C: K
+        B: K * p
+        """
+
+        c = A.shape[1]
+        if C is not None:
+            for i in range(c):
+                A[:,i] *= C[i]
+        return - LA.norm(self.M - np.dot(A, B))
+
+
+class SparseMF:
+    # M ~ A B'
+
+    def __init__(self, M, threshold=0.01):
+        self.M = M
+        self.threshold = threshold
+
+    @staticmethod
+    def random(N=500, p=100):
+        M = np.random.rand(N, p) * 10
+        s = M.sum(axis=1)
+        for k in range(N):
+            M[k] /= s[k]
+        return NMF(M=M)
 
     def __call__(self, A, B, C=None):
         """A: N * K
         C: K
         B: K * p
         """
+        
         c = A.shape[1]
         if C is not None:
             for i in range(c):
-                A[:,i] = A[:,i] * C[i]
-        return - LA.norm(self.M - np.dot(A, B)) / self.norm
+                A[:,i] *= C[i]
+        T = B > threshold
+        B *= T
+        return - LA.norm(self.M - np.dot(A, B)) - np.sum(T)
diff --git a/pyrimidine/benchmarks/optimization.py b/pyrimidine/benchmarks/optimization.py
index 6fae410..2428c51 100755
--- a/pyrimidine/benchmarks/optimization.py
+++ b/pyrimidine/benchmarks/optimization.py
@@ -4,7 +4,6 @@
 from .benchmarks import Problem
 
 
-
 class Knapsack(BaseProblem):
     """Knapsack Problem
 
@@ -67,7 +66,7 @@ def __call__(self, x):
             return - 1/(1 + np.exp(-v)) * M
 
 
-class MultiKnapsack(Problem):
+class MultiKnapsack(BaseProblem):
     """Multi Choice Knapsack Problem
 
     max sum_ij cij xij
@@ -126,7 +125,7 @@ def __call__(self, x):
             return - 1/(1 + np.exp(-v)) * M
 
 
-class MLE:
+class MLE(BaseProblem):
     # max likelihood estimate
     def __init__(self, pdf, x):
         self.pdf = logpdf
@@ -142,7 +141,7 @@ def __call__(self, t):
         return np.sum([self.logpdf(xi, *t) for xi in self.x])
 
 
-class MixMLE:
+class MixMLE(MLE):
     # mix version of max likelihood estimate
     # x|k ~ pk
     def __init__(self, pdfs, x):
@@ -167,7 +166,7 @@ def __call__(self, t, a):
 
 from scipy.spatial.distance import pdist, squareform
 
-class ShortestPath:
+class ShortestPath(BaseProblem):
     def __init__(self, points):
         """TSP
         
@@ -204,7 +203,7 @@ def _heart(t, a=1):
 heart_path = CurvePath(x, y)
 
 
-class MinSpanningTree:
+class MinSpanningTree(BaseProblem):
     def __init__(self, nodes, edges=[]):
         self.nodes = nodes
         self.edges = edges
@@ -224,7 +223,8 @@ def prufer_decode(self, x):
         edges.append(tuple(Q))
         return edges
 
-class FacilityLayout(object):
+
+class FacilityLayout(BaseProblem):
     '''
     F: F
     D: D
diff --git a/pyrimidine/chromosome.py b/pyrimidine/chromosome.py
index d5f2405..6f8e5f5 100755
--- a/pyrimidine/chromosome.py
+++ b/pyrimidine/chromosome.py
@@ -7,7 +7,6 @@
 
 from .base import BaseChromosome, BaseGene
 from .gene import *
-from .deco import clear_cache
 
 
 def _asarray(out):
@@ -29,8 +28,7 @@ def __new__(cls, array=None, element_class=None):
         if array is None:
             array = []
 
-        o = np.asarray(array, dtype=element_class).view(cls)
-        return o
+        return np.asarray(array, dtype=element_class).view(cls)
 
     def __array_finalize__(self, obj):
         if obj is None:
@@ -87,12 +85,12 @@ def cross(self, other):
         k = randint(1, len(self)-1)
         return self.__class__(np.concatenate((self[:k], other[k:]), axis=0))
 
-    def clone(self, *args, **kwargs):
-        o = self.copy()
-        o.set_cache(**self._cache)
-        return o
+    def clone(self, type_=None):
+        if type_ is None:
+            return self.__class__(np.copy(self))
+        else:
+            return type_(np.copy(self))
 
-    @clear_cache
     def mutate(self, indep_prob=0.1):
         for i in range(len(self)):
             if random() < indep_prob:
@@ -106,7 +104,6 @@ class VectorChromosome(NumpyArrayChromosome):
 
 class MatrixChromosome(NumpyArrayChromosome):
     
-    @clear_cache
     def mutate(self, indep_prob=0.1):
         r, c = self.shape
         for i in range(r):
@@ -126,24 +123,26 @@ class NaturalChromosome(VectorChromosome):
 
     element_class = NaturalGene
 
-    @clear_cache
     def mutate(self, indep_prob=0.1):
         for i in range(len(self)):
             if random()< indep_prob:
                 self[i] = NaturalGene.random()
 
-    def __str__(self):
-        return "".join(map(str, self))
-
     def dual(self):
         return self.__class__(self.element_class.ub - self)
 
+class DigitChromosome(NaturalChromosome):
+
+    element_class = DigitGene
+
+    def __str__(self):
+        return "".join(map(str, self))
+
 
 class BinaryChromosome(NaturalChromosome):
 
     element_class = BinaryGene
 
-    @clear_cache
     def mutate(self, indep_prob=0.5):
         for i in range(len(self)):
             if random() < indep_prob:
@@ -171,7 +170,6 @@ def move_toward(self, other):
         r = choice(rotations(self, other))
         rotate(self, r)
 
-    @clear_cache
     def mutate(self):
         i, j = randint2(0, self.default_size-1)
         self[[i,j]] = self[[j,i]]
@@ -194,7 +192,9 @@ class FloatChromosome(NumpyArrayChromosome):
     element_class = FloatGene
     sigma = 0.05
 
-    @clear_cache
+    def __str__(self):
+        return "|".join(format(c, '.4') for c in self)
+
     def mutate(self, indep_prob=0.1, mu=0, sigma=None):
         sigma = sigma or self.sigma
         for i in range(len(self)):
@@ -325,7 +325,7 @@ def normalize(self):
 
 
 class QuantumChromosome(CircleChromosome):
-    
+
     measure_result = None
 
     def decode(self):
@@ -372,10 +372,9 @@ def cross(self, other):
         k = randint(1, len(self)-1)
         return self[:k] + other[k:]
 
-    def clone(self, *args, **kwargs):
+    def clone(self, type_=None):
         return copy.copy(self)
 
-    @clear_cache
     def mutate(self, indep_prob=0.1):
         a = self.random()
         for k in range(len(self)):
diff --git a/pyrimidine/deco.py b/pyrimidine/deco.py
new file mode 100755
index 0000000..1ace938
--- /dev/null
+++ b/pyrimidine/deco.py
@@ -0,0 +1,198 @@
+#!/usr/bin/env python3
+
+"""
+Decorators
+"""
+
+from types import MethodType
+
+import copy
+
+
+def clear_cache(func):
+    def mthd(obj, *args, **kwargs):
+        result = func(obj, *args, **kwargs)
+        obj.clear_cache()
+        return result
+    return mthd
+
+def side_effect(func):
+    """Decorator for methods with side effect
+
+    Apply the decorator to methods with side effects.
+    If all the methods called by a particular method have the decorator applied,
+    it is not advisable to include the decorator in that method.
+    
+    Args:
+        func (TYPE): a method
+    
+    Returns:
+        Decorated method
+    """
+
+    def mthd(obj, *args, **kwargs):
+        result = func(obj, *args, **kwargs)
+        if hasattr(obj, '_cache'):
+            obj.clear_cache()
+        return result
+    return mthd
+
+
+def clear_fitness(func):
+    def mthd(obj, *args, **kwargs):
+        result = func(obj, *args, **kwargs)
+        obj.clear_cache('fitness')
+        return result
+    return mthd
+
+
+class add_memory:
+
+    def __init__(self, memory={}):
+        self._memory = memory
+
+    def __call__(self, cls):
+
+        cls._memory = self._memory
+
+        def memory(obj):
+            return obj._memory
+
+        cls.memory = property(memory)
+
+        def fitness(obj):
+            if obj.memory['fitness'] is not None:
+                return obj.memory['fitness']
+            else:
+                return obj._fitness()
+
+        cls.fitness = property(fitness)
+
+        cls_clone = cls.clone
+        def _clone(obj, *args, **kwargs):
+            cpy = cls_clone(obj, *args, **kwargs)
+            cpy._memory = obj._memory
+            return cpy
+        cls.clone = _clone
+
+        cls_new = cls.__new__
+        def _new(cls, *args, **kwargs):
+            obj = cls_new(cls, *args, **kwargs)
+            obj._memory = copy.copy(cls._memory)
+            return obj
+
+        cls.__new__ = _new
+
+        return cls
+
+
+usual_side_effect = ['mutate', 'extend', 'pop', 'remove']
+
+
+class add_cache:
+
+    """Handle with cache for class
+    
+    Attributes:
+        attrs (tuple[str]): a tuple of attributes
+        methods (tuple[str]): a tuple of method names
+    """
+
+    def __init__(self, attrs, methods=()):
+        self.methods = methods
+        self.attrs = attrs
+
+    def __call__(self, cls):
+        cls._cache = {a: None for a in self.attrs}
+
+        @property
+        def cache(obj):
+            return obj._cache
+
+        cls.cache = property(cache)
+
+        def _clear_cache(obj, k=None):
+            if k is None:
+                obj._cache = {k: None for k in obj._cache.keys()}
+            elif k in obj._cache:
+                obj._cache[k] = None
+
+        def _cleared(obj, k=None):
+            if k is None:
+                return all(v == None for v in obj._cache.values())
+            elif k in obj._cache:
+                return obj._cache[k] == None
+
+        def _set_cache(obj, **d):
+            obj._cache.update(d)
+
+        cls_clone = cls.clone
+        def _clone(obj, cache=True, *args, **kwargs):
+            cpy = cls_clone(obj, *args, **kwargs)
+            if cache:
+                cpy.set_cache(**obj._cache)
+            return cpy
+
+        cls.cleared = _cleared
+        cls.clear_cache = _clear_cache
+        cls.set_cache = _set_cache
+        cls.clone = _clone
+
+        for a in self.attrs:
+            def f(obj):
+                if obj._cache[a] is None:
+                    f = getattr(obj, '_'+a)()
+                    obj._cache[a] = f
+                    return f
+                else:
+                    return obj._cache[a]
+            setattr(cls, a, property(f))
+
+        def _after_setter(obj):
+            obj.clear_cache()
+
+        cls.after_setter = _after_setter
+
+        for name in self.methods:
+            if hasattr(obj, name):
+                setattr(cls, name, clear_cache(getattr(cls, name)))
+
+        cls_new = cls.__new__
+        def _new(cls, *args, **kwargs):
+            obj = cls_new(cls, *args, **kwargs)
+            obj._cache = copy.copy(cls._cache)
+            return obj
+
+        cls.__new__ = _new
+
+        return cls
+
+
+fitness_cache = add_cache(('fitness',))
+
+
+class Regester:
+    # regerster operators, used in the future version
+
+    def __init__(name, key=None):
+        self.name = name
+        self.key = key
+
+    def __call__(cls):
+
+        def _regester_operator(self, name, key=None):
+            if hasattr(self, name):
+                raise AttributeError(f'"{name}" is an attribute of "{self.__name__}", and would not be regestered.')
+            self._operators.append(key)
+            setattr(self, name, MethodType(key, self))
+
+        def _element_regester(self, name, e):
+            if hasattr(self, e):
+                raise AttributeError(f'"{e}" is an attribute of "{self.__name__}", would not be regestered.')
+            self._elements.append(e)
+            setattr(self, name, e)
+
+        cls.regester_operator = _regester_operator
+        cls.regester_element = _regester_element
+
+        return cls
diff --git a/pyrimidine/gene.py b/pyrimidine/gene.py
index 593b779..4287ebb 100755
--- a/pyrimidine/gene.py
+++ b/pyrimidine/gene.py
@@ -1,5 +1,9 @@
 #!/usr/bin/env python3
 
+"""
+Gene classes
+"""
+
 import numpy as np
 from . import BaseGene
 
@@ -11,6 +15,9 @@ class NaturalGene(np.int_, BaseGene):
     def random(cls, *args, **kwargs):
         return np.random.randint(cls.ub, dtype=cls, *args, **kwargs)
 
+class DigitGene(NaturalGene):
+    pass
+
 class IntegerGene(np.int_, BaseGene):
     lb, ub = -10, 10
 
@@ -49,7 +56,8 @@ def period(self):
 
 
 class CircleGene(PeriodicGene):
-    lb, ub = -np.pi, np.pi
+    lb, ub = 0, 2*np.pi
+    period = 2*np.pi
 
     
 class SemiCircleGene(CircleGene):
diff --git a/pyrimidine/gsa.py b/pyrimidine/gsa.py
index dd478d2..e11da16 100644
--- a/pyrimidine/gsa.py
+++ b/pyrimidine/gsa.py
@@ -24,7 +24,6 @@ class Particle(PolyIndividual):
     default_size = 2
     accelerate = 0
 
-
     @property
     def position(self):
         return self.chromosomes[0]
@@ -32,7 +31,7 @@ def position(self):
     @position.setter
     def position(self, x):
         self.chromosomes[0] = x
-        self.__fitness = None
+        self._cache['fitness'] = None
 
     @property
     def velocity(self):
@@ -51,14 +50,14 @@ def move(self):
         D = cpy.fitness - self.fitness
         if flag:
             self.chromosomes = cpy.chromosomes
-            self.__fitness = cpy.fitness
+            self._cache['fitness'] = cpy.fitness
 
 
 class GravitySearch(PopulationMixin):
     """Standard GSA
     
     Extends:
-        BaseFitnessMixin
+        PopulationMixin
     """
 
     element_class = Particle
@@ -74,7 +73,6 @@ def compute_mass(self):
         m = (fitnesses - worst_fitness + epsilon) / (best_fitness - worst_fitness + epsilon)
         return m / m.sum()
 
-
     def compute_accelerate(self):
         # compute force
         D = np.array([[pj.position - pi.position for pi in self] for pj in self])
@@ -92,9 +90,8 @@ def compute_accelerate(self):
 
         # set accelerate
         for i, particle in enumerate(self):
-            particle.accelerate = A[i,:]
+            particle.accelerate = A[i, :]
 
-    
     def transition(self, k):
         """
         Transitation of the states of particles
diff --git a/pyrimidine/individual.py b/pyrimidine/individual.py
index ca8c7ec..a1a7a7e 100755
--- a/pyrimidine/individual.py
+++ b/pyrimidine/individual.py
@@ -85,17 +85,14 @@ def gender(self):
         raise NotImplementedError
 
 
+from .deco import add_memory, add_cache
+
+@add_memory(memory={"fitness": None})
 class MemoryIndividual(BaseIndividual):
     # Individual with memory, used in PSO
 
-    _memory = {"fitness": None}
-
-    @property
-    def memory(self):
-        return self._memory
-
     def init(self, fitness=True, type_=None):
-        self._memory = {k: None for k in self.__class__.memory.keys()}
+        self._memory = {k: None for k in self.__class__._memory.keys()}
         self.backup(check=False)
 
     def backup(self, check=False):
@@ -105,7 +102,7 @@ def backup(self, check=False):
             check (bool, optional): check whether the fitness increases.
         """
 
-        f = self._fitness()
+        f = super().fitness
         if not check or (self.memory['fitness'] is None or f > self.memory['fitness']):
             def _map(k):
                 if k == 'fitness':
@@ -118,31 +115,22 @@ def _map(k):
                     return getattr(self, k)
             self._memory = {k: _map(k) for k in self.memory.keys()}
 
-    @property
-    def fitness(self):
-        if 'fitness' in self.memory and self.memory['fitness'] is not None:
-            return self.memory['fitness']
-        else:
-            return super().fitness
-
-    def clone(self, *args, **kwargs):
-        cpy = super().clone(*args, **kwargs)
-        cpy._memory = self.memory
-        return cpy
-
 
+@add_cache(('fitness',))
 class PhantomIndividual(BaseIndividual):
     # Another implementation of the individual class with memory
 
+    _cache = {'fitness': None}
+
     phantom = None
 
-    def init(self, fitness=True, type_=None):
-        self.phantom = self.clone(fitness=fitness, type_=type_)
+    def init(self):
+        self.phantom = self.clone()
 
     def backup(self):
         if self.fitness < self.phantom.fitness:
             self.chromosomes = self.phantom.chromosomes
-            self.cache_fitness(self.phantom.fitness)
+            self.set_cache(fitness=self.phantom.fitness)
 
 
 # Following are functions to create individuals
diff --git a/pyrimidine/local_search/simulated_annealing.py b/pyrimidine/local_search/simulated_annealing.py
index 51be3ae..e5aebf4 100755
--- a/pyrimidine/local_search/simulated_annealing.py
+++ b/pyrimidine/local_search/simulated_annealing.py
@@ -32,7 +32,7 @@ class SimulatedAnnealing(PhantomIndividual):
 
     def init(self):
         # initialize phantom solution
-        self.phantom = self.clone(fitness=None)
+        self.phantom = self.clone()
 
     def transition(self, *args, **kwargs):
         T = self.initT
diff --git a/pyrimidine/meta.py b/pyrimidine/meta.py
index 8a6b977..db5fdd9 100755
--- a/pyrimidine/meta.py
+++ b/pyrimidine/meta.py
@@ -49,9 +49,7 @@ def __new__(cls, name, bases=(), attrs={}):
         attrs = inherit(attrs, 'params', bases)
 
         def _getattr(self, key):
-            if key in self.__dict__:
-                return self.__dict__[key]
-            elif key in self.params:
+            if key in self.params:
                 return self.params[key]
             elif key in self.alias:
                 return getattr(self, self.alias[key])
@@ -110,6 +108,11 @@ class cls(self, other):
             pass
         return cls
 
+    def __call__(self, *args, **kwargs):
+        obj = super().__call__(*args, **kwargs)
+        obj.params = copy.deepcopy(self.params)
+        return obj
+
 
 class MetaContainer(ParamType):
     """Meta class of containers
@@ -171,6 +174,7 @@ def _getitem(self, k):
         def _setitem(self, k, v):
             # print(DeprecationWarning('get item directly is not recommended now.'))
             self.__elements[k] =v
+            self.after_setter()
 
         def _iter(self):
             return iter(self.__elements)
@@ -340,7 +344,7 @@ def __getitem__(self, class_):
 
     def __call__(self, *args, **kwargs):
         o = super().__call__(*args, **kwargs)
-        for e in o.__elements:  # consider in future
+        for e in o:  # consider in future
             if not isinstance(e, self.element_class):
                 raise TypeError(f'"{e}" is not an instance of type "{self.element_class}"')
         return o
@@ -362,7 +366,7 @@ def __floordiv__(self, n):
 
     def __call__(self, *args, **kwargs):
         o = super().__call__(*args, **kwargs)
-        for e, t in (o.__elements, self.element_class):  # consider in future
+        for e, t in zip(o, self.element_class):  # consider in future
             if not isinstance(e, t):
                 raise TypeError(f'"{e}" is not an instance of type "{t}"')
         return o
@@ -452,9 +456,8 @@ def __floordiv__(self, n):
         return self.set(default_size=n)
 
     def __call__(self, *args, **kwargs):
-        o = super().__call__(*args, **kwargs)
-        o.params = copy.deepcopy(self.params)
+        obj = super().__call__(*args, **kwargs)
         if hasattr(self, '_cache'):
-            o._cache = copy.copy(self._cache)
-        return o
+            obj._cache = copy.copy(self._cache)
+        return obj
 
diff --git a/pyrimidine/mixin.py b/pyrimidine/mixin.py
index 2d45445..14b43a3 100755
--- a/pyrimidine/mixin.py
+++ b/pyrimidine/mixin.py
@@ -26,27 +26,13 @@
 from ezstat import Statistics
 from .errors import *
 
-from .deco import clear_cache
+from .deco import side_effect
 
 
 class IterativeMixin:
     # Mixin class for iterative algrithms
 
     params = {'n_iter': 100}
-    _cache = {}
-
-    @property
-    def cache(self):
-        return self._cache
-
-    def clear_cache(self, k=None):
-        if k is None:
-            self._cache = {k: None for k in self._cache.keys()}
-        elif k in self._cache:
-            self._cache[k] = None
-
-    def set_cache(self, **d):
-        self._cache.update(d)
 
     @property
     def solution(self):
@@ -167,6 +153,9 @@ def perf(self, n_repeats=10, *args, **kwargs):
     def clone(self, type_=None, *args, **kwargs):
         raise NotImplementedError
 
+    def copy(self):
+        raise NotImplementedError
+
     def encode(self):
         raise NotImplementedError
  
@@ -190,6 +179,10 @@ def load(filename='model.pkl'):
         else:
             raise FileNotFoundError(f'Could not find the file {filename}!')
 
+    def after_setter(self):
+        if hasattr(self, '_cache'):
+            self.clear_cache()
+
 
 class FitnessMixin(IterativeMixin):
     """Iterative models drived by the fitness/objective function
@@ -197,21 +190,9 @@ class FitnessMixin(IterativeMixin):
     The fitness should be stored until the the state of the model is changed.
     
     Extends:
-        BaseIterativeMixin
+        IterativeMixin
     """
 
-    _cache = {'fitness': None}
-
-    def cache_fitness(self, v):
-        self._cache['fitness'] = v
-
-    @property
-    def fitness(self):
-        if self._cache['fitness'] is None:
-            f = self._fitness()
-            self.cache_fitness(f)
-        return self._cache['fitness']
-
     def get_fitness(self):
         raise NotImplementedError
 
@@ -219,6 +200,10 @@ def _fitness(self):
         # the alias of the fitness
         return self.get_fitness()
 
+    @property
+    def fitness(self):
+        return self._fitness()
+
     @classmethod
     def set_fitness(cls, f=None):
         if f is None:
@@ -226,22 +211,8 @@ def set_fitness(cls, f=None):
                 f = globals()['_fitness']
             else:
                 raise Exception('Function `_fitness` is not defined before setting fitness. You may forget to create the class in the context of environment.')
-        class C(cls):
-            def _fitness(self):
-                return f(self)
-        return C
-
-    def clone(self, type_=None, fitness=True):
-        if type_ is None:
-            type_ = self.__class__
-        if fitness is True:
-            fitness = self.fitness
-        cpy = type_(list(map(methodcaller('clone', type_=type_.element_class, fitness=True), self)))
-        if fitness is True:
-            cpy.cache_fitness(self.fitness)
-        else:
-            cpy.cache_fitness(fitness)
-        return cpy
+        cls._fitness = f
+        return cls
 
     def evolve(self, stat=None, attrs=('solution',), *args, **kwargs):
         """Get the history of solution and its fitness by default.
@@ -251,9 +222,9 @@ def evolve(self, stat=None, attrs=('solution',), *args, **kwargs):
             stat = {'Fitness':'fitness'}
         return super().evolve(stat=stat, attrs=attrs, *args, **kwargs)
 
-    def after_setter(self):
-        # clean up the fitness after updating the chromosome
-        self.clear_cache()
+    @property
+    def solution(self):
+        return self
 
 
 class ContainerMixin(IterativeMixin):
@@ -264,26 +235,25 @@ def init(self, *args, **kwargs):
         for element in self:
             element.init(*args, **kwargs)
 
-    @clear_cache
     def transition(self, *args, **kwargs):
         for element in self:
             element.transition(*args, **kwargs)
 
-    @clear_cache
+    @side_effect
     def remove(self, individual):
         self.elements.remove(individual)
 
-    @clear_cache
+    @side_effect
     def pop(self, k=-1):
         self.elements.pop(k)
 
-    @clear_cache
+    @side_effect
     def extend(self, inds):
         self.elements.extend(inds)
 
-    @clear_cache
-    def add_individuals(self, inds):
-        self.elements.extend(inds)
+    @side_effect
+    def append(self, ind):
+        self.elements.append(ind)
 
 
 class PopulationMixin(FitnessMixin, ContainerMixin):
@@ -301,9 +271,6 @@ def evolve(self, stat=None, *args, **kwargs):
             'STD Fitness': 'std_fitness', 'Population': 'n_elements'}
         return super().evolve(stat=stat, *args, **kwargs)
 
-    def after_setter(self):
-        self.clear_cache()
-
     @classmethod
     def set_fitness(cls, *args, **kwargs):
         # set fitness for the element_class.
@@ -426,3 +393,13 @@ def drop(self, n=1):
             n = int(n)
         ks = self.argsort()
         self.elements = [self[k] for k in ks[n:]]
+
+    def clone(self, type_=None, fitness=True):
+        if type_ is None:
+            type_ = self.__class__
+        cpy = type_(list(map(methodcaller('clone', type_=type_.element_class, fitness=True), self)))
+        if fitness is True:
+            cpy.cache_fitness(self.fitness)
+        else:
+            cpy.cache_fitness(fitness)
+        return cpy
diff --git a/pyrimidine/multipopulation.py b/pyrimidine/multipopulation.py
index 1abebe8..2850a38 100755
--- a/pyrimidine/multipopulation.py
+++ b/pyrimidine/multipopulation.py
@@ -59,16 +59,16 @@ def _target(male, female):
         for p in ps:
             p.join()
 
-        self.populations[0].add_individuals(children[::2])
-        self.populations[1].add_individuals(children[1::2])
+        self.populations[0].extend(children[::2])
+        self.populations[1].extend(children[1::2])
 
     def match(self, male, female):
         return True
 
     def transition(self, *args, **kwargs):
         elder = self.__class__([
-            self.populations[0].clone_best_individuals(self.n_elders * self.populations[0].default_size),
-            self.populations[1].clone_best_individuals(self.n_elders * self.populations[1].default_size)
+            self.populations[0].get_best_individuals(self.n_elders * self.populations[0].default_size, copy=True),
+            self.populations[1].get_best_individuals(self.n_elders * self.populations[1].default_size, copy=True)
             ])
         self.select()
         self.mate()
diff --git a/pyrimidine/population.py b/pyrimidine/population.py
index 4404b79..d72873a 100755
--- a/pyrimidine/population.py
+++ b/pyrimidine/population.py
@@ -16,7 +16,7 @@
 from .chromosome import BinaryChromosome
 from .meta import MetaList
 
-from .deco import clear_cache
+from .deco import side_effect
 
 
 class StandardPopulation(BasePopulation):
@@ -69,7 +69,7 @@ def transition(self, *args, **kwargs):
 
         super().transition(*args, **kwargs)
         self.update_hall_of_fame()
-        self.add_individuals(self.hall_of_fame)
+        self.extend(self.hall_of_fame)
 
     def update_hall_of_fame(self):
         for ind in self:
@@ -165,7 +165,6 @@ def transition(self, *args, **kwargs):
         self.eliminate()
         self.merge(elder)
 
-    @clear_cache
     def eliminate(self):
         # remove some individuals randomly from the population
         for individual in self:
@@ -180,7 +179,6 @@ def transition(self, *args, **kwargs):
             individual.age += 1
         super().transition(*args, **kwargs)
 
-    @clear_cache
     def eliminate(self):
         # remove some old individuals
         for individual in self:
@@ -192,6 +190,12 @@ class LocalSearchPopulation(StandardPopulation):
     """Population with `local_search` method
     """
 
+    params = {'n_local_iter': 10}
+
+    def init(self):
+        for i in self:
+            i.n_iter = self.n_local_iter
+
     def transition(self, *args, **kwargs):
         """Transitation of the states of population
 
@@ -205,7 +209,6 @@ class ModifiedPopulation(StandardPopulation):
 
     params = {'mutate_prob_ub':0.5, 'mutate_prob_lb':0.1}
 
-    @clear_cache
     def mutate(self):
         fm = self.best_fitness
         fa = self.mean_fitness
@@ -219,11 +222,21 @@ def mutate(self):
                 individual.mutate()
 
 
-def makeBinaryPopulation(cls=None, n_populations=20, size=8, as_chromosome=True):
-    # make a binary population
+def makeBinaryPopulation(n_individuals=20, size=8, as_chromosome=True, cls=None):
+    """Make a binary population
+    
+    Args:
+        n_individuals (int, optional): the number of the individuals in the population
+        size (int, optional): the size of an individual or a chromosome
+        as_chromosome (bool, optional): take chromosomes as the individuals of the population
+        cls (None, optional): the type of the population (HOFPopulation by default)
+    
+    Returns:
+        subclass of BasePopulation
+    """
 
     cls = cls or HOFPopulation
     if as_chromosome:
-        return cls[BinaryChromosome // size] // n_populations
+        return cls[BinaryChromosome // size] // n_individuals
     else:
-        return cls[binaryIndividual(size)] // n_populations
+        return cls[binaryIndividual(size)] // n_individuals
diff --git a/pyrimidine/saga.py b/pyrimidine/saga.py
index 5f7554c..fd33f88 100755
--- a/pyrimidine/saga.py
+++ b/pyrimidine/saga.py
@@ -126,7 +126,7 @@ def mate(self):
             if random() < min(individual.cross_prob, other.cross_prob):
                 if self.match(individual, other):
                     children.append(individual.cross(other))
-        self.add_individuals(children)
+        self.extend(children)
 
     @classmethod
     def match(cls, individual, other):
diff --git a/tests/test_de.py b/tests/test_de.py
index 1bf3a15..15cb847 100644
--- a/tests/test_de.py
+++ b/tests/test_de.py
@@ -1,6 +1,5 @@
 #!/usr/bin/env python3
 
-import unittest
 
 from pyrimidine.individual import MonoIndividual
 from pyrimidine.population import HOFPopulation, BasePopulation
@@ -9,41 +8,50 @@
 
 from pyrimidine.benchmarks.special import rosenbrock
 
+import pytest
 
-class TestDE(unittest.TestCase):
+@pytest.fixture(scope="class")
+def pop():
+    n = 20
+    f = rosenbrock
 
-    def setUp(self):
+    class MyIndividual(MonoIndividual):
+        element_class = FloatChromosome // n
 
-        n = 20
-        f = rosenbrock
+        def _fitness(self):
+            return -f(self.chromosome)
 
-        class MyIndividual(MonoIndividual):
-            element_class = FloatChromosome // n
+    class _Population1(DifferentialEvolution, BasePopulation):
+        element_class = MyIndividual
+        default_size = 10
 
-            def _fitness(self):
-                return -f(self.chromosome)
+    class _Population2(HOFPopulation):
+        element_class = MyIndividual
+        default_size = 10
 
-        class _Population1(DifferentialEvolution, BasePopulation):
-            element_class = MyIndividual
-            default_size = 10
+    # _Population2 = HOFPopulation[MyIndividual] // 10
 
-        class _Population2(HOFPopulation):
-            element_class = MyIndividual
-            default_size = 10
+    return _Population1, _Population2
 
-        # _Population2 = HOFPopulation[MyIndividual] // 10
-        self.Population1 = _Population1
-        self.population1 = Population1.random()
-        self.Population2 = _Population2
-        self.population2 = Population2.random()
 
+class TestDE:
+
+    @classmethod
+    def setup_class(cls):
+        cls.Populations = cls.request.getfixturevalue("pop")
+    
     def test_clone(self):
-        self.population2 = self.population1.clone(type_=self.Population2) # population 2 with the same initial values to population 1
-        assert isinstance(self.population2, self.Population2)
+        P1, P2 = cls.Populations
+        p2 = P1.random().clone(type_=P2)
+        assert isinstance(p2, P2)
 
     def test_evolve(self):
+        P1, P2 = cls.Populations
+        self.population1 = P1.random()
+        self.population2 = P2.random()
         stat={'Mean Fitness':'mean_fitness', 'Best Fitness':'best_fitness'}
-        data1 = self.population1.evolve(stat=stat, n_iter=10, history=True)
-        data2 = self.population2.evolve(stat=stat, n_iter=10, history=True)
-        assert True
+        data1 = self.population1.evolve(stat=stat, n_iter=5, history=True)
+        data2 = self.population2.evolve(stat=stat, n_iter=5, history=True)
+        assert ('Mean Fitness' in data1.columns and 'Best Fitness' in data1.columns and
+            'Mean Fitness' in data2.columns and 'Best Fitness' in data2.columns)
 
diff --git a/tests/test_deco.py b/tests/test_deco.py
new file mode 100644
index 0000000..b96eab0
--- /dev/null
+++ b/tests/test_deco.py
@@ -0,0 +1,28 @@
+#!/usr/bin/env python3
+
+
+from pyrimidine.deco import *
+
+
+class TestDeco:
+
+    def test_cache(self):
+
+        @add_cache(attrs=('f',), methods=('bar',))
+        class C:
+            v = 1
+            def _f(self):
+                return self.v
+
+            def bar(self):
+                pass
+
+        c = C()
+
+        assert c._f() == 1
+        assert c.f == 1
+        assert c._cache['f'] == 1
+        c.v = 2
+        assert c._f() == 2
+        c.bar()
+        assert c.f == 2
diff --git a/tests/test_pso.py b/tests/test_pso.py
index b06cdd1..3da072c 100644
--- a/tests/test_pso.py
+++ b/tests/test_pso.py
@@ -6,12 +6,14 @@
 from pyrimidine.pso import Particle, ParticleSwarm
 from pyrimidine.benchmarks.special import rosenbrock
 
+
 def evaluate(x):
     return - rosenbrock(x)
 
-class TestPSO(unittest.TestCase):
 
-    def setUp(self):
+class TestPSO():
+
+    def test_pso(self):
         # generate a knapsack problem randomly
 
         class _Particle(Particle):
@@ -25,8 +27,9 @@ class MyParticleSwarm(ParticleSwarm, BasePopulation):
             element_class = _Particle
             default_size = 10
 
-        self.ParticleSwarm = MyParticleSwarm
+        ParticleSwarm = MyParticleSwarm
+
+        pop = ParticleSwarm.random()
+        pop.transition()
+        assert isinstance(pop, ParticleSwarm)
 
-    def test_pso(self):
-        pop = self.ParticleSwarm.random()
-        data = pop.transition()
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 2f97494..48ab936 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -4,7 +4,7 @@
 from pyrimidine.utils import *
 
 
-class TestUtils(unittest.TestCase):
+class TestUtils:
 
     def test_pattern(self):
         assert pattern([[0,1,1],[1,1,0]]) == '*1*'
diff --git a/travis.yml b/travis.yml
new file mode 100644
index 0000000..2af55be
--- /dev/null
+++ b/travis.yml
@@ -0,0 +1,20 @@
+language: python
+
+python:
+  - "3.8"
+  - "3.9"
+  - "3.10"
+  - "3.11"
+
+install: pip install numpy
+  
+before_script: 
+  - pip install pytest
+  - pip install poetry
+
+script:
+  - poetry build
+  - pytest test.py
+
+notifications:
+  email: false