Repository: ICB-DCM/pyABC

Issue:
Allow for initial weights in AdaptiveAggregatedDistance
Use case is, if one has a reasonable set of first weights but would want to have them changed adaptively over the following populations. 
Initializing with weights will also skip the prepopulation.

Discussion:
codecov-io: # [Codecov](https://codecov.io/gh/ICB-DCM/pyABC/pull/279?src=pr&el=h1) Report
> Merging [#279](https://codecov.io/gh/ICB-DCM/pyABC/pull/279?src=pr&el=desc) into [develop](https://codecov.io/gh/ICB-DCM/pyABC/commit/f03104337cff23cfdbcdaf3c7bc31cf8a9c7ec8b?src=pr&el=desc) will **not change** coverage.
> The diff coverage is `n/a`.

[![Impacted file tree graph](https://codecov.io/gh/ICB-DCM/pyABC/pull/279/graphs/tree.svg?width=650&token=GcyGdqSygr&height=150&src=pr)](https://codecov.io/gh/ICB-DCM/pyABC/pull/279?src=pr&el=tree)

```diff
@@           Coverage Diff            @@
##           develop     #279   +/-   ##
========================================
  Coverage    83.76%   83.76%           
========================================
  Files           88       88           
  Lines         4871     4871           
========================================
  Hits          4080     4080           
  Misses         791      791
```



------

[Continue to review full report at Codecov](https://codecov.io/gh/ICB-DCM/pyABC/pull/279?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/ICB-DCM/pyABC/pull/279?src=pr&el=footer). Last update [f031043...70afb26](https://codecov.io/gh/ICB-DCM/pyABC/pull/279?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/pyabc/distance/distance.py b/pyabc/distance/distance.py
index b9cb65f..245de38 100644
--- a/pyabc/distance/distance.py
+++ b/pyabc/distance/distance.py
@@ -143,6 +143,9 @@ class AdaptivePNormDistance(PNormDistance):
 
     p: float, optional (default = 2)
         p for p-norm. Required p >= 1, p = np.inf allowed (infinity-norm).
+    initial_weights: dict, optional
+        Weights to be used in the initial iteration. Dictionary with
+        observables as keys and weights as values.
     factors: dict, optional
         As in PNormDistance.
     adaptive: bool, optional (default = True)
@@ -165,13 +168,13 @@ class AdaptivePNormDistance(PNormDistance):
         smallest non-zero absolute weight. In practice usually not necessary,
         it is theoretically required to ensure convergence.
 
-
     .. [#prangle] Prangle, Dennis. "Adapting the ABC Distance Function".
                 Bayesian Analysis, 2017. doi:10.1214/16-BA1002.
     """
 
     def __init__(self,
                  p: float = 2,
+                 initial_weights: dict = None,
                  factors: dict = None,
                  adaptive: bool = True,
                  scale_function=None,
@@ -180,6 +183,7 @@ class AdaptivePNormDistance(PNormDistance):
         # call p-norm constructor
         super().__init__(p=p, weights=None, factors=factors)
 
+        self.initial_weights = initial_weights
         self.factors = factors
         self.adaptive = adaptive
 
@@ -218,6 +222,11 @@ class AdaptivePNormDistance(PNormDistance):
         super().initialize(t, get_all_sum_stats, x_0)
         self.x_0 = x_0
 
+        # initial weights pre-defined
+        if self.initial_weights is not None:
+            self.weights[t] = self.initial_weights
+            return
+
         # execute function
         all_sum_stats = get_all_sum_stats()
 
@@ -488,6 +497,9 @@ class AdaptiveAggregatedDistance(AggregatedDistance):
 
     distances: List[Distance]
         As in AggregatedDistance.
+    initial_weights: list, optional
+        Weights to be used in the initial iteration. List with
+        a weight for each distance function.
     factors: Union[List, dict]
         As in AggregatedDistance.
     adaptive: bool, optional (default = True)
@@ -504,10 +516,12 @@ class AdaptiveAggregatedDistance(AggregatedDistance):
     def __init__(
             self,
             distances: List[Distance],
+            initial_weights: List = None,
             factors: Union[List, dict] = None,
             adaptive: bool = True,
             scale_function: Callable = None):
         super().__init__(distances=distances)
+        self.initial_weights = initial_weights
         self.factors = factors
         self.adaptive = adaptive
         self.x_0 = None
@@ -525,6 +539,10 @@ class AdaptiveAggregatedDistance(AggregatedDistance):
         super().initialize(t, get_all_sum_stats, x_0)
         self.x_0 = x_0
 
+        if self.initial_weights is not None:
+            self.weights[t] = self.initial_weights
+            return
+
         # execute function
         all_sum_stats = get_all_sum_stats()
 


Test patch:
diff --git a/test/test_distance_function.py b/test/test_distance_function.py
index d3f8f1f..e40c2fc 100644
--- a/test/test_distance_function.py
+++ b/test/test_distance_function.py
@@ -84,7 +84,7 @@ def test_pnormdistance():
     # call distance function, also to initialize w
     d = dist_f(abc.sample_from_prior()[0], abc.sample_from_prior()[1], t=0)
 
-    expected = pow(1**2 + 2**2, 1 / 2)
+    expected = pow(1**2 + 2**2, 1/2)
     assert expected == d
 
     assert sum(abs(a - b) for a, b in
@@ -117,6 +117,7 @@ def test_adaptivepnormdistance():
             scale_function=scale_function)
         dist_f.initialize(0, abc.sample_from_prior, x_0=x_0)
         dist_f(abc.sample_from_prior()[0], abc.sample_from_prior()[1], t=0)
+        assert dist_f.weights[0] != {'s1': 1, 's2': 1, 's3': 1}
 
     # test max weight ratio
     for scale_function in scale_functions:
@@ -125,6 +126,28 @@ def test_adaptivepnormdistance():
             max_weight_ratio=20)
         dist_f.initialize(0, abc.sample_from_prior, x_0=x_0)
         dist_f(abc.sample_from_prior()[0], abc.sample_from_prior()[1], t=0)
+        assert dist_f.weights[0] != {'s1': 1, 's2': 1, 's3': 1}
+
+
+def test_adaptivepnormdistance_initial_weights():
+    abc = MockABC([{'s1': -1, 's2': -1, 's3': -1},
+                   {'s1': -1, 's2': 0, 's3': 1}])
+    x_0 = {'s1': 0, 's2': 0, 's3': 1}
+
+    # first test that for PNormDistance, the weights stay constant
+    initial_weights = {'s1': 1, 's2': 2, 's3': 3}
+    dist_f = AdaptivePNormDistance(initial_weights=initial_weights)
+    dist_f.initialize(0, abc.sample_from_prior, x_0=x_0)
+    assert dist_f.weights[0] == initial_weights
+
+    # call distance function
+    d = dist_f(abc.sample_from_prior()[0], abc.sample_from_prior()[1], t=0)
+    expected = pow(sum([(2*1)**2, (3*2)**2]), 1/2)
+    assert expected == d
+
+    # check updating works
+    dist_f.update(1, abc.sample_from_prior)
+    assert dist_f.weights[1] != dist_f.weights[0]
 
 
 def test_aggregateddistance():
@@ -165,6 +188,33 @@ def test_adaptiveaggregateddistance():
         val = distance(
             abc.sample_from_prior()[0], abc.sample_from_prior()[1], t=0)
         assert isinstance(val, float)
+        assert (distance.weights[0] != [1, 1]).any()
+
+
+def test_adaptiveaggregateddistance_calibration():
+    abc = MockABC([{'s0': -1, 's1': -1},
+                   {'s0': -1, 's1': 0}])
+    x_0 = {'s0': 0, 's1': 0}
+
+    def distance0(x, x_0):
+        return abs(x['s0'] - x_0['s0'])
+
+    def distance1(x, x_0):
+        return np.sqrt((x['s1'] - x_0['s1'])**2)
+
+    scale_functions = [span, mean, median]
+    initial_weights = np.array([2, 3])
+    for scale_function in scale_functions:
+        distance = AdaptiveAggregatedDistance(
+            [distance0, distance1], scale_function=scale_function,
+            initial_weights=initial_weights)
+        distance.initialize(0, abc.sample_from_prior, x_0=x_0)
+        val = distance(
+            abc.sample_from_prior()[0], abc.sample_from_prior()[1], t=0)
+        assert isinstance(val, float)
+        assert (distance.weights[0] == initial_weights).all()
+        distance.update(1, abc.sample_from_prior)
+        assert (distance.weights[1] != distance.weights[0]).all()
 
 
 def test_normalkernel():