diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index f9e62d9d..176f163e 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -7,7 +7,7 @@ Release Notes
 0.12 Series
 ...........
 
-0.12.19 (2026-07-27)
+0.12.19 (2026-07-28)
 --------------------
 
 General:
@@ -24,12 +24,16 @@ Visualization:
   per-generation walltimes and no longer include the idle time that passed
   between a stored analysis and a later resumed run. Resolves #667.
 
-Storage:
+Storage (breaking):
 
 * Store the per-generation walltime in the database (new ``wall_time`` column,
   database version 2). Databases created with older pyABC versions must be
   migrated via ``abc-migrate`` before they can be resumed; for such databases
   the walltime plots fall back to the previous, end-time-based behavior.
+* Particle weights are now stored using the global normalization (weights sum
+  to 1 across all particles of all models, matching the in-memory
+  ``Population`` representation). This is bundled into database version 2 and
+  handled by ``abc-migrate`` for existing databases. Resolves #47.
 
 0.12.18 (2026-04-14)
 --------------------
diff --git a/README.md b/README.md
index a5647e90..5f769d7c 100644
--- a/README.md
+++ b/README.md
@@ -26,4 +26,5 @@
 - 📄 **Cite**: [https://pyabc.rtfd.io/en/latest/cite.html](https://pyabc.rtfd.io/en/latest/cite.html)
 
 #### Related Projects
-- 🧠 **Neural Posterior Estimation**: [BayesFlow](https://bayesflow.org/main/_examples/From_ABC_to_BayesFlow.html)
+- **Parameter Estimation with Likelihoods**: [pyPESTO](https://github.com/ICB-DCM/pyPESTO)
+- **Neural Posterior Estimation**: [BayesFlow](https://bayesflow.org/main/_examples/From_ABC_to_BayesFlow.html)
diff --git a/pyabc/acceptor/pdf_norm.py b/pyabc/acceptor/pdf_norm.py
index f9afbf29..be3c9607 100644
--- a/pyabc/acceptor/pdf_norm.py
+++ b/pyabc/acceptor/pdf_norm.py
@@ -32,7 +32,10 @@ def pdf_norm_max_found(
         prev_pdf_norm = -np.inf
 
     # take maximum over all normalizations
-    pdf_norm = max(prev_pdf_norm, *pdfs)
+    if len(pdfs) == 0:
+        pdf_norm = prev_pdf_norm
+    else:
+        pdf_norm = max(prev_pdf_norm, float(pdfs.max()))
 
     return pdf_norm
 
diff --git a/pyabc/distance/distance.py b/pyabc/distance/distance.py
index 18b08c44..14f55f12 100644
--- a/pyabc/distance/distance.py
+++ b/pyabc/distance/distance.py
@@ -105,7 +105,9 @@ def _dict_to_vect(self, x):
 
     def _calculate_whitening_transformation_matrix(self, sum_stats):
         # create data matrix, shape (n_sample, n_y)
-        x = np.asarray([self._dict_to_vect(x) for x in sum_stats])
+        # force float dtype so in-place centering below works also for
+        # integer-valued summary statistics (e.g. counts)
+        x = np.asarray([self._dict_to_vect(x) for x in sum_stats], dtype=float)
         # center
         mean = np.mean(x, axis=0)
         x -= mean
diff --git a/pyabc/distance/pnorm.py b/pyabc/distance/pnorm.py
index 9276d50d..c761dad3 100644
--- a/pyabc/distance/pnorm.py
+++ b/pyabc/distance/pnorm.py
@@ -148,7 +148,7 @@ def get_weights(self, t: int) -> np.ndarray:
 
     @staticmethod
     def format_dict(
-        vals: dict[str, float] | dict[int, dict[str, float]],
+        vals: dict[str, float] | dict[int, dict[str, float]] | None,
         t: int,
         s_ids: list[str],
     ) -> dict[int, float | np.ndarray]:
@@ -172,8 +172,7 @@ def format_dict(
             vals = {t: vals}
 
         # convert dicts to arrays
-        for _t, dct in vals.items():
-            vals[_t] = dict2arr(dct, keys=s_ids)
+        vals = {_t: dict2arr(dct, keys=s_ids) for _t, dct in vals.items()}
 
         return vals
 
@@ -322,7 +321,9 @@ def __init__(
         # call p-norm constructor
         super().__init__(p=p, fixed_weights=fixed_weights, sumstat=sumstat)
 
-        self.initial_scale_weights: dict[str, float] = initial_scale_weights
+        self.initial_scale_weights: dict[str, float] | None = (
+            initial_scale_weights
+        )
 
         self.scale_weights: dict[int, np.ndarray] = {}
 
diff --git a/pyabc/distance/scale.py b/pyabc/distance/scale.py
index ef4b1ca2..13760b19 100644
--- a/pyabc/distance/scale.py
+++ b/pyabc/distance/scale.py
@@ -141,7 +141,7 @@ def root_mean_square_deviation(
     rmse = np.sqrt(mse)
 
     # debugging
-    warn_obs_off(off_ixs=np.flatnonzero(bs > 2 * std), s_ids=s_ids)
+    warn_obs_off(off_ixs=np.flatnonzero(np.abs(bs) > 2 * std), s_ids=s_ids)
 
     return rmse
 
@@ -161,7 +161,7 @@ def std_or_rmsd(
     bs = bias(samples=samples, s0=s0)
     std = standard_deviation(samples=samples)
 
-    if sum(bs > 2 * std) > 1 / 3 * len(std):
+    if sum(np.abs(bs) > 2 * std) > 1 / 3 * len(std):
         logger.info('Too many high-bias values, correcting only for scale.')
         return std
 
@@ -169,7 +169,7 @@ def std_or_rmsd(
     rmse = np.sqrt(mse)
 
     # debugging
-    warn_obs_off(off_ixs=np.flatnonzero(bs > 2 * std), s_ids=s_ids)
+    warn_obs_off(off_ixs=np.flatnonzero(np.abs(bs) > 2 * std), s_ids=s_ids)
 
     return rmse
 
diff --git a/pyabc/epsilon/temperature.py b/pyabc/epsilon/temperature.py
index d28167fb..7e0e84ba 100644
--- a/pyabc/epsilon/temperature.py
+++ b/pyabc/epsilon/temperature.py
@@ -86,6 +86,9 @@ def __init__(
         enforce_exact_final_temperature: bool = True,
         log_file: str | None = None,
     ):
+        # normalize a single callable to a list, as all consumers iterate
+        if schemes is not None and callable(schemes):
+            schemes = [schemes]
         self.schemes = schemes
 
         if aggregate_fun is None:
diff --git a/pyabc/inference_util/inference_util.py b/pyabc/inference_util/inference_util.py
index a1df6862..f255c980 100644
--- a/pyabc/inference_util/inference_util.py
+++ b/pyabc/inference_util/inference_util.py
@@ -690,10 +690,10 @@ def create_analysis_id():
     return str(uuid.uuid4())
 
 
-def eps_from_hist(history: History, t: int | None = None) -> float:
+def eps_from_hist(history: History, t: int | None = None) -> float | None:
     """Read epsilon value for time `t` from `history`. Defaults to latest."""
     pops = history.get_all_populations()
-    if len(pops) == 0 or (t is not None and t not in pops.t):
+    if len(pops) == 0 or (t is not None and t not in pops.t.values):
         return None
     if t is None:
         return pops.epsilon.to_numpy()[-1]
diff --git a/pyabc/predictor/predictor.py b/pyabc/predictor/predictor.py
index 897c2f27..1e718396 100644
--- a/pyabc/predictor/predictor.py
+++ b/pyabc/predictor/predictor.py
@@ -732,6 +732,8 @@ def __init__(
 
         if f_score is None:
             self.f_score = root_mean_square_error
+        else:
+            self.f_score = f_score
 
         # holds the chosen predictor model
         self.chosen_one: Predictor | None = None
diff --git a/pyabc/sampler/multicore.py b/pyabc/sampler/multicore.py
index 9b2c3d9f..44660299 100644
--- a/pyabc/sampler/multicore.py
+++ b/pyabc/sampler/multicore.py
@@ -36,7 +36,10 @@ def work(feed_q, result_q, simulate_one, max_eval, single_core_sampler):
             break
 
         res = single_core_sampler.sample_until_n_accepted(
-            1, simulate_one, max_eval
+            1,
+            simulate_one,
+            t=0,  # t is not used in this sampler
+            max_eval=max_eval,
         )
         result_q.put((res, single_core_sampler.nr_evaluations_))
 
diff --git a/pyabc/sampler/redis_eps/cli.py b/pyabc/sampler/redis_eps/cli.py
index a48739cc..a82cf324 100644
--- a/pyabc/sampler/redis_eps/cli.py
+++ b/pyabc/sampler/redis_eps/cli.py
@@ -199,7 +199,8 @@ def _work(
             logger.info('Received stop signal. Shutdown redis worker.')
             return
 
-        # TODO other messages (some integers?) are ignored
+        # Any other messages on the channel (e.g. redis subscription-count
+        # notifications) are ignored.
 
         # check total time condition
         elapsed_time = time() - start_time
diff --git a/pyabc/sge/db.py b/pyabc/sge/db.py
index 1463f608..b8d6db1a 100644
--- a/pyabc/sge/db.py
+++ b/pyabc/sge/db.py
@@ -49,11 +49,9 @@ def wait_for_job(self, ID, max_run_time_h):
         Return true if we should still wait for the job.
         Return false otherwise
         """
-        # TODO Possible SQL injection error should be fixed, e.g. via
-        #  pre-calculated expressions
         with self.connection:
             results = self.connection.execute(
-                'SELECT status, time from status WHERE ID=' + str(ID)
+                'SELECT status, time from status WHERE ID=?', (ID,)
             ).fetchall()
             nr_rows = len(results)
 
diff --git a/pyabc/storage/bytes_storage.py b/pyabc/storage/bytes_storage.py
index e4d006b3..5f95b4b6 100644
--- a/pyabc/storage/bytes_storage.py
+++ b/pyabc/storage/bytes_storage.py
@@ -16,7 +16,7 @@ def r_to_py(object_):
                 py_object_ = conv.rpy2py(object_)
             # Ensure factor columns are converted to strings
             for col in py_object_.columns:
-                if isinstance(py_object_[col], pd.CategoricalDtype):
+                if isinstance(py_object_[col].dtype, pd.CategoricalDtype):
                     py_object_[col] = py_object_[col].astype(str)
             return py_object_
         return object_
diff --git a/pyabc/storage/db_model.py b/pyabc/storage/db_model.py
index 1b58d324..d5eae1bb 100644
--- a/pyabc/storage/db_model.py
+++ b/pyabc/storage/db_model.py
@@ -169,6 +169,8 @@ class SummaryStatistic(Base):
     value = Column(BytesStorage)
 
 
-def datetime2str(datetime: datetime.datetime) -> str:
+def datetime2str(datetime: datetime.datetime | None) -> str:
     """Format print datetime."""
+    if datetime is None:
+        return 'None'
     return datetime.strftime('%Y-%m-%d %H:%M:%S')
diff --git a/pyabc/storage/history.py b/pyabc/storage/history.py
index d08d95a1..3a0024a1 100644
--- a/pyabc/storage/history.py
+++ b/pyabc/storage/history.py
@@ -344,10 +344,10 @@ def get_distribution(
         ).sort_index()
         w = df[['id', 'w']].drop_duplicates().set_index('id').sort_index()
         w_arr = w.w.values
-        if w_arr.size > 0 and not np.isclose(w_arr.sum(), 1):
-            raise AssertionError(
-                f'Weight not close to 1, w.sum()={w_arr.sum()}'
-            )
+        # Stored weights are global (summing to 1 across all models); within a
+        # single model they sum to the model probability.
+        if w_arr.size > 0:
+            w_arr = w_arr / w_arr.sum()
         return pars, w_arr
 
     @with_session
@@ -735,18 +735,16 @@ def _save_to_population_db(
             # append model
             population.models.append(model)
 
-            # TODO This normalization is different than in the in-memory
-            #  population. It would be cleaner to update the db too.
-            total_model_weight = sum(p.weight for p in model_population)
-
             # iterate over model population of particles
             for py_particle in model_population:
                 # a store_item is a Particle
                 py_parameter = py_particle.parameter
 
-                # create new particle
+                # create new particle. The stored weight is the global
+                # particle weight (normalized across all particles of all
+                # models so that they sum to 1)
                 particle = Particle(
-                    w=py_particle.weight / total_model_weight,
+                    w=py_particle.weight,
                     proposal_id=py_particle.proposal_id,
                 )
                 # append particle to model
@@ -858,14 +856,13 @@ def get_model_probabilities(self, t: int | None = None) -> pd.DataFrame:
             .order_by(Model.m)
             .all()
         )
-        # TODO this is a mess
+        # Two return shapes: for a single t, a frame indexed by model id `m`
+        # with a `p` column; for t=None, a t-by-model pivot table (below).
         if t is not None:
             p_models_df = pd.DataFrame(
                 [p[:2] for p in p_models], columns=['p', 'm']
             ).set_index('m')
-            # TODO the following line is redundant
-            # only models with no-zero weight are stored for each population
-            p_models_df = p_models_df[p_models_df.p >= 0]
+            # Only models with non-zero weight are stored per population
             return p_models_df
         else:
             p_models_df = (
@@ -933,7 +930,7 @@ def get_weighted_distances(self, t: int | None = None) -> pd.DataFrame:
         distances = []
         for model in models:
             for particle in model.particles:
-                weight = particle.w * model.p_model
+                weight = particle.w
                 for sample in particle.samples:
                     weights.append(weight)
                     distances.append(sample.distance)
@@ -1018,14 +1015,31 @@ def get_weighted_sum_stats_for_model(
             .filter(ABCSMC.id == self.id)
             .filter(Population.t == t)
             .filter(Model.m == m)
+            .options(
+                subqueryload(Particle.samples).subqueryload(
+                    Sample.summary_statistics
+                )
+            )
             .all()
         )
 
+        # model probability, used to renormalize the stored global weights
+        # back to the within-model posterior (weights summing to 1)
+        p_model = (
+            self._session.query(Model.p_model)
+            .join(Population)
+            .join(ABCSMC)
+            .filter(ABCSMC.id == self.id)
+            .filter(Population.t == t)
+            .filter(Model.m == m)
+            .scalar()
+        )
+
         results = []
         weights = []
         for particle in particles:
             for sample in particle.samples:
-                weights.append(particle.w)
+                weights.append(particle.w / p_model if p_model else particle.w)
                 sum_stats = {}
                 for ss in sample.summary_statistics:
                     sum_stats[ss.name] = ss.value
@@ -1075,7 +1089,7 @@ def get_weighted_sum_stats(
 
         for model in models:
             for particle in model.particles:
-                weight = particle.w * model.p_model
+                weight = particle.w
                 for sample in particle.samples:
                     # extract sum stats
                     sum_stats = {}
@@ -1125,7 +1139,7 @@ def get_population(self, t: int | None = None):
             py_m = model.m
             for particle in model.particles:
                 # weight
-                py_weight = particle.w * model.p_model
+                py_weight = particle.w
 
                 # parameter
                 py_parameter = {}
@@ -1134,7 +1148,9 @@ def get_population(self, t: int | None = None):
                 py_parameter = PyParameter(**py_parameter)
 
                 # simulations
-                # TODO this is legacy from when there were multiple
+                # NOTE: samples is a one-to-many relationship for legacy
+                # reasons (a particle could once store multiple samples);
+                # today exactly one is expected.
                 if len(particle.samples) != 1:
                     raise AssertionError('There should be exactly one sample.')
                 sample = particle.samples[0]
diff --git a/pyabc/storage/json.py b/pyabc/storage/json.py
index d2f91fdd..e90afcae 100644
--- a/pyabc/storage/json.py
+++ b/pyabc/storage/json.py
@@ -19,7 +19,7 @@ def save_dict_to_json(dct: dict, file_: str):
     for key, val in dct.items():
         # cannot handle ndarrays
         if isinstance(val, np.ndarray):
-            dct[key] = list(val)
+            dct[key] = val.tolist()
     with open(file_, 'w') as f:
         json.dump(dct, f)
 
diff --git a/pyabc/storage/migrate.py b/pyabc/storage/migrate.py
index a5ddda71..1a22ced6 100644
--- a/pyabc/storage/migrate.py
+++ b/pyabc/storage/migrate.py
@@ -14,6 +14,57 @@
 SQLITE_STR = 'sqlite:///'
 
 
+def _to_db_file(db: str) -> str:
+    """Normalize a database identifier to a file name.
+
+    Parameters
+    ----------
+    db: Database file name, or sqlite URL ``sqlite:///<file>``.
+
+    Returns
+    -------
+    db_file: The database file name.
+
+    Raises
+    ------
+    ValueError: If a URL of a dialect other than sqlite is passed.
+    """
+    if db.startswith(SQLITE_STR):
+        return db[len(SQLITE_STR) :]
+    if '://' in db:
+        raise ValueError(
+            f'Cannot handle database identifier {db}: migration currently '
+            f'only supports sqlite databases, i.e. either a file name, or a '
+            f'URL of the form {SQLITE_STR}<file>.'
+        )
+    return db
+
+
+def _alembic_config(db: str) -> 'Config':
+    """Create the alembic configuration operating on a database.
+
+    Parameters
+    ----------
+    db: Database the migrations are applied to, either as a file name or as a
+        sqlite URL ``sqlite:///<file>``.
+
+    Returns
+    -------
+    cfg: The alembic configuration.
+    """
+    # config base path
+    base_path = os.path.dirname(os.path.abspath(__file__))
+    # read configuration file
+    cfg = Config(os.path.join(base_path, 'alembic.ini'))
+    # set absolute script location path
+    cfg.set_main_option(
+        'script_location', os.path.join(base_path, 'migrations')
+    )
+    # set target database file
+    cfg.set_main_option('sqlalchemy.url', SQLITE_STR + _to_db_file(db))
+    return cfg
+
+
 @click.command(
     help='**Migrate pyABC database**\n\n'
     "Sometimes, changes to pyABC's storage format are unavoidable. "
@@ -36,8 +87,8 @@ def migrate(src: str, dst: str, version: str) -> None:
 
     Parameters
     ----------
-    src: Source
-    dst: Destination
+    src: Source, either a file name or a sqlite URL
+    dst: Destination, either a file name or a sqlite URL
     version: Version to migrate to
     """
     if Config is None or command is None:
@@ -48,10 +99,11 @@ def migrate(src: str, dst: str, version: str) -> None:
         return
 
     # to file paths if URLs
-    if src.startswith(SQLITE_STR):
-        src = src[len(SQLITE_STR) :]
-    if dst.startswith(SQLITE_STR):
-        dst = dst[len(SQLITE_STR) :]
+    try:
+        src, dst = _to_db_file(src), _to_db_file(dst)
+    except ValueError as e:
+        print(f'Error: {e}')
+        return
 
     # copy file
     if src != dst:
@@ -61,16 +113,5 @@ def migrate(src: str, dst: str, version: str) -> None:
         # copy source to destination
         shutil.copyfile(src=src, dst=dst)
 
-    # config base path
-    base_path = os.path.dirname(os.path.abspath(__file__))
-    # read configuration file
-    cfg = Config(os.path.join(base_path, 'alembic.ini'))
-    # set absolute script location path
-    cfg.set_main_option(
-        'script_location', os.path.join(base_path, 'migrations')
-    )
-    # set target database file
-    cfg.set_main_option('sqlalchemy.url', SQLITE_STR + dst)
-
     # run the actual upgrade
-    command.upgrade(cfg, version)
+    command.upgrade(_alembic_config(dst), version)
diff --git a/pyabc/storage/migrations/versions/2_20260724_add_populations_wall_time.py b/pyabc/storage/migrations/versions/2_20260724_add_populations_wall_time.py
index 0f3d5af6..c618ee19 100644
--- a/pyabc/storage/migrations/versions/2_20260724_add_populations_wall_time.py
+++ b/pyabc/storage/migrations/versions/2_20260724_add_populations_wall_time.py
@@ -1,9 +1,17 @@
-"""add populations.wall_time
+"""add populations.wall_time and store global particle weights
 
 Revision ID: 2
 Revises: 1
 Create Date: 2026-07-24 00:00:00.000000
 
+This revision bundles two v2 storage-format changes:
+
+* adds the ``populations.wall_time`` column (wall-time tracking), and
+* converts particle weights from the old per-model normalization
+  (``w = g_i / p_model``, summing to 1 within each model) to the global
+  convention (``w = g_i``, summing to 1 across all particles of all models),
+  matching the in-memory ``Population`` representation. The transform is
+  ``w := w * p_model``; within-model normalization is recovered on read.
 """
 
 import sqlalchemy as sa
@@ -17,11 +25,24 @@
 
 
 def upgrade():
-    op.add_column(
-        table_name='populations',
-        column=sa.Column('wall_time', sa.FLOAT, nullable=True),
+    # the column may exist already, as the downgrade does not remove it
+    inspector = sa.inspect(op.get_bind())
+    columns = [col['name'] for col in inspector.get_columns('populations')]
+    if 'wall_time' not in columns:
+        op.add_column(
+            table_name='populations',
+            column=sa.Column('wall_time', sa.FLOAT, nullable=True),
+        )
+    # per-model-normalized weights -> global weights (w := w * p_model)
+    op.execute(
+        'UPDATE particles SET w = w * ('
+        'SELECT p_model FROM models WHERE models.id = particles.model_id)'
     )
 
 
 def downgrade():
-    pass
+    # global weights -> per-model-normalized weights (w := w / p_model)
+    op.execute(
+        'UPDATE particles SET w = w / ('
+        'SELECT p_model FROM models WHERE models.id = particles.model_id)'
+    )
diff --git a/pyabc/storage/numpy_bytes_storage.py b/pyabc/storage/numpy_bytes_storage.py
index c1553135..7d7c7257 100644
--- a/pyabc/storage/numpy_bytes_storage.py
+++ b/pyabc/storage/numpy_bytes_storage.py
@@ -52,6 +52,12 @@ def np_from_bytes(arr_bytes):
 
         # try to convert to primitive types
         for type_ in _primitive_types:
+            if (
+                type_ is int
+                and not np.issubdtype(arr.dtype, np.integer)
+                and arr.dtype != np.bool_
+            ):
+                continue
             try:
                 if type_(arr) == arr:
                     return type_(arr)
diff --git a/pyabc/sumstat/subset.py b/pyabc/sumstat/subset.py
index c47da004..24c668a5 100644
--- a/pyabc/sumstat/subset.py
+++ b/pyabc/sumstat/subset.py
@@ -135,7 +135,9 @@ def select(
         """Select based on GMM clusters."""
         # normalize
         if self.normalize_labels:
-            y_norm = (y - np.mean(y, axis=0)) / np.std(y, axis=0)
+            std = np.std(y, axis=0)
+            # avoid division by zero for constant (zero-variance) columns
+            y_norm = (y - np.mean(y, axis=0)) / np.where(std == 0, 1.0, std)
         else:
             y_norm = y
 
@@ -242,13 +244,18 @@ def get_augmented_subset(
     # sort remaining values by distance to reference point
     y_left = y[~in_cluster]
     distances: np.ndarray = np.linalg.norm(y_left - ref, ord=2, axis=1)
-    # indices of the required closest parameters
-    ixs_nearest: np.ndarray = np.argpartition(distances, required)[:required]
+    # indices of the required closest parameters (np.argpartition requires
+    # kth < len, so select all remaining when we need at least all of them)
+    n_left = len(distances)
+    if required >= n_left:
+        ixs_nearest = np.arange(n_left)
+    else:
+        ixs_nearest = np.argpartition(distances, required)[:required]
 
     ixs_not_in_cluster: np.ndarray = np.flatnonzero(~in_cluster)
     in_cluster[ixs_not_in_cluster[ixs_nearest]] = True
 
-    if sum(in_cluster) != desired:
+    if sum(in_cluster) != min(desired, len(y)):
         raise AssertionError('Unexpected number of entries.')
 
     return in_cluster
diff --git a/pyabc/transition/jump.py b/pyabc/transition/jump.py
index ab9f0d2c..19d679ed 100644
--- a/pyabc/transition/jump.py
+++ b/pyabc/transition/jump.py
@@ -30,7 +30,12 @@ def __init__(self, domain: np.ndarray, p_stay: float = 0.7):
         if not 0 <= p_stay <= 1:
             raise ValueError('p_stay must be in [0, 1].')
         self.p_stay = p_stay
-        self.p_move = (1 - p_stay) / (len(self.domain) - 1)
+        # guard against a single-value domain (no other value to move to)
+        self.p_move = (
+            0.0
+            if len(self.domain) == 1
+            else (1 - p_stay) / (len(self.domain) - 1)
+        )
 
         # cache a random variable (later the start index and 0 must be swapped)
         indices = np.arange(len(domain))
diff --git a/pyabc/transition/transitionmeta.py b/pyabc/transition/transitionmeta.py
index db4fa76a..7dfe5b52 100644
--- a/pyabc/transition/transitionmeta.py
+++ b/pyabc/transition/transitionmeta.py
@@ -9,13 +9,15 @@ def wrap_fit(f):
     @functools.wraps(f)
     def fit(self, X: pd.DataFrame, w: np.ndarray):
         self.X = X
-        self.w = w
         if len(X.columns) == 0:
+            self.w = w
             self.no_parameters = True
             return
         self.no_parameters = False
         if w.size > 0 and not np.isclose(w.sum(), 1):
-            w /= w.sum()
+            # normalize out-of-place so the caller's array is not mutated
+            w = w / w.sum()
+        self.w = w
         f(self, X, w)
 
     return fit
diff --git a/pyabc/util/dict2arr.py b/pyabc/util/dict2arr.py
index 4d0602f7..1bc5205e 100644
--- a/pyabc/util/dict2arr.py
+++ b/pyabc/util/dict2arr.py
@@ -41,8 +41,7 @@ def dict2arr(dct: dict | np.ndarray, keys: list) -> np.ndarray:
     if len(arr) == 1:
         return np.asarray(arr[0])
     # flatten
-    arr = [val for sub_arr in arr for val in sub_arr]
-    return np.asarray(arr)
+    return np.concatenate([np.asarray(sub_arr) for sub_arr in arr])
 
 
 def dict2arrlabels(dct: dict, keys: list) -> list[str]:
diff --git a/pyabc/visualization/credible.py b/pyabc/visualization/credible.py
index 776aa345..b570685c 100644
--- a/pyabc/visualization/credible.py
+++ b/pyabc/visualization/credible.py
@@ -14,14 +14,14 @@
 def _prepare_credible_intervals(
     history: History,
     m: int,
-    ts: list[int] | int,
-    par_names: list,
-    levels: list,
+    ts: list[int] | int | None,
+    par_names: list | None,
+    levels: list | None,
     show_mean: bool,
     show_kde_max: bool,
     show_kde_max_1d: bool,
-    kde: Transition,
-    kde_1d: Transition,
+    kde: Transition | None,
+    kde_1d: Transition | None,
 ):
     if levels is None:
         levels = [0.95]
@@ -337,8 +337,8 @@ def plot_credible_intervals_plotly(
                     error_y={
                         'type': 'data',
                         'symmetric': False,
-                        'array': cis[i_par, :, i_c] - median[i_par],
-                        'arrayminus': median[i_par] - cis[i_par, :, -1 - i_c],
+                        'array': cis[i_par, :, -1 - i_c] - median[i_par],
+                        'arrayminus': median[i_par] - cis[i_par, :, i_c],
                     },
                     mode='lines+markers',
                     marker={'color': colors[i_c]},
@@ -427,7 +427,8 @@ def plot_credible_intervals_for_time(
     if ms is None:
         ms = [0] * n_run
     elif not isinstance(ms, list) or len(ms) == 1:
-        ms = [ms] * n_run
+        # broadcast a single model id (int, or length-1 list) across runs
+        ms = [ms[0] if isinstance(ms, list) else ms] * n_run
     if levels is None:
         levels = [0.95]
     levels = sorted(levels)
@@ -512,7 +513,7 @@ def plot_credible_intervals_for_time(
                     color=f'C{i_c}',
                 )
             # reference value
-            if refvals[i_run] is not None:
+            if refvals is not None and refvals[i_run] is not None:
                 ax.plot([i_run], [refvals[i_run][par]], 'x', color='black')
         ax.set_title(f'Parameter {par}')
         # mean
diff --git a/pyabc/visualization/kde.py b/pyabc/visualization/kde.py
index 894f38b5..c0ea0b0b 100644
--- a/pyabc/visualization/kde.py
+++ b/pyabc/visualization/kde.py
@@ -246,9 +246,11 @@ def plot_kde_1d(
         xname = x
     if ax is None:
         _, ax = plt.subplots()
-    ax.plot(x_vals, pdf, **kwargs)
-    # TODO This fixes the upper bound inadequately
-    # ax.set_ylim(bottom=min(ax.get_ylim()[0], 0))
+    (line,) = ax.plot(x_vals, pdf, **kwargs)
+    # a density is non-negative, but `set_ylim` would switch off autoscaling
+    line.sticky_edges.y.append(0.0)
+    ax.update_datalim([(x_vals[0], 0.0)])
+    ax.autoscale_view()
     ax.set_xlabel(xname)
     ax.set_ylabel('Posterior')
     ax.set_xlim(xmin, xmax)
@@ -307,14 +309,11 @@ def plot_kde_1d_plotly(
         row=row,
         col=col,
     )
-    # set trace color to blue
-    # fig.update_traces(marker_color="blue", row=row, col=col)
-    # fig.add_trace(
-    #     go.Scatter(x=x_vals, y=pdf, name=xname, **kwargs),
-    #     row=row,
-    #     col=col,
-    # )
     fig.update_xaxes(title_text=xname, range=[xmin, xmax], row=row, col=col)
+    # a density is non-negative, plotly's default rangemode ignores zero
+    fig.update_yaxes(
+        title_text='Posterior', rangemode='tozero', row=row, col=col
+    )
 
     # add vertical line for reference value
     if refval is not None:
diff --git a/pyabc/weighted_statistics/weighted_statistics.py b/pyabc/weighted_statistics/weighted_statistics.py
index 4a6e1dc5..42b64d10 100644
--- a/pyabc/weighted_statistics/weighted_statistics.py
+++ b/pyabc/weighted_statistics/weighted_statistics.py
@@ -127,8 +127,8 @@ def resample(points, weights, n):
         A total of `n` points sampled from `points` with putting back
         according to `weights`.
     """
-    weights = np.asarray(weights)
-    weights /= np.sum(weights)
+    weights = np.asarray(weights, dtype=float)
+    weights = weights / np.sum(weights)
     indices = np.random.choice(
         points.shape[0], size=n, p=weights
     )  # sample index from multi-dimensional sample
diff --git a/test/base/test_epsilon.py b/test/base/test_epsilon.py
index 4b4faacf..80c8adf6 100644
--- a/test/base/test_epsilon.py
+++ b/test/base/test_epsilon.py
@@ -253,3 +253,23 @@ def model(p):
         < 3
         < pyabc.weighted_quantile(df.theta.to_numpy(), w, alpha=0.75)
     )
+
+
+def test_temperature_single_scheme_normalized():
+    """Regression: a single callable ``schemes`` argument is normalized to a
+    list, so all consumers (which iterate over ``self.schemes``) work. A bare
+    callable previously broke iteration (e.g. in ``is_adaptive``)."""
+    from pyabc.epsilon.temperature import PolynomialDecayFixedIterScheme
+
+    scheme = PolynomialDecayFixedIterScheme()
+    temp = pyabc.Temperature(schemes=scheme)
+    assert temp.schemes == [scheme]
+    # iterating over the schemes must not raise on a scalar callable
+    temp.is_adaptive()
+
+    # a list of schemes is preserved unchanged
+    schemes = [
+        PolynomialDecayFixedIterScheme(),
+        PolynomialDecayFixedIterScheme(),
+    ]
+    assert pyabc.Temperature(schemes=schemes).schemes == schemes
diff --git a/test/base/test_predictor.py b/test/base/test_predictor.py
index b5c13531..24423428 100644
--- a/test/base/test_predictor.py
+++ b/test/base/test_predictor.py
@@ -174,3 +174,24 @@ def test_wrong_input():
     """Test all kinds of wrong inputs."""
     with pytest.raises(ValueError):
         HiddenLayerHandle(method='potato')(n_in=10, n_out=10, n_sample=100)
+
+
+def test_model_selection_custom_f_score():
+    """Regression: a custom ``f_score`` passed to ``ModelSelectionPredictor``
+    must be stored and used. Previously only the default branch set
+    ``self.f_score``, so a custom scorer was silently dropped and ``fit()``
+    raised ``AttributeError``."""
+
+    def my_score(y1, y2, sigma):
+        return float(np.mean(np.abs(y1 - y2)))
+
+    msp = ModelSelectionPredictor(
+        predictors=[LinearPredictor()], f_score=my_score
+    )
+    assert msp.f_score is my_score
+
+    rng = np.random.RandomState(0)
+    x = rng.normal(size=(40, 2))
+    y = x @ np.array([[1.0], [2.0]]) + 0.01 * rng.normal(size=(40, 1))
+    msp.fit(x, y)  # would previously raise AttributeError
+    assert msp.chosen_one is not None
diff --git a/test/base/test_storage.py b/test/base/test_storage.py
index f5fc25b4..04e98bf5 100644
--- a/test/base/test_storage.py
+++ b/test/base/test_storage.py
@@ -282,6 +282,86 @@ def test_sum_stats_save_load(history: History):
     assert (sum_stats[1]['ss33'] == example_df()).all().all()
 
 
+def test_global_particle_weight_convention(history: History):
+    """Regression: particle weights are stored with the global convention
+    (summing to 1 across all particles of all models, matching the in-memory
+    ``Population``), while within-model read paths renormalize via the model
+    probability. Previously the DB stored per-model-normalized weights.
+    """
+    # two models: p_model(0)=0.4, p_model(1)=0.6; global weights sum to 1
+    particle_list = [
+        Particle(
+            m=0,
+            parameter=Parameter({'a': 1.0}),
+            weight=0.3,
+            sum_stat={'ss': 1.0},
+            distance=0.1,
+        ),
+        Particle(
+            m=0,
+            parameter=Parameter({'a': 2.0}),
+            weight=0.1,
+            sum_stat={'ss': 2.0},
+            distance=0.2,
+        ),
+        Particle(
+            m=1,
+            parameter=Parameter({'a': 3.0}),
+            weight=0.2,
+            sum_stat={'ss': 3.0},
+            distance=0.3,
+        ),
+        Particle(
+            m=1,
+            parameter=Parameter({'a': 4.0}),
+            weight=0.4,
+            sum_stat={'ss': 4.0},
+            distance=0.4,
+        ),
+    ]
+    history.append_population(
+        0, 0.5, Population(particle_list), 4, ['m0', 'm1']
+    )
+
+    # model probabilities are the per-model sums of the global weights
+    mp = history.get_model_probabilities(t=0)
+    assert np.isclose(mp.loc[0, 'p'], 0.4)
+    assert np.isclose(mp.loc[1, 'p'], 0.6)
+
+    # the RAW stored weights are global and sum to 1 (per-model storage would
+    # have summed to the number of models, i.e. 2)
+    ext = history.get_population_extended(t=0, tidy=False)
+    raw_w = ext[['particle_id', 'w']].drop_duplicates().w.values
+    assert np.isclose(raw_w.sum(), 1.0)
+    assert np.allclose(sorted(raw_w), [0.1, 0.2, 0.3, 0.4])
+
+    # within-model posterior weights sum to 1 (= g_i / p_model)
+    _, w0 = history.get_distribution(m=0, t=0)
+    _, w1 = history.get_distribution(m=1, t=0)
+    assert np.isclose(w0.sum(), 1.0) and np.isclose(w1.sum(), 1.0)
+    assert np.allclose(sorted(w0), sorted([0.3 / 0.4, 0.1 / 0.4]))
+    assert np.allclose(sorted(w1), sorted([0.2 / 0.6, 0.4 / 0.6]))
+
+    # global weighted distances / summary statistics sum to 1
+    wd = history.get_weighted_distances(t=0)
+    assert np.isclose(wd.w.sum(), 1.0)
+    assert np.allclose(sorted(wd.w.values), [0.1, 0.2, 0.3, 0.4])
+    w_ss, _ = history.get_weighted_sum_stats(t=0)
+    assert np.isclose(sum(w_ss), 1.0)
+
+    # within-model summary statistics for a single model sum to 1
+    w_m0, _ = history.get_weighted_sum_stats_for_model(m=0, t=0)
+    assert np.isclose(w_m0.sum(), 1.0)
+    assert np.allclose(sorted(w_m0), sorted([0.3 / 0.4, 0.1 / 0.4]))
+
+    # get_population reconstructs the global weights, summing to 1
+    pop = history.get_population(t=0)
+    assert np.isclose(sum(p.weight for p in pop.particles), 1.0)
+    assert np.allclose(
+        sorted(p.weight for p in pop.particles), [0.1, 0.2, 0.3, 0.4]
+    )
+
+
 def test_total_nr_samples(history: History):
     particle_list = [
         Particle(
@@ -595,3 +675,35 @@ def model(p):
     finally:
         if os.path.exists(db_file):
             os.remove(db_file)
+
+
+def test_save_dict_to_json_numpy_arrays():
+    """Regression: `save_dict_to_json` handles integer and multi-dimensional
+    numpy arrays."""
+    from pyabc.storage.json import load_dict_from_json, save_dict_to_json
+
+    f = tempfile.mkstemp(suffix='.json')[1]
+    try:
+        save_dict_to_json(
+            {1: np.array([1, 2, 3]), 2: np.array([[1.0, 2.0], [3.0, 4.0]])}, f
+        )
+        got = load_dict_from_json(f)
+    finally:
+        if os.path.exists(f):
+            os.remove(f)
+    assert got[1] == [1, 2, 3]
+    assert got[2] == [[1.0, 2.0], [3.0, 4.0]]
+
+
+def test_abcsmc_repr_with_unfinished_run():
+    """Regression: `ABCSMC.__repr__` must not crash for a run that has been
+    started but not finished (`end_time is None`)."""
+    from pyabc.storage.db_model import ABCSMC, datetime2str
+
+    row = ABCSMC(
+        id=1,
+        start_time=datetime.datetime(2024, 1, 1, 12, 0, 0),
+        end_time=None,
+    )
+    assert 'end_time=None' in repr(row)
+    assert datetime2str(None) == 'None'
diff --git a/test/base/test_sumstat.py b/test/base/test_sumstat.py
index 9acd6f9a..048d25ae 100644
--- a/test/base/test_sumstat.py
+++ b/test/base/test_sumstat.py
@@ -297,3 +297,53 @@ def model(p):
     df_info, w_info = h.get_distribution()
     off_info = abs(pyabc.weighted_mean(df_info.p0, w_info) - 0.1)
     assert off_comp > off_info
+
+
+def test_dict2arr_multi_key_concatenation():
+    """Regression/efficiency: `dict2arr` concatenates the values of several
+    keys into a single flat 1d array."""
+    dct = {
+        'a': np.array([1.0, 2.0]),
+        'b': 3.0,
+        'c': np.array([4.0, 5.0, 6.0]),
+    }
+    out = dict2arr(dct, keys=['a', 'b', 'c'])
+    assert np.array_equal(out, np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0]))
+
+
+def test_gmm_subsetter_zero_variance_labels():
+    """Regression: GMMSubsetter z-score normalization must guard against a
+    constant (zero-variance) label column instead of dividing by zero (which
+    produced NaNs that broke the subsequent GMM fit)."""
+    rng = np.random.RandomState(0)
+    n = 200
+    x = rng.normal(size=(n, 2))
+    # second label column is constant -> zero standard deviation
+    y = np.column_stack([rng.normal(size=n), np.full(n, 7.0)])
+    w = np.full(
+        (n, 1), 1.0 / n
+    )  # weights shaped (n_sample, 1), as read_sample
+    x_new, y_new, w_new = GMMSubsetter().select(x, y, w)
+    assert np.all(np.isfinite(x_new))
+    assert len(w_new) == len(x_new)
+
+
+def test_get_augmented_subset_full_fraction():
+    """Regression: `get_augmented_subset` must not pass an out-of-bounds `kth`
+    to `np.argpartition` when the required count equals the number of
+    remaining samples (`min_fraction=1.0`)."""
+    from pyabc.sumstat.subset import get_augmented_subset
+
+    y = np.arange(10, dtype=float).reshape(-1, 1)
+    ref = np.array([0.0])
+
+    in_cluster = np.zeros(10, dtype=bool)
+    in_cluster[:3] = True
+    res = get_augmented_subset(y, ref, in_cluster, min_fraction=1.0)
+    assert res.sum() == 10
+
+    # partial augmentation still selects exactly the desired count
+    ic2 = np.zeros(10, dtype=bool)
+    ic2[:2] = True
+    res2 = get_augmented_subset(y, ref, ic2, min_fraction=0.5)
+    assert res2.sum() == 5
diff --git a/test/migrate/test_migrate.py b/test/migrate/test_migrate.py
index 752a405c..94903cda 100644
--- a/test/migrate/test_migrate.py
+++ b/test/migrate/test_migrate.py
@@ -1,11 +1,37 @@
 """Migration tests."""
 
 import os
+import sqlite3
 import tempfile
 
+import numpy as np
+import pandas as pd
 import pytest
 
 import pyabc
+from pyabc.parameters import Parameter
+from pyabc.population import Particle, Population
+from pyabc.storage.version import __db_version__
+
+# model names of the test database created below
+MODEL_NAMES = ['m0', 'm1']
+
+# global particle weights per population index and model, i.e. as stored by
+#  the current format: they sum to 1 across all particles of all models,
+#  and within a model to that model's probability
+WEIGHTS = {
+    0: {0: [0.3, 0.1], 1: [0.2, 0.4]},
+    1: {0: [0.1, 0.15], 1: [0.5, 0.25]},
+}
+
+# wall times per population index
+WALL_TIMES = {0: 12.5, 1: 7.25}
+
+# SQLite can only drop table columns from version 3.35 on
+requires_drop_column = pytest.mark.skipif(
+    sqlite3.sqlite_version_info < (3, 35),
+    reason='Dropping a table column requires SQLite>=3.35',
+)
 
 
 def test_db_import(script_runner):
@@ -29,3 +55,253 @@ def test_db_import(script_runner):
 
     # remove file
     os.remove(db_file)
+
+
+def create_current_db(db_file: str) -> None:
+    """Create a database in the current format, holding two models."""
+    h = pyabc.History('sqlite:///' + db_file)
+    h.store_initial_data(None, {}, {'ss': 0.0}, {}, MODEL_NAMES, '', '', '{}')
+    for t, weights in WEIGHTS.items():
+        particles = [
+            Particle(
+                m=m,
+                parameter=Parameter({'a': float(m), 'b': float(ix)}),
+                weight=w,
+                sum_stat={'ss': float(ix)},
+                distance=0.1 * (ix + 1),
+            )
+            for m, ws in weights.items()
+            for ix, w in enumerate(ws)
+        ]
+        h.append_population(
+            t,
+            0.5 / (t + 1),
+            Population(particles),
+            10,
+            MODEL_NAMES,
+            wall_time=WALL_TIMES[t],
+        )
+
+
+def query(db_file: str, sql: str) -> list:
+    """Run a raw SQL query on the database file."""
+    con = sqlite3.connect(db_file)
+    try:
+        return con.execute(sql).fetchall()
+    finally:
+        con.close()
+
+
+def db_version(db_file: str) -> str:
+    """Storage format version of the database."""
+    return str(query(db_file, 'SELECT version_num FROM version')[0][0])
+
+
+def columns(db_file: str, table: str) -> list[str]:
+    """Column names of a database table."""
+    return [row[1] for row in query(db_file, f'PRAGMA table_info({table})')]
+
+
+def stored_weights(db_file: str) -> dict:
+    """Stored particle weights as ``{t: {m: sorted weights}}``.
+
+    The pre-population is not included.
+    """
+    rows = query(
+        db_file,
+        'SELECT populations.t, models.m, particles.w FROM particles '
+        'JOIN models ON models.id = particles.model_id '
+        'JOIN populations ON populations.id = models.population_id '
+        'WHERE populations.t >= 0',
+    )
+    weights = {}
+    for t, m, w in rows:
+        weights.setdefault(t, {}).setdefault(m, []).append(w)
+    return {
+        t: {m: sorted(ws) for m, ws in per_model.items()}
+        for t, per_model in weights.items()
+    }
+
+
+def global_weights() -> dict:
+    """Expected weights in the current format, summing to 1 over all models."""
+    return {
+        t: {m: sorted(ws) for m, ws in per_model.items()}
+        for t, per_model in WEIGHTS.items()
+    }
+
+
+def per_model_weights() -> dict:
+    """Expected weights in version 1, summing to 1 within each model."""
+    return {
+        t: {m: sorted(np.asarray(ws) / sum(ws)) for m, ws in per_model.items()}
+        for t, per_model in WEIGHTS.items()
+    }
+
+
+def assert_weights_close(actual: dict, expected: dict) -> None:
+    """Assert that two ``{t: {m: weights}}`` dictionaries match."""
+    assert actual.keys() == expected.keys()
+    for t, per_model in expected.items():
+        assert actual[t].keys() == per_model.keys()
+        for m, ws in per_model.items():
+            assert np.allclose(actual[t][m], ws)
+
+
+def alembic_config(db: str):
+    """Alembic configuration, skipping the test if alembic is missing."""
+    pytest.importorskip('alembic')
+    from pyabc.storage.migrate import _alembic_config
+
+    return _alembic_config(db)
+
+
+def to_v1(db_file: str) -> None:
+    """Turn a current-format database into a version 1 database.
+
+    Applies the version 2 downgrade, which reverts the weight normalization
+    (from global back to within-model, ``w = g_i / p_model``), and in addition
+    drops the wall time column, which did not exist in version 1 but is kept
+    by the downgrade.
+    """
+    command = pytest.importorskip('alembic.command')
+    command.downgrade(alembic_config(db_file), '1')
+
+    con = sqlite3.connect(db_file)
+    with con:
+        con.execute('ALTER TABLE populations DROP COLUMN wall_time')
+    con.close()
+
+
+@requires_drop_column
+def test_migrate_v1_to_v2(script_runner, tmp_path):
+    """Migrating a version 1 database to the current format.
+
+    Checks that the wall time column is added and that particle weights are
+    converted from the per-model to the global normalization, on a database
+    with two models, i.e. with model probabilities != 1.
+    """
+    src = str(tmp_path / 'v1.db')
+    dst = str(tmp_path / 'v2.db')
+
+    # create a database in the current format and record reference values
+    create_current_db(src)
+    h = pyabc.History('sqlite:///' + src)
+    p_models = {t: h.get_model_probabilities(t=t) for t in WEIGHTS}
+    distributions = {
+        (t, m): h.get_distribution(m=m, t=t) for t in WEIGHTS for m in [0, 1]
+    }
+    assert_weights_close(stored_weights(src), global_weights())
+
+    # turn it into a version 1 database
+    to_v1(src)
+    assert db_version(src) == '1'
+    assert 'wall_time' not in columns(src, 'populations')
+    assert_weights_close(stored_weights(src), per_model_weights())
+
+    # an outdated database cannot be imported
+    with pytest.raises(AssertionError, match='Database has version 1'):
+        pyabc.History('sqlite:///' + src)
+
+    # call the migration script
+    ret = script_runner.run(['abc-migrate', '--src', src, '--dst', dst])
+    assert ret.success
+
+    # the source database is left untouched
+    assert db_version(src) == '1'
+    assert 'wall_time' not in columns(src, 'populations')
+    assert_weights_close(stored_weights(src), per_model_weights())
+
+    # the destination database is up-to-date and has the new column
+    assert db_version(dst) == __db_version__ == '2'
+    assert 'wall_time' in columns(dst, 'populations')
+
+    # weights are back to the global normalization
+    assert_weights_close(stored_weights(dst), global_weights())
+
+    # the pre-population's dummy particle is not affected
+    assert query(
+        dst,
+        'SELECT particles.w FROM particles '
+        'JOIN models ON models.id = particles.model_id '
+        'JOIN populations ON populations.id = models.population_id '
+        'WHERE populations.t = -1',
+    ) == [(1.0,)]
+
+    # the migrated database gives the same results as the original one
+    h = pyabc.History('sqlite:///' + dst)
+    for t in WEIGHTS:
+        assert np.allclose(
+            h.get_model_probabilities(t=t).p.values, p_models[t].p.values
+        )
+        for m in [0, 1]:
+            df, w = h.get_distribution(m=m, t=t)
+            df_expected, w_expected = distributions[(t, m)]
+            pd.testing.assert_frame_equal(df, df_expected)
+            # within-model weights sum to 1
+            assert np.allclose(w, w_expected)
+            assert np.isclose(w.sum(), 1.0)
+
+    # wall times are unknown for migrated populations
+    populations = h.get_all_populations()
+    assert 'wall_time' in populations.columns
+    assert populations.wall_time.isna().all()
+
+
+def test_migrate_v2_downgrade(tmp_path):
+    """The version 2 revision can be reverted and then applied again.
+
+    The downgrade inverts the weight conversion, but keeps the wall time
+    column, so that the upgrade must tolerate an existing column.
+    """
+    command = pytest.importorskip('alembic.command')
+
+    db_file = str(tmp_path / 'db.db')
+    create_current_db(db_file)
+    cfg = alembic_config(db_file)
+
+    # revert to version 1
+    command.downgrade(cfg, '1')
+    assert db_version(db_file) == '1'
+    assert 'wall_time' in columns(db_file, 'populations')
+    assert_weights_close(stored_weights(db_file), per_model_weights())
+
+    # and migrate back to the current version
+    command.upgrade(cfg, 'head')
+    assert db_version(db_file) == __db_version__
+    assert 'wall_time' in columns(db_file, 'populations')
+    assert_weights_close(stored_weights(db_file), global_weights())
+
+
+def test_db_identifier(tmp_path):
+    """Databases can be specified as file names or as sqlite URLs."""
+    pytest.importorskip('alembic')
+    from pyabc.storage.migrate import _alembic_config, _to_db_file
+
+    db_file = str(tmp_path / 'db.db')
+
+    # file names and URLs are equivalent
+    assert _to_db_file(db_file) == _to_db_file('sqlite:///' + db_file)
+    assert _alembic_config(db_file).get_main_option(
+        'sqlalchemy.url'
+    ) == _alembic_config('sqlite:///' + db_file).get_main_option(
+        'sqlalchemy.url'
+    )
+
+    # other dialects are not supported
+    with pytest.raises(ValueError, match='only supports sqlite'):
+        _to_db_file('postgresql://user@localhost/db')
+
+
+def test_migrate_unsupported_dialect(script_runner, tmp_path):
+    """Migrating a non-sqlite database gives an error message."""
+    ret = script_runner.run(
+        [
+            'abc-migrate',
+            '--src',
+            'postgresql://user@localhost/db',
+            '--dst',
+            str(tmp_path / 'db.db'),
+        ]
+    )
+    assert 'only supports sqlite' in ret.stdout