Enhance service type validation and demand calculation

blocksnet/config/service_types/schemas.py

@@ -20,11 +20,34 @@ class UnitsSchema(DfSchema):
     site_area: Series[float] = Field(ge=0)
     build_floor_area: Series[float] = Field(ge=0)
 
-    # @classmethod
-    # def _before_validate(cls, df: pd.DataFrame):
-    #     if "parking_area" in df:
-    #         df["site_area"] += df["parking_area"]
-    #     return df
+    @classmethod
+    def _before_validate(cls, df: pd.DataFrame):
+        # Reject units that occupy neither site area nor build floor area:
+        # the optimizer treats them as a degenerate case (sort key divides by
+        # ``build_floor_area`` for ``site_area == 0`` units). Validating at
+        # load time keeps the run-time invariant intact.
+        bad = df[(df.get("site_area", 0) == 0) & (df.get("build_floor_area", 0) == 0)]
+        if not bad.empty:
+            offenders = bad.assign(_label=bad.index.astype(str))._label.tolist()
+            raise ValueError(
+                "Service units with both site_area=0 and build_floor_area=0 are not allowed: "
+                f"{offenders}"
+            )
+        # Embedded units (site_area == 0) cannot absorb a parking footprint in
+        # the standard ``site_area`` channel; fold ``parking_area`` into the
+        # build floor area as a deterministic, capacity-neutral surcharge so
+        # the parking demand participates in BFA constraints instead of being
+        # silently dropped. For units with a real site footprint, parking
+        # extends the on-site footprint.
+        if "parking_area" in df.columns:
+            embedded_mask = df["site_area"] == 0
+            df.loc[embedded_mask, "build_floor_area"] = (
+                df.loc[embedded_mask, "build_floor_area"] + df.loc[embedded_mask, "parking_area"]
+            )
+            df.loc[~embedded_mask, "site_area"] = (
+                df.loc[~embedded_mask, "site_area"] + df.loc[~embedded_mask, "parking_area"]
+            )
+        return df
 
 
 class LandUseSchema(DfSchema):

blocksnet/optimization/services/acl/facade.py

@@ -79,7 +79,6 @@ def __init__(
         self._converter.update_cap(self._capacity_checker._demands)
         self.num_params = 0
         self._last_provisions: Dict[str, float] = {}
-        self._last_blocks_services_demand: Dict[int, Dict[str, int]] | None = None
 
     @property
     def start_provisions(self) -> Dict[str, float]:
@@ -185,6 +184,10 @@ def check_constraints(self, x: ArrayLike) -> bool:
         """
         Check if a solution satisfies all constraints.
 
+        The demand added by the current candidate (new residents created by
+        the residual build floor area) is computed from ``x`` itself, so
+        that feasibility does not depend on the previous trial's state.
+
         Parameters
         ----------
         x : ArrayLike
@@ -196,8 +199,9 @@ def check_constraints(self, x: ArrayLike) -> bool:
             True if all constraints are satisfied, False otherwise.
         """
         X = self._converter(x)
+        blocks_services_demand = self.get_delta_demand(x)
         return self._area_checker.check_constraints(X) and self._capacity_checker.check_constraints(
-            X, self._last_blocks_services_demand
+            X, blocks_services_demand
         )
 
     def get_X(self, x: ArrayLike) -> List[Variable]:
@@ -249,22 +253,25 @@ def get_max_capacity(self, block_id: int, service: str) -> float:
         """
         Get the maximum capacity for a service type in a block.
 
+        Only the static catchment demand from :class:`CapacityChecker` is used
+        here. The new residents' demand is evaluated by ``check_constraints``
+        against the current ``x``; mixing it into upper bounds would cause
+        double-counting across overlapping catchments and would couple bounds
+        to the previous trial's state.
+
         Parameters
         ----------
         block_id : int
             ID of the block to query.
         service : str
-            Service type identifier.
+            The service type identifier.
 
         Returns
         -------
         float
             Maximum capacity for the specified service in the given block.
         """
-        demand = 0
-        if self._last_blocks_services_demand is not None:
-            demand = self._last_blocks_services_demand[block_id][service]
-        return self._capacity_checker.get_demand(block_id, service, demand)
+        return self._capacity_checker.get_demand(block_id, service, 0)
 
     def get_max_capacities(self, block_id: int) -> Dict[str, float]:
         """
@@ -421,10 +428,6 @@ def get_provisions(self, x_last: ArrayLike, x: ArrayLike) -> tuple[Dict[str, flo
                 provisions[st] = self._last_provisions[st]
 
         self._last_provisions = provisions
-        if self._last_blocks_services_demand is not None:
-            self._last_blocks_services_demand = None
-        else:
-            self._last_blocks_services_demand = self.get_delta_demand(x)
         self._converter.update_cap({block_id: self.get_max_capacities(block_id) for block_id in self._blocks_lu.keys()})
 
         return provisions, changed_services
@@ -461,6 +464,10 @@ def get_upper_bound_var(self, var_num: int) -> int:
         """
         Get the upper bound for a variable.
 
+        Bounds come from the static area and catchment-demand limits. The
+        delta-demand generated by the current solution is enforced by
+        ``check_constraints`` instead of being baked into bounds here.
+
         Parameters
         ----------
         var_num : int
@@ -473,10 +480,7 @@ def get_upper_bound_var(self, var_num: int) -> int:
         """
         var = self._converter.X[var_num]
         area_ub = self._area_checker.get_ub_var(var)
-        demand = 0
-        if self._last_blocks_services_demand is not None:
-            demand = self._last_blocks_services_demand[var.block_id][var.service_type]
-        demand_ub = self._capacity_checker.get_ub_var(var, demand)
+        demand_ub = self._capacity_checker.get_ub_var(var, 0)
         if demand_ub == -1:
             if area_ub == -1:
                 return 1e9
@@ -535,35 +539,34 @@ def get_delta_demand(self, x: ArrayLike) -> Dict[int, Dict[str, int]]:
             {"total_build_floor_area": "sum"}
         )  # BEFORE cutting variables_df
         for service_type in self._chosen_service_types:
-            variables_df = vars_df[vars_df.service_type == service_type]
-
-            # Aggregate capacity updates by block
-            delta_df = variables_df.groupby("block_id").agg({"total_capacity": "sum"})
-
-            # Get service type parameters
+            # Service type parameters
             _, demand, _ = service_types_config[service_type].values()
 
-            delta_df["max_population"] = 0
-
-            # 3' BFA Refill
+            # 3' BFA Refill — population is rounded ONCE per block; per-service
+            # demand is derived from that integer to keep population→demand
+            # arithmetic consistent across services.
             for block_id in self._blocks_lu.keys():
                 block_id = int(block_id)
-                if self._blocks_lu[block_id].name == "RESIDENTIAL":
-                    bfa_unit = (
-                        self._area_checker.build_floor_areas[block_id]
-                        - agg_total_build_area.loc[block_id, "total_build_floor_area"]
-                    )
-                    if block_id not in blocks_demand.keys():
-                        blocks_demand.update({block_id: {}})
-                    if service_type not in blocks_demand[block_id].keys():
-                        blocks_demand[block_id].update({service_type: 0})
-                    blocks_demand[block_id][service_type] = round_floor(
-                        (bfa_unit * (1 / BFA_COEF - 1) * demand) / (LIVING_DEMAND * 1000)
+                if self._blocks_lu[block_id].name != "RESIDENTIAL":
+                    continue
+                used_bfa = (
+                    float(agg_total_build_area.loc[block_id, "total_build_floor_area"])
+                    if block_id in agg_total_build_area.index
+                    else 0.0
+                )
+                bfa_unit = self._area_checker.build_floor_areas[block_id] - used_bfa
+                if bfa_unit <= 0:
+                    continue
+                population = round_floor((bfa_unit * (1 / BFA_COEF - 1)) / LIVING_DEMAND)
+                if population < 0:
+                    raise ValueError(
+                        f"Negative population for block {block_id}. BFA unit: {bfa_unit}"
                     )
-                    if blocks_demand[block_id][service_type] < 0:
-                        raise ValueError(
-                            f"Negative population for block {block_id} and service {service_type}. BFA unit: {bfa_unit}, demand: {demand}"
-                        )
+                if block_id not in blocks_demand.keys():
+                    blocks_demand.update({block_id: {}})
+                if service_type not in blocks_demand[block_id].keys():
+                    blocks_demand[block_id].update({service_type: 0})
+                blocks_demand[block_id][service_type] = round_floor(population / 1000 * demand)
         return blocks_demand
 
     def get_delta_population(self, x: ArrayLike) -> Dict[int, int]:

blocksnet/optimization/services/acl/provision_adapter.py

@@ -189,19 +189,30 @@ def calculate_provision(
             old_provision_df["demand"] = self.get_start_provision_df(service_type)[DEMAND_LEFT_COLUMN]
 
             if build_floor_areas is not None:
-                # 3' BFA Refill
-                for block_id in delta_df.index.unique():
+                # 3' BFA Refill — demand must include every residential block that
+                # got new residents, not only those hosting this service type.
+                # Otherwise blocks with no local unit contribute population but
+                # no demand, and the catchment LP under-weights their needs.
+                extra_population = {}
+                for block_id, land_use in self._blocks_lus.items():
+                    if land_use.name != "RESIDENTIAL":
+                        continue
                     block_id = int(block_id)
-                    if self._blocks_lus[block_id].name == "RESIDENTIAL":
-                        bfa_unit = (
-                            build_floor_areas[block_id] - agg_total_build_area.loc[block_id, "total_build_floor_area"]
-                        )
-                        delta_df.loc[block_id, "max_population"] = round_floor(
-                            (bfa_unit * (1 / BFA_COEF - 1) * demand) / (LIVING_DEMAND * 1000)
-                        )
-
-                # Update demand
-                old_provision_df.loc[delta_df.index, "demand"] += delta_df["max_population"]
+                    used_bfa = (
+                        agg_total_build_area.loc[block_id, "total_build_floor_area"]
+                        if block_id in agg_total_build_area.index
+                        else 0
+                    )
+                    bfa_unit = build_floor_areas[block_id] - used_bfa
+                    extra_population[block_id] = round_floor(
+                        (bfa_unit * (1 / BFA_COEF - 1) * demand) / (LIVING_DEMAND * 1000)
+                    )
+
+                # Update demand for every residential block with new residents
+                for block_id, pop in extra_population.items():
+                    old_provision_df.loc[block_id, "demand"] += pop
+                    if block_id in delta_df.index:
+                        delta_df.loc[block_id, "max_population"] = pop
 
             if old_provision_df["demand"].sum() == 0:
                 if old_provision_df["capacity"].sum() > 0:

blocksnet/optimization/services/core/optimizer.py

@@ -176,6 +176,7 @@ def __init__(
         constraints: Constraints,
         vars_order: VariablesOrder = None,
         vars_chooser: VariableChooser = None,
+        n_ei_candidates: int = 200,
     ):
         """
         Initialize the TPE optimizer with objective, constraints, and variable ordering.
@@ -190,6 +191,10 @@ def __init__(
             The order in which variables are optimized (default is NeutralOrder).
         vars_chooser : VariableChooser, optional
             Strategy for choosing variables to optimize (default is WeightChooser).
+        n_ei_candidates : int, optional
+            Number of EI candidates considered by TPE per trial. The default
+            is ``200`` (production sampler); verification runs may lower this
+            to keep wall time tractable.
         """
         super().__init__(objective, constraints)
 
@@ -198,7 +203,7 @@ def __init__(
             prior_weight=2,
             consider_endpoints=True,
             n_startup_trials=0,
-            n_ei_candidates=200,
+            n_ei_candidates=n_ei_candidates,
         )
         self._study = optuna.create_study(direction=optuna.study.StudyDirection.MAXIMIZE, sampler=sampler)