Merge pull request #219 from daveshah1/placer_heap

HeAP-based analytical placer

41 files changed, 2565 insertions, 215 deletions

diff --git a/.cirrus/Dockerfile.ubuntu16.04 b/.cirrus/Dockerfile.ubuntu16.04
index 1b93cfb8..0c8201b8 100644
--- a/.cirrus/Dockerfile.ubuntu16.04
+++ b/.cirrus/Dockerfile.ubuntu16.04
@@ -8,7 +8,7 @@ RUN set -e -x ;\
     apt-get -y install \
         build-essential autoconf cmake clang bison wget flex gperf \
         libreadline-dev gawk tcl-dev libffi-dev graphviz xdot python3-dev \
-        libboost-all-dev qt5-default git libftdi-dev pkg-config
+        libboost-all-dev qt5-default git libftdi-dev pkg-config libeigen3-dev
 
 RUN set -e -x ;\
     mkdir -p /usr/local/src ;\
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4f29d132..ade76d60 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,6 +5,8 @@ project(nextpnr)
 option(BUILD_GUI "Build GUI" ON)
 option(BUILD_PYTHON "Build Python Integration" ON)
 option(BUILD_TESTS "Build GUI" OFF)
+option(BUILD_HEAP "Build HeAP analytic placer" ON)
+option(USE_OPENMP "Use OpenMP to accelerate analytic placer" OFF)
 option(COVERAGE "Add code coverage info" OFF)
 option(STATIC_BUILD "Create static build" OFF)
 option(EXTERNAL_CHIPDB "Create build with pre-built chipdb binaries" OFF)
@@ -53,12 +55,16 @@ endforeach()
 
 set(CMAKE_CXX_STANDARD 11)
 if (MSVC)
-set(CMAKE_CONFIGURATION_TYPES "Debug;Release" CACHE STRING "" FORCE)
-set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} /D_DEBUG /W4 /wd4100 /wd4244 /wd4125 /wd4800 /wd4456 /wd4458 /wd4305 /wd4459 /wd4121 /wd4996")
-set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} /W4 /wd4100 /wd4244 /wd4125 /wd4800 /wd4456 /wd4458 /wd4305 /wd4459 /wd4121 /wd4996 /wd4127")
+    set(CMAKE_CONFIGURATION_TYPES "Debug;Release" CACHE STRING "" FORCE)
+    set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} /D_DEBUG /W4 /wd4100 /wd4244 /wd4125 /wd4800 /wd4456 /wd4458 /wd4305 /wd4459 /wd4121 /wd4996")
+    set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} /W4 /wd4100 /wd4244 /wd4125 /wd4800 /wd4456 /wd4458 /wd4305 /wd4459 /wd4121 /wd4996 /wd4127")
 else()
-set(CMAKE_CXX_FLAGS_DEBUG "-Wall -fPIC -ggdb -pipe")
-set(CMAKE_CXX_FLAGS_RELEASE "-Wall -fPIC -O3 -g -pipe")
+    set(CMAKE_CXX_FLAGS_DEBUG "-Wall -fPIC -ggdb -pipe")
+    if (USE_OPENMP)
+        set(CMAKE_CXX_FLAGS_RELEASE "-Wall -fPIC -O3 -g -pipe -fopenmp")
+    else()
+        set(CMAKE_CXX_FLAGS_RELEASE "-Wall -fPIC -O3 -g -pipe")
+    endif()
 endif()
 set(CMAKE_DEFIN)
 
@@ -181,6 +187,14 @@ if (BUILD_PYTHON)
 endif()
 
 include_directories(common/ json/ ${Boost_INCLUDE_DIRS} ${PYTHON_INCLUDE_DIRS})
+
+if(BUILD_HEAP)
+    find_package (Eigen3 REQUIRED NO_MODULE)
+    include_directories(${EIGEN3_INCLUDE_DIRS})
+    add_definitions(${EIGEN3_DEFINITIONS})
+    add_definitions(-DWITH_HEAP)
+endif()
+
 aux_source_directory(common/ COMMON_SRC_FILES)
 aux_source_directory(json/ JSON_PARSER_FILES)
 set(COMMON_FILES ${COMMON_SRC_FILES} ${JSON_PARSER_FILES})
diff --git a/README.md b/README.md
index 5b79d1fb..0417e6f7 100644
--- a/README.md
+++ b/README.md
@@ -36,6 +36,7 @@ of the selected architecture:
 - Python 3.5 or later, including development libraries (`python3-dev` for Ubuntu)
   - on Windows make sure to install same version as supported by [vcpkg](https://github.com/Microsoft/vcpkg/blob/master/ports/python3/CONTROL)
 - Boost libraries (`libboost-dev libboost-filesystem-dev libboost-thread-dev libboost-program-options-dev libboost-python-dev libboost-dev` or `libboost-all-dev` for Ubuntu)
+- Eigen3 (`libeigen3-dev` for Ubuntu) is required to build the analytic placer
 - Latest git Yosys is required to synthesise the demo design
 - For building on Windows with MSVC, usage of vcpkg is advised for dependency installation.
   - For 32 bit builds: `vcpkg install boost-filesystem boost-program-options boost-thread boost-python qt5-base`
@@ -119,11 +120,11 @@ Use cmake `-D` options to specify which version of nextpnr you want to build.
 Use `-DARCH=...` to set the architecture. It is a semicolon separated list.
 Use `cmake . -DARCH=all` to build all supported architectures.
 
-The following runs a debug build of the iCE40 architecture without GUI
-and without Python support and only HX1K support:
+The following runs a debug build of the iCE40 architecture without GUI,
+ without Python support, without the HeAP analytic placer and only HX1K support:
 
 ```
-cmake -DARCH=ice40 -DCMAKE_BUILD_TYPE=Debug -DBUILD_PYTHON=OFF -DBUILD_GUI=OFF -DICE40_HX1K_ONLY=1 .
+cmake -DARCH=ice40 -DCMAKE_BUILD_TYPE=Debug -DBUILD_PYTHON=OFF -DBUILD_GUI=OFF -DBUILD_HEAP=OFF -DICE40_HX1K_ONLY=1 .
 make -j$(nproc)
 ```
 
@@ -134,6 +135,9 @@ cmake -DARCH=ice40 -DBUILD_PYTHON=OFF -DBUILD_GUI=OFF -DSTATIC_BUILD=ON .
 make -j$(nproc)
 ```
 
+The HeAP placer's solver can optionally use OpenMP for a speedup on very large designs. Enable this by passing
+`-DUSE_OPENMP=yes` to cmake (compiler support may vary).
+
 You can change the location where nextpnr will be installed (this will usually default to `/usr/local`) by using
 `-DCMAKE_INSTALL_PREFIX=/install/prefix`.
 
diff --git a/common/command.cc b/common/command.cc
index 1399efdb..6f4137fe 100644
--- a/common/command.cc
+++ b/common/command.cc
@@ -27,6 +27,7 @@
 #include "pybindings.h"
 #endif
 
+#include <boost/algorithm/string/join.hpp>
 #include <boost/filesystem/convenience.hpp>
 #include <boost/program_options.hpp>
 #include <fstream>
@@ -120,8 +121,18 @@ po::options_description CommandHandler::getGeneralOptions()
     general.add_options()("json", po::value<std::string>(), "JSON design file to ingest");
     general.add_options()("seed", po::value<int>(), "seed value for random number generator");
     general.add_options()("randomize-seed,r", "randomize seed value for random number generator");
+
+    general.add_options()(
+            "placer", po::value<std::string>(),
+            std::string("placer algorithm to use; available: " + boost::algorithm::join(Arch::availablePlacers, ", ") +
+                        "; default: " + Arch::defaultPlacer)
+                    .c_str());
+
     general.add_options()("slack_redist_iter", po::value<int>(), "number of iterations between slack redistribution");
     general.add_options()("cstrweight", po::value<float>(), "placer weighting for relative constraint satisfaction");
+    general.add_options()("starttemp", po::value<float>(), "placer SA start temperature");
+    general.add_options()("placer-budgets", "use budget rather than criticality in placer timing weights");
+
     general.add_options()("pack-only", "pack design only without placement or routing");
 
     general.add_options()("ignore-loops", "ignore combinational loops in timing analysis");
@@ -183,10 +194,27 @@ void CommandHandler::setupContext(Context *ctx)
         settings->set("timing/allowFail", true);
     }
 
+    if (vm.count("placer")) {
+        std::string placer = vm["placer"].as<std::string>();
+        if (std::find(Arch::availablePlacers.begin(), Arch::availablePlacers.end(), placer) ==
+            Arch::availablePlacers.end())
+            log_error("Placer algorithm '%s' is not supported (available options: %s)\n", placer.c_str(),
+                      boost::algorithm::join(Arch::availablePlacers, ", ").c_str());
+        settings->set("placer", placer);
+    } else {
+        settings->set("placer", Arch::defaultPlacer);
+    }
+
     if (vm.count("cstrweight")) {
         settings->set("placer1/constraintWeight", vm["cstrweight"].as<float>());
     }
+    if (vm.count("starttemp")) {
+        settings->set("placer1/startTemp", vm["starttemp"].as<float>());
+    }
 
+    if (vm.count("placer-budgets")) {
+        settings->set("placer1/budgetBased", true);
+    }
     if (vm.count("freq")) {
         auto freq = vm["freq"].as<double>();
         if (freq > 0)
diff --git a/common/nextpnr.cc b/common/nextpnr.cc
index bb941d3d..54333b15 100644
--- a/common/nextpnr.cc
+++ b/common/nextpnr.cc
@@ -221,6 +221,9 @@ delay_t Context::getNetinfoRouteDelay(const NetInfo *net_info, const PortRef &us
         return 0;
 #endif
 
+    if (net_info->wires.empty())
+        return predictDelay(net_info, user_info);
+
     WireId src_wire = getNetinfoSourceWire(net_info);
     if (src_wire == WireId())
         return 0;
@@ -421,4 +424,25 @@ void BaseCtx::addClock(IdString net, float freq)
     }
 }
 
+void BaseCtx::createRectangularRegion(IdString name, int x0, int y0, int x1, int y1)
+{
+    std::unique_ptr<Region> new_region(new Region());
+    new_region->name = name;
+    new_region->constr_bels = true;
+    new_region->constr_pips = false;
+    new_region->constr_wires = false;
+    for (int x = x0; x <= x1; x++) {
+        for (int y = y0; y <= y1; y++) {
+            for (auto bel : getCtx()->getBelsByTile(x, y))
+                new_region->bels.insert(bel);
+        }
+    }
+    region[name] = std::move(new_region);
+}
+void BaseCtx::addBelToRegion(IdString name, BelId bel) { region[name]->bels.insert(bel); }
+void BaseCtx::constrainCellToRegion(IdString cell, IdString region_name)
+{
+    cells[cell]->region = region[region_name].get();
+}
+
 NEXTPNR_NAMESPACE_END
diff --git a/common/nextpnr.h b/common/nextpnr.h
index d58ae529..5967ecee 100644
--- a/common/nextpnr.h
+++ b/common/nextpnr.h
@@ -637,6 +637,9 @@ struct BaseCtx
 
     // Intended to simplify Python API
     void addClock(IdString net, float freq);
+    void createRectangularRegion(IdString name, int x0, int y0, int x1, int y1);
+    void addBelToRegion(IdString name, BelId bel);
+    void constrainCellToRegion(IdString cell, IdString region_name);
 };
 
 NEXTPNR_NAMESPACE_END
diff --git a/common/place_common.cc b/common/place_common.cc
index b3eb4267..73a320d0 100644
--- a/common/place_common.cc
+++ b/common/place_common.cc
@@ -304,7 +304,7 @@ class ConstraintLegaliseWorker
     // Set the strength to locked on all cells in chain
     void lockdown_chain(CellInfo *root)
     {
-        root->belStrength = STRENGTH_LOCKED;
+        root->belStrength = STRENGTH_STRONG;
         for (auto child : root->constr_children)
             lockdown_chain(child);
     }
@@ -380,7 +380,7 @@ class ConstraintLegaliseWorker
                                 rippedCells.insert(confl_cell->name);
                             }
                         }
-                        ctx->bindBel(target, ctx->cells.at(cp.first).get(), STRENGTH_LOCKED);
+                        ctx->bindBel(target, ctx->cells.at(cp.first).get(), STRENGTH_STRONG);
                         rippedCells.erase(cp.first);
                     }
                     for (auto cp : solution) {
@@ -529,4 +529,12 @@ int get_constraints_distance(const Context *ctx, const CellInfo *cell)
     return dist;
 }
 
+bool check_cell_bel_region(const CellInfo *cell, BelId bel)
+{
+    if (cell->region != nullptr && cell->region->constr_bels && !cell->region->bels.count(bel))
+        return false;
+    else
+        return true;
+}
+
 NEXTPNR_NAMESPACE_END
diff --git a/common/place_common.h b/common/place_common.h
index 79dec067..fa5ce4c2 100644
--- a/common/place_common.h
+++ b/common/place_common.h
@@ -49,6 +49,10 @@ bool legalise_relative_constraints(Context *ctx);
 
 // Get the total distance from satisfied constraints for a cell
 int get_constraints_distance(const Context *ctx, const CellInfo *cell);
+
+// Check that a Bel is within the region for a cell
+bool check_cell_bel_region(const CellInfo *cell, BelId bel);
+
 NEXTPNR_NAMESPACE_END
 
 #endif
diff --git a/common/placer1.cc b/common/placer1.cc
index 5b72602f..98251627 100644
--- a/common/placer1.cc
+++ b/common/placer1.cc
@@ -24,6 +24,7 @@
 #include "placer1.h"
 #include <algorithm>
 #include <boost/lexical_cast.hpp>
+#include <boost/range/adaptor/reversed.hpp>
 #include <chrono>
 #include <cmath>
 #include <iostream>
@@ -43,10 +44,32 @@
 #include "timing.h"
 #include "util.h"
 
+namespace std {
+template <> struct hash<std::pair<NEXTPNR_NAMESPACE_PREFIX IdString, std::size_t>>
+{
+    std::size_t operator()(const std::pair<NEXTPNR_NAMESPACE_PREFIX IdString, std::size_t> &idp) const noexcept
+    {
+        std::size_t seed = 0;
+        boost::hash_combine(seed, hash<NEXTPNR_NAMESPACE_PREFIX IdString>()(idp.first));
+        boost::hash_combine(seed, hash<std::size_t>()(idp.second));
+        return seed;
+    }
+};
+} // namespace std
+
 NEXTPNR_NAMESPACE_BEGIN
 
 class SAPlacer
 {
+  private:
+    struct BoundingBox
+    {
+        int x0 = 0, x1 = 0, y0 = 0, y1 = 0;
+        bool is_inside_inc(int x, int y) const { return x >= x0 && x <= x1 && y >= y0 && y <= y1; }
+        bool touches_bounds(int x, int y) const { return x == x0 || x == x1 || y == y0 || y == y1; }
+        wirelen_t hpwl() const { return wirelen_t((x1 - x0) + (y1 - y0)); }
+    };
+
   public:
     SAPlacer(Context *ctx, Placer1Cfg cfg) : ctx(ctx), cfg(cfg)
     {
@@ -78,13 +101,44 @@ class SAPlacer
         }
         diameter = std::max(max_x, max_y) + 1;
 
-        costs.resize(ctx->nets.size());
+        net_bounds.resize(ctx->nets.size());
+        net_arc_tcost.resize(ctx->nets.size());
+        moveChange.already_bounds_changed.resize(ctx->nets.size());
+        moveChange.already_changed_arcs.resize(ctx->nets.size());
         old_udata.reserve(ctx->nets.size());
+        net_by_udata.reserve(ctx->nets.size());
         decltype(NetInfo::udata) n = 0;
         for (auto &net : ctx->nets) {
             old_udata.emplace_back(net.second->udata);
+            net_arc_tcost.at(n).resize(net.second->users.size());
+            moveChange.already_changed_arcs.at(n).resize(net.second->users.size());
             net.second->udata = n++;
+            net_by_udata.push_back(net.second.get());
+        }
+        for (auto &region : sorted(ctx->region)) {
+            Region *r = region.second;
+            BoundingBox bb;
+            if (r->constr_bels) {
+                bb.x0 = std::numeric_limits<int>::max();
+                bb.x1 = std::numeric_limits<int>::min();
+                bb.y0 = std::numeric_limits<int>::max();
+                bb.y1 = std::numeric_limits<int>::min();
+                for (auto bel : r->bels) {
+                    Loc loc = ctx->getBelLocation(bel);
+                    bb.x0 = std::min(bb.x0, loc.x);
+                    bb.x1 = std::max(bb.x1, loc.x);
+                    bb.y0 = std::min(bb.y0, loc.y);
+                    bb.y1 = std::max(bb.y1, loc.y);
+                }
+            } else {
+                bb.x0 = 0;
+                bb.y0 = 0;
+                bb.x1 = max_x;
+                bb.y1 = max_y;
+            }
+            region_bounds[r->name] = bb;
         }
+        build_port_index();
     }
 
     ~SAPlacer()
@@ -93,97 +147,121 @@ class SAPlacer
             net.second->udata = old_udata[net.second->udata];
     }
 
-    bool place()
+    bool place(bool refine = false)
     {
         log_break();
         ctx->lock();
 
         size_t placed_cells = 0;
-        // Initial constraints placer
-        for (auto &cell_entry : ctx->cells) {
-            CellInfo *cell = cell_entry.second.get();
-            auto loc = cell->attrs.find(ctx->id("BEL"));
-            if (loc != cell->attrs.end()) {
-                std::string loc_name = loc->second;
-                BelId bel = ctx->getBelByName(ctx->id(loc_name));
-                if (bel == BelId()) {
-                    log_error("No Bel named \'%s\' located for "
-                              "this chip (processing BEL attribute on \'%s\')\n",
-                              loc_name.c_str(), cell->name.c_str(ctx));
-                }
+        std::vector<CellInfo *> autoplaced;
+        std::vector<CellInfo *> chain_basis;
+        if (!refine) {
+            // Initial constraints placer
+            for (auto &cell_entry : ctx->cells) {
+                CellInfo *cell = cell_entry.second.get();
+                auto loc = cell->attrs.find(ctx->id("BEL"));
+                if (loc != cell->attrs.end()) {
+                    std::string loc_name = loc->second;
+                    BelId bel = ctx->getBelByName(ctx->id(loc_name));
+                    if (bel == BelId()) {
+                        log_error("No Bel named \'%s\' located for "
+                                  "this chip (processing BEL attribute on \'%s\')\n",
+                                  loc_name.c_str(), cell->name.c_str(ctx));
+                    }
 
-                IdString bel_type = ctx->getBelType(bel);
-                if (bel_type != cell->type) {
-                    log_error("Bel \'%s\' of type \'%s\' does not match cell "
-                              "\'%s\' of type \'%s\'\n",
-                              loc_name.c_str(), bel_type.c_str(ctx), cell->name.c_str(ctx), cell->type.c_str(ctx));
-                }
-                if (!ctx->isValidBelForCell(cell, bel)) {
-                    log_error("Bel \'%s\' of type \'%s\' is not valid for cell "
-                              "\'%s\' of type \'%s\'\n",
-                              loc_name.c_str(), bel_type.c_str(ctx), cell->name.c_str(ctx), cell->type.c_str(ctx));
-                }
+                    IdString bel_type = ctx->getBelType(bel);
+                    if (bel_type != cell->type) {
+                        log_error("Bel \'%s\' of type \'%s\' does not match cell "
+                                  "\'%s\' of type \'%s\'\n",
+                                  loc_name.c_str(), bel_type.c_str(ctx), cell->name.c_str(ctx), cell->type.c_str(ctx));
+                    }
+                    if (!ctx->isValidBelForCell(cell, bel)) {
+                        log_error("Bel \'%s\' of type \'%s\' is not valid for cell "
+                                  "\'%s\' of type \'%s\'\n",
+                                  loc_name.c_str(), bel_type.c_str(ctx), cell->name.c_str(ctx), cell->type.c_str(ctx));
+                    }
 
-                auto bound_cell = ctx->getBoundBelCell(bel);
-                if (bound_cell) {
-                    log_error("Cell \'%s\' cannot be bound to bel \'%s\' since it is already bound to cell \'%s\'\n",
-                              cell->name.c_str(ctx), loc_name.c_str(), bound_cell->name.c_str(ctx));
-                }
+                    auto bound_cell = ctx->getBoundBelCell(bel);
+                    if (bound_cell) {
+                        log_error(
+                                "Cell \'%s\' cannot be bound to bel \'%s\' since it is already bound to cell \'%s\'\n",
+                                cell->name.c_str(ctx), loc_name.c_str(), bound_cell->name.c_str(ctx));
+                    }
 
-                ctx->bindBel(bel, cell, STRENGTH_USER);
-                locked_bels.insert(bel);
-                placed_cells++;
+                    ctx->bindBel(bel, cell, STRENGTH_USER);
+                    locked_bels.insert(bel);
+                    placed_cells++;
+                }
             }
-        }
-        int constr_placed_cells = placed_cells;
-        log_info("Placed %d cells based on constraints.\n", int(placed_cells));
-        ctx->yield();
+            int constr_placed_cells = placed_cells;
+            log_info("Placed %d cells based on constraints.\n", int(placed_cells));
+            ctx->yield();
 
-        // Sort to-place cells for deterministic initial placement
-        std::vector<CellInfo *> autoplaced;
-        for (auto &cell : ctx->cells) {
-            CellInfo *ci = cell.second.get();
-            if (ci->bel == BelId()) {
-                autoplaced.push_back(cell.second.get());
+            // Sort to-place cells for deterministic initial placement
+
+            for (auto &cell : ctx->cells) {
+                CellInfo *ci = cell.second.get();
+                if (ci->bel == BelId()) {
+                    autoplaced.push_back(cell.second.get());
+                }
             }
-        }
-        std::sort(autoplaced.begin(), autoplaced.end(), [](CellInfo *a, CellInfo *b) { return a->name < b->name; });
-        ctx->shuffle(autoplaced);
-        auto iplace_start = std::chrono::high_resolution_clock::now();
-        // Place cells randomly initially
-        log_info("Creating initial placement for remaining %d cells.\n", int(autoplaced.size()));
-
-        for (auto cell : autoplaced) {
-            place_initial(cell);
-            placed_cells++;
-            if ((placed_cells - constr_placed_cells) % 500 == 0)
+            std::sort(autoplaced.begin(), autoplaced.end(), [](CellInfo *a, CellInfo *b) { return a->name < b->name; });
+            ctx->shuffle(autoplaced);
+            auto iplace_start = std::chrono::high_resolution_clock::now();
+            // Place cells randomly initially
+            log_info("Creating initial placement for remaining %d cells.\n", int(autoplaced.size()));
+
+            for (auto cell : autoplaced) {
+                place_initial(cell);
+                placed_cells++;
+                if ((placed_cells - constr_placed_cells) % 500 == 0)
+                    log_info("  initial placement placed %d/%d cells\n", int(placed_cells - constr_placed_cells),
+                             int(autoplaced.size()));
+            }
+            if ((placed_cells - constr_placed_cells) % 500 != 0)
                 log_info("  initial placement placed %d/%d cells\n", int(placed_cells - constr_placed_cells),
                          int(autoplaced.size()));
+            if (cfg.budgetBased && ctx->slack_redist_iter > 0)
+                assign_budget(ctx);
+            ctx->yield();
+            auto iplace_end = std::chrono::high_resolution_clock::now();
+            log_info("Initial placement time %.02fs\n",
+                     std::chrono::duration<float>(iplace_end - iplace_start).count());
+            log_info("Running simulated annealing placer.\n");
+        } else {
+            for (auto &cell : ctx->cells) {
+                CellInfo *ci = cell.second.get();
+                if (ci->belStrength > STRENGTH_STRONG)
+                    continue;
+                else if (ci->constr_parent != nullptr)
+                    continue;
+                else if (!ci->constr_children.empty() || ci->constr_z != ci->UNCONSTR)
+                    chain_basis.push_back(ci);
+                else
+                    autoplaced.push_back(ci);
+            }
+            require_legal = false;
+            diameter = 3;
+            log_info("Running simulated annealing placer for refinement.\n");
         }
-        if ((placed_cells - constr_placed_cells) % 500 != 0)
-            log_info("  initial placement placed %d/%d cells\n", int(placed_cells - constr_placed_cells),
-                     int(autoplaced.size()));
-        if (ctx->slack_redist_iter > 0)
-            assign_budget(ctx);
-        ctx->yield();
-        auto iplace_end = std::chrono::high_resolution_clock::now();
-        log_info("Initial placement time %.02fs\n", std::chrono::duration<float>(iplace_end - iplace_start).count());
         auto saplace_start = std::chrono::high_resolution_clock::now();
-        log_info("Running simulated annealing placer.\n");
 
-        // Calculate metric after initial placement
-        curr_metric = 0;
-        curr_tns = 0;
-        for (auto &net : ctx->nets) {
-            wirelen_t wl = get_net_metric(ctx, net.second.get(), MetricType::COST, curr_tns);
-            costs[net.second->udata] = CostChange{wl, -1};
-            curr_metric += wl;
-        }
+        // Invoke timing analysis to obtain criticalities
+        if (!cfg.budgetBased)
+            get_criticalities(ctx, &net_crit);
+
+        // Calculate costs after initial placement
+        setup_costs();
+        curr_wirelen_cost = total_wirelen_cost();
+        curr_timing_cost = total_timing_cost();
+        last_wirelen_cost = curr_wirelen_cost;
+        last_timing_cost = curr_timing_cost;
+
+        wirelen_t avg_wirelen = curr_wirelen_cost;
+        wirelen_t min_wirelen = curr_wirelen_cost;
 
         int n_no_progress = 0;
-        wirelen_t min_metric = curr_metric;
-        double avg_metric = curr_metric;
-        temp = 10000;
+        temp = refine ? 1e-7 : cfg.startTemp;
 
         // Main simulated annealing loop
         for (int iter = 1;; iter++) {
@@ -191,9 +269,9 @@ class SAPlacer
             improved = false;
 
             if (iter % 5 == 0 || iter == 1)
-                log_info("  at iteration #%d: temp = %f, cost = "
-                         "%.0f, est tns = %.02fns\n",
-                         iter, temp, double(curr_metric), curr_tns);
+                log_info("  at iteration #%d: temp = %f, timing cost = "
+                         "%.0f, wirelen = %.0f\n",
+                         iter, temp, double(curr_timing_cost), double(curr_wirelen_cost));
 
             for (int m = 0; m < 15; ++m) {
                 // Loop through all automatically placed cells
@@ -205,10 +283,17 @@ class SAPlacer
                     if (try_bel != BelId() && try_bel != cell->bel)
                         try_swap_position(cell, try_bel);
                 }
+                // Also try swapping chains, if applicable
+                for (auto cb : chain_basis) {
+                    Loc chain_base_loc = ctx->getBelLocation(cb->bel);
+                    BelId try_base = random_bel_for_cell(cb, chain_base_loc.z);
+                    if (try_base != BelId() && try_base != cb->bel)
+                        try_swap_chain(cb, try_base);
+                }
             }
 
-            if (curr_metric < min_metric) {
-                min_metric = curr_metric;
+            if (curr_wirelen_cost < min_wirelen) {
+                min_wirelen = curr_wirelen_cost;
                 improved = true;
             }
 
@@ -218,9 +303,10 @@ class SAPlacer
             else
                 n_no_progress++;
 
-            if (temp <= 1e-3 && n_no_progress >= 5) {
-                if (iter % 5 != 0)
-                    log_info("  at iteration #%d: temp = %f, cost = %f\n", iter, temp, double(curr_metric));
+            if (temp <= 1e-7 && n_no_progress >= (refine ? 1 : 5)) {
+                log_info("  at iteration #%d: temp = %f, timing cost = "
+                         "%.0f, wirelen = %.0f \n",
+                         iter, temp, double(curr_timing_cost), double(curr_wirelen_cost));
                 break;
             }
 
@@ -228,61 +314,64 @@ class SAPlacer
 
             int M = std::max(max_x, max_y) + 1;
 
-            double upper = 0.6, lower = 0.4;
+            if (ctx->verbose)
+                log("iter #%d: temp = %f, timing cost = "
+                    "%.0f, wirelen = %.0f, dia = %d, Ra = %.02f \n",
+                    iter, temp, double(curr_timing_cost), double(curr_wirelen_cost), diameter, Raccept);
 
-            if (curr_metric < 0.95 * avg_metric && curr_metric > 0) {
-                avg_metric = 0.8 * avg_metric + 0.2 * curr_metric;
+            if (curr_wirelen_cost < 0.95 * avg_wirelen && curr_wirelen_cost > 0) {
+                avg_wirelen = 0.8 * avg_wirelen + 0.2 * curr_wirelen_cost;
             } else {
-                if (Raccept >= 0.8) {
-                    temp *= 0.7;
-                } else if (Raccept > upper) {
-                    if (diameter < M)
-                        diameter++;
-                    else
-                        temp *= 0.9;
-                } else if (Raccept > lower) {
+                double diam_next = diameter * (1.0 - 0.44 + Raccept);
+                diameter = std::max<int>(1, std::min<int>(M, int(diam_next + 0.5)));
+                if (Raccept > 0.96) {
+                    temp *= 0.5;
+                } else if (Raccept > 0.8) {
+                    temp *= 0.9;
+                } else if (Raccept > 0.15 && diameter > 1) {
                     temp *= 0.95;
                 } else {
-                    // Raccept < 0.3
-                    if (diameter > 1)
-                        diameter--;
-                    else
-                        temp *= 0.8;
+                    temp *= 0.8;
                 }
             }
             // Once cooled below legalise threshold, run legalisation and start requiring
             // legal moves only
-            if (temp < legalise_temp && require_legal) {
+            if (diameter < legalise_dia && require_legal) {
                 if (legalise_relative_constraints(ctx)) {
                     // Only increase temperature if something was moved
                     autoplaced.clear();
+                    chain_basis.clear();
                     for (auto cell : sorted(ctx->cells)) {
-                        if (cell.second->belStrength < STRENGTH_STRONG)
+                        if (cell.second->belStrength <= STRENGTH_STRONG && cell.second->constr_parent == nullptr &&
+                            !cell.second->constr_children.empty())
+                            chain_basis.push_back(cell.second);
+                        else if (cell.second->belStrength < STRENGTH_STRONG)
                             autoplaced.push_back(cell.second);
                     }
-                    temp = post_legalise_temp;
-                    diameter *= post_legalise_dia_scale;
+                    // temp = post_legalise_temp;
+                    // diameter = std::min<int>(M, diameter * post_legalise_dia_scale);
                     ctx->shuffle(autoplaced);
 
                     // Legalisation is a big change so force a slack redistribution here
-                    if (ctx->slack_redist_iter > 0)
+                    if (ctx->slack_redist_iter > 0 && cfg.budgetBased)
                         assign_budget(ctx, true /* quiet */);
                 }
                 require_legal = false;
-            } else if (ctx->slack_redist_iter > 0 && iter % ctx->slack_redist_iter == 0) {
+            } else if (cfg.budgetBased && ctx->slack_redist_iter > 0 && iter % ctx->slack_redist_iter == 0) {
                 assign_budget(ctx, true /* quiet */);
             }
 
+            // Invoke timing analysis to obtain criticalities
+            if (!cfg.budgetBased && ctx->timing_driven)
+                get_criticalities(ctx, &net_crit);
+            // Need to rebuild costs after criticalities change
+            setup_costs();
             // Recalculate total metric entirely to avoid rounding errors
             // accumulating over time
-            curr_metric = 0;
-            curr_tns = 0;
-            for (auto &net : ctx->nets) {
-                wirelen_t wl = get_net_metric(ctx, net.second.get(), MetricType::COST, curr_tns);
-                costs[net.second->udata] = CostChange{wl, -1};
-                curr_metric += wl;
-            }
-
+            curr_wirelen_cost = total_wirelen_cost();
+            curr_timing_cost = total_timing_cost();
+            last_wirelen_cost = curr_wirelen_cost;
+            last_timing_cost = curr_timing_cost;
             // Let the UI show visualization updates.
             ctx->yield();
         }
@@ -334,7 +423,8 @@ class SAPlacer
                 ctx->unbindBel(cell->bel);
             }
             IdString targetType = cell->type;
-            for (auto bel : ctx->getBels()) {
+
+            auto proc_bel = [&](BelId bel) {
                 if (ctx->getBelType(bel) == targetType && ctx->isValidBelForCell(cell, bel)) {
                     if (ctx->checkBelAvail(bel)) {
                         uint64_t score = ctx->rng64();
@@ -352,7 +442,18 @@ class SAPlacer
                         }
                     }
                 }
+            };
+
+            if (cell->region != nullptr && cell->region->constr_bels) {
+                for (auto bel : cell->region->bels) {
+                    proc_bel(bel);
+                }
+            } else {
+                for (auto bel : ctx->getBels()) {
+                    proc_bel(bel);
+                }
             }
+
             if (best_bel == BelId()) {
                 if (iters == 0 || ripup_bel == BelId())
                     log_error("failed to place cell '%s' of type '%s'\n", cell->name.c_str(ctx), cell->type.c_str(ctx));
@@ -373,49 +474,38 @@ class SAPlacer
     // Attempt a SA position swap, return true on success or false on failure
     bool try_swap_position(CellInfo *cell, BelId newBel)
     {
-        static std::vector<NetInfo *> updates;
-        updates.clear();
+        static const double epsilon = 1e-20;
+        moveChange.reset();
+        if (!require_legal && is_constrained(cell))
+            return false;
         BelId oldBel = cell->bel;
         CellInfo *other_cell = ctx->getBoundBelCell(newBel);
-        if (other_cell != nullptr && other_cell->belStrength > STRENGTH_WEAK) {
+        if (!require_legal && other_cell != nullptr &&
+            (is_constrained(other_cell) || other_cell->belStrength > STRENGTH_WEAK)) {
             return false;
         }
         int old_dist = get_constraints_distance(ctx, cell);
         int new_dist;
         if (other_cell != nullptr)
             old_dist += get_constraints_distance(ctx, other_cell);
-        wirelen_t new_metric = 0, delta;
+        double delta = 0;
         ctx->unbindBel(oldBel);
         if (other_cell != nullptr) {
             ctx->unbindBel(newBel);
         }
 
-        for (const auto &port : cell->ports) {
-            if (port.second.net != nullptr) {
-                auto &cost = costs[port.second.net->udata];
-                if (cost.new_cost == 0)
-                    continue;
-                cost.new_cost = 0;
-                updates.emplace_back(port.second.net);
-            }
-        }
+        ctx->bindBel(newBel, cell, STRENGTH_WEAK);
 
         if (other_cell != nullptr) {
-            for (const auto &port : other_cell->ports)
-                if (port.second.net != nullptr) {
-                    auto &cost = costs[port.second.net->udata];
-                    if (cost.new_cost == 0)
-                        continue;
-                    cost.new_cost = 0;
-                    updates.emplace_back(port.second.net);
-                }
+            ctx->bindBel(oldBel, other_cell, STRENGTH_WEAK);
         }
 
-        ctx->bindBel(newBel, cell, STRENGTH_WEAK);
+        add_move_cell(moveChange, cell, oldBel);
 
         if (other_cell != nullptr) {
-            ctx->bindBel(oldBel, other_cell, STRENGTH_WEAK);
+            add_move_cell(moveChange, other_cell, newBel);
         }
+
         if (!ctx->isBelLocationValid(newBel) || ((other_cell != nullptr && !ctx->isBelLocationValid(oldBel)))) {
             ctx->unbindBel(newBel);
             if (other_cell != nullptr)
@@ -423,26 +513,18 @@ class SAPlacer
             goto swap_fail;
         }
 
-        new_metric = curr_metric;
-
         // Recalculate metrics for all nets touched by the peturbation
-        for (const auto &net : updates) {
-            auto &c = costs[net->udata];
-            new_metric -= c.curr_cost;
-            float temp_tns = 0;
-            wirelen_t net_new_wl = get_net_metric(ctx, net, MetricType::COST, temp_tns);
-            new_metric += net_new_wl;
-            c.new_cost = net_new_wl;
-        }
+        compute_cost_changes(moveChange);
 
         new_dist = get_constraints_distance(ctx, cell);
         if (other_cell != nullptr)
             new_dist += get_constraints_distance(ctx, other_cell);
-        delta = new_metric - curr_metric;
-        delta += (cfg.constraintWeight / temp) * (new_dist - old_dist);
+        delta = lambda * (moveChange.timing_delta / std::max<double>(last_timing_cost, epsilon)) +
+                (1 - lambda) * (double(moveChange.wirelen_delta) / std::max<double>(last_wirelen_cost, epsilon));
+        delta += (cfg.constraintWeight / temp) * (new_dist - old_dist) / last_wirelen_cost;
         n_move++;
         // SA acceptance criterea
-        if (delta < 0 || (temp > 1e-6 && (ctx->rng() / float(0x3fffffff)) <= std::exp(-delta / temp))) {
+        if (delta < 0 || (temp > 1e-8 && (ctx->rng() / float(0x3fffffff)) <= std::exp(-delta / temp))) {
             n_accept++;
         } else {
             if (other_cell != nullptr)
@@ -450,32 +532,148 @@ class SAPlacer
             ctx->unbindBel(newBel);
             goto swap_fail;
         }
-        curr_metric = new_metric;
-        for (const auto &net : updates) {
-            auto &c = costs[net->udata];
-            c = CostChange{c.new_cost, -1};
-        }
-
+        commit_cost_changes(moveChange);
+#if 0
+        log_info("swap %s -> %s\n", cell->name.c_str(ctx), ctx->getBelName(newBel).c_str(ctx));
+        if (other_cell != nullptr)
+            log_info("swap %s -> %s\n", other_cell->name.c_str(ctx), ctx->getBelName(oldBel).c_str(ctx));
+#endif
         return true;
     swap_fail:
         ctx->bindBel(oldBel, cell, STRENGTH_WEAK);
         if (other_cell != nullptr) {
             ctx->bindBel(newBel, other_cell, STRENGTH_WEAK);
         }
-        for (const auto &net : updates)
-            costs[net->udata].new_cost = -1;
+        return false;
+    }
+
+    inline bool is_constrained(CellInfo *cell)
+    {
+        return cell->constr_parent != nullptr || !cell->constr_children.empty();
+    }
+
+    // Swap the Bel of a cell with another, return the original location
+    BelId swap_cell_bels(CellInfo *cell, BelId newBel)
+    {
+        BelId oldBel = cell->bel;
+#if 0
+        log_info("%s old: %s new: %s\n", cell->name.c_str(ctx), ctx->getBelName(cell->bel).c_str(ctx), ctx->getBelName(newBel).c_str(ctx));
+#endif
+        CellInfo *bound = ctx->getBoundBelCell(newBel);
+        if (bound != nullptr)
+            ctx->unbindBel(newBel);
+        ctx->unbindBel(oldBel);
+        ctx->bindBel(newBel, cell, is_constrained(cell) ? STRENGTH_STRONG : STRENGTH_WEAK);
+        if (bound != nullptr)
+            ctx->bindBel(oldBel, bound, is_constrained(bound) ? STRENGTH_STRONG : STRENGTH_WEAK);
+        return oldBel;
+    }
+
+    // Discover the relative positions of all cells in a chain
+    void discover_chain(Loc baseLoc, CellInfo *cell, std::vector<std::pair<CellInfo *, Loc>> &cell_rel)
+    {
+        Loc cellLoc = ctx->getBelLocation(cell->bel);
+        Loc rel{cellLoc.x - baseLoc.x, cellLoc.y - baseLoc.y, cellLoc.z};
+        cell_rel.emplace_back(std::make_pair(cell, rel));
+        for (auto child : cell->constr_children)
+            discover_chain(baseLoc, child, cell_rel);
+    }
+
+    // Attempt to swap a chain with a non-chain
+    bool try_swap_chain(CellInfo *cell, BelId newBase)
+    {
+        std::vector<std::pair<CellInfo *, Loc>> cell_rel;
+        std::unordered_set<IdString> cells;
+        std::vector<std::pair<CellInfo *, BelId>> moves_made;
+        std::vector<std::pair<CellInfo *, BelId>> dest_bels;
+        double delta = 0;
+        moveChange.reset();
+        if (ctx->debug)
+            log_info("finding cells for chain swap %s\n", cell->name.c_str(ctx));
+
+        Loc baseLoc = ctx->getBelLocation(cell->bel);
+        discover_chain(baseLoc, cell, cell_rel);
+        Loc newBaseLoc = ctx->getBelLocation(newBase);
+        NPNR_ASSERT(newBaseLoc.z == baseLoc.z);
+        for (const auto &cr : cell_rel)
+            cells.insert(cr.first->name);
+
+        for (const auto &cr : cell_rel) {
+            Loc targetLoc = {newBaseLoc.x + cr.second.x, newBaseLoc.y + cr.second.y, cr.second.z};
+            BelId targetBel = ctx->getBelByLocation(targetLoc);
+            if (targetBel == BelId())
+                return false;
+            if (ctx->getBelType(targetBel) != cell->type)
+                return false;
+            CellInfo *bound = ctx->getBoundBelCell(targetBel);
+            // We don't consider swapping chains with other chains, at least for the time being - unless it is
+            // part of this chain
+            if (bound != nullptr && !cells.count(bound->name) &&
+                (bound->belStrength >= STRENGTH_STRONG || is_constrained(bound)))
+                return false;
+            dest_bels.emplace_back(std::make_pair(cr.first, targetBel));
+        }
+        if (ctx->debug)
+            log_info("trying chain swap %s\n", cell->name.c_str(ctx));
+        // <cell, oldBel>
+        for (const auto &db : dest_bels) {
+            BelId oldBel = swap_cell_bels(db.first, db.second);
+            moves_made.emplace_back(std::make_pair(db.first, oldBel));
+        }
+        for (const auto &mm : moves_made) {
+            if (!ctx->isBelLocationValid(mm.first->bel) || !check_cell_bel_region(mm.first, mm.first->bel))
+                goto swap_fail;
+            if (!ctx->isBelLocationValid(mm.second))
+                goto swap_fail;
+            CellInfo *bound = ctx->getBoundBelCell(mm.second);
+            if (bound && !check_cell_bel_region(bound, bound->bel))
+                goto swap_fail;
+            add_move_cell(moveChange, mm.first, mm.second);
+            if (bound != nullptr)
+                add_move_cell(moveChange, bound, mm.first->bel);
+        }
+        compute_cost_changes(moveChange);
+        delta = lambda * (moveChange.timing_delta / last_timing_cost) +
+                (1 - lambda) * (double(moveChange.wirelen_delta) / last_wirelen_cost);
+        n_move++;
+        // SA acceptance criterea
+        if (delta < 0 || (temp > 1e-9 && (ctx->rng() / float(0x3fffffff)) <= std::exp(-delta / temp))) {
+            n_accept++;
+            if (ctx->debug)
+                log_info("accepted chain swap %s\n", cell->name.c_str(ctx));
+        } else {
+            goto swap_fail;
+        }
+        commit_cost_changes(moveChange);
+        return true;
+    swap_fail:
+        for (const auto &entry : boost::adaptors::reverse(moves_made))
+            swap_cell_bels(entry.first, entry.second);
         return false;
     }
 
     // Find a random Bel of the correct type for a cell, within the specified
     // diameter
-    BelId random_bel_for_cell(CellInfo *cell)
+    BelId random_bel_for_cell(CellInfo *cell, int force_z = -1)
     {
         IdString targetType = cell->type;
         Loc curr_loc = ctx->getBelLocation(cell->bel);
+        int count = 0;
+
+        int dx = diameter, dy = diameter;
+        if (cell->region != nullptr && cell->region->constr_bels) {
+            dx = std::min(diameter, (region_bounds[cell->region->name].x1 - region_bounds[cell->region->name].x0) + 1);
+            dy = std::min(diameter, (region_bounds[cell->region->name].y1 - region_bounds[cell->region->name].y0) + 1);
+            // Clamp location to within bounds
+            curr_loc.x = std::max(region_bounds[cell->region->name].x0, curr_loc.x);
+            curr_loc.x = std::min(region_bounds[cell->region->name].x1, curr_loc.x);
+            curr_loc.y = std::max(region_bounds[cell->region->name].y0, curr_loc.y);
+            curr_loc.y = std::min(region_bounds[cell->region->name].y1, curr_loc.y);
+        }
+
         while (true) {
-            int nx = ctx->rng(2 * diameter + 1) + std::max(curr_loc.x - diameter, 0);
-            int ny = ctx->rng(2 * diameter + 1) + std::max(curr_loc.y - diameter, 0);
+            int nx = ctx->rng(2 * dx + 1) + std::max(curr_loc.x - dx, 0);
+            int ny = ctx->rng(2 * dy + 1) + std::max(curr_loc.y - dy, 0);
             int beltype_idx, beltype_cnt;
             std::tie(beltype_idx, beltype_cnt) = bel_types.at(targetType);
             if (beltype_cnt < cfg.minBelsForGridPick)
@@ -488,41 +686,262 @@ class SAPlacer
             if (fb.size() == 0)
                 continue;
             BelId bel = fb.at(ctx->rng(int(fb.size())));
+            if (force_z != -1) {
+                Loc loc = ctx->getBelLocation(bel);
+                if (loc.z != force_z)
+                    continue;
+            }
+            if (!check_cell_bel_region(cell, bel))
+                continue;
             if (locked_bels.find(bel) != locked_bels.end())
                 continue;
+            count++;
             return bel;
         }
     }
 
+    // Return true if a net is to be entirely ignored
+    inline bool ignore_net(NetInfo *net)
+    {
+        return net->driver.cell == nullptr || net->driver.cell->bel == BelId() ||
+               ctx->getBelGlobalBuf(net->driver.cell->bel);
+    }
+
+    // Get the bounding box for a net
+    inline BoundingBox get_net_bounds(NetInfo *net)
+    {
+        BoundingBox bb;
+        NPNR_ASSERT(net->driver.cell != nullptr);
+        Loc dloc = ctx->getBelLocation(net->driver.cell->bel);
+        bb.x0 = dloc.x;
+        bb.x1 = dloc.x;
+        bb.y0 = dloc.y;
+        bb.y1 = dloc.y;
+
+        for (auto user : net->users) {
+            if (user.cell->bel == BelId())
+                continue;
+            Loc uloc = ctx->getBelLocation(user.cell->bel);
+            bb.x0 = std::min(bb.x0, uloc.x);
+            bb.x1 = std::max(bb.x1, uloc.x);
+            bb.y0 = std::min(bb.y0, uloc.y);
+            bb.y1 = std::max(bb.y1, uloc.y);
+        }
+
+        return bb;
+    }
+
+    // Get the timing cost for an arc of a net
+    inline double get_timing_cost(NetInfo *net, size_t user)
+    {
+        int cc;
+        if (net->driver.cell == nullptr)
+            return 0;
+        if (ctx->getPortTimingClass(net->driver.cell, net->driver.port, cc) == TMG_IGNORE)
+            return 0;
+        if (cfg.budgetBased) {
+            double delay = ctx->getDelayNS(ctx->predictDelay(net, net->users.at(user)));
+            return std::min(10.0, std::exp(delay - ctx->getDelayNS(net->users.at(user).budget) / 10));
+        } else {
+            auto crit = net_crit.find(net->name);
+            if (crit == net_crit.end() || crit->second.criticality.empty())
+                return 0;
+            double delay = ctx->getDelayNS(ctx->predictDelay(net, net->users.at(user)));
+            return delay * std::pow(crit->second.criticality.at(user), crit_exp);
+        }
+    }
+
+    // Set up the cost maps
+    void setup_costs()
+    {
+        for (auto net : sorted(ctx->nets)) {
+            NetInfo *ni = net.second;
+            if (ignore_net(ni))
+                continue;
+            net_bounds[ni->udata] = get_net_bounds(ni);
+            if (ctx->timing_driven && int(ni->users.size()) < cfg.timingFanoutThresh)
+                for (size_t i = 0; i < ni->users.size(); i++)
+                    net_arc_tcost[ni->udata][i] = get_timing_cost(ni, i);
+        }
+    }
+
+    // Get the total wiring cost for the design
+    wirelen_t total_wirelen_cost()
+    {
+        wirelen_t cost = 0;
+        for (const auto &net : net_bounds)
+            cost += net.hpwl();
+        return cost;
+    }
+
+    // Get the total timing cost for the design
+    double total_timing_cost()
+    {
+        double cost = 0;
+        for (const auto &net : net_arc_tcost) {
+            for (auto arc_cost : net) {
+                cost += arc_cost;
+            }
+        }
+        return cost;
+    }
+
+    // Cost-change-related data for a move
+    struct MoveChangeData
+    {
+        std::vector<decltype(NetInfo::udata)> bounds_changed_nets;
+        std::vector<std::pair<decltype(NetInfo::udata), size_t>> changed_arcs;
+
+        std::vector<bool> already_bounds_changed;
+        std::vector<std::vector<bool>> already_changed_arcs;
+
+        std::vector<std::pair<decltype(NetInfo::udata), BoundingBox>> new_net_bounds;
+        std::vector<std::pair<std::pair<decltype(NetInfo::udata), size_t>, double>> new_arc_costs;
+
+        wirelen_t wirelen_delta = 0;
+        double timing_delta = 0;
+
+        void reset()
+        {
+            for (auto bc : bounds_changed_nets)
+                already_bounds_changed[bc] = false;
+            for (const auto &tc : changed_arcs)
+                already_changed_arcs[tc.first][tc.second] = false;
+            bounds_changed_nets.clear();
+            changed_arcs.clear();
+            new_net_bounds.clear();
+            new_arc_costs.clear();
+            wirelen_delta = 0;
+            timing_delta = 0;
+        }
+
+    } moveChange;
+
+    void add_move_cell(MoveChangeData &mc, CellInfo *cell, BelId old_bel)
+    {
+        Loc curr_loc = ctx->getBelLocation(cell->bel);
+        Loc old_loc = ctx->getBelLocation(old_bel);
+        // Check net bounds
+        for (const auto &port : cell->ports) {
+            NetInfo *pn = port.second.net;
+            if (pn == nullptr)
+                continue;
+            if (ignore_net(pn))
+                continue;
+            const BoundingBox &curr_bounds = net_bounds[pn->udata];
+            // If the old location was at the edge of the bounds, or the new location exceeds the bounds,
+            // an update is needed
+            if (curr_bounds.touches_bounds(old_loc.x, old_loc.y) || !curr_bounds.is_inside_inc(curr_loc.x, curr_loc.y))
+                if (!mc.already_bounds_changed[pn->udata]) {
+                    mc.bounds_changed_nets.push_back(pn->udata);
+                    mc.already_bounds_changed[pn->udata] = true;
+                }
+            if (ctx->timing_driven && int(pn->users.size()) < cfg.timingFanoutThresh) {
+                // Output ports - all arcs change timing
+                if (port.second.type == PORT_OUT) {
+                    int cc;
+                    TimingPortClass cls = ctx->getPortTimingClass(cell, port.first, cc);
+                    if (cls != TMG_IGNORE)
+                        for (size_t i = 0; i < pn->users.size(); i++)
+                            if (!mc.already_changed_arcs[pn->udata][i]) {
+                                mc.changed_arcs.emplace_back(std::make_pair(pn->udata, i));
+                                mc.already_changed_arcs[pn->udata][i] = true;
+                            }
+                } else if (port.second.type == PORT_IN) {
+                    auto usr = fast_port_to_user.at(&port.second);
+                    if (!mc.already_changed_arcs[pn->udata][usr]) {
+                        mc.changed_arcs.emplace_back(std::make_pair(pn->udata, usr));
+                        mc.already_changed_arcs[pn->udata][usr] = true;
+                    }
+                }
+            }
+        }
+    }
+
+    void compute_cost_changes(MoveChangeData &md)
+    {
+        for (const auto &bc : md.bounds_changed_nets) {
+            wirelen_t old_hpwl = net_bounds.at(bc).hpwl();
+            auto bounds = get_net_bounds(net_by_udata.at(bc));
+            md.new_net_bounds.emplace_back(std::make_pair(bc, bounds));
+            md.wirelen_delta += (bounds.hpwl() - old_hpwl);
+            md.already_bounds_changed[bc] = false;
+        }
+        if (ctx->timing_driven) {
+            for (const auto &tc : md.changed_arcs) {
+                double old_cost = net_arc_tcost.at(tc.first).at(tc.second);
+                double new_cost = get_timing_cost(net_by_udata.at(tc.first), tc.second);
+                md.new_arc_costs.emplace_back(std::make_pair(tc, new_cost));
+                md.timing_delta += (new_cost - old_cost);
+                md.already_changed_arcs[tc.first][tc.second] = false;
+            }
+        }
+    }
+
+    void commit_cost_changes(MoveChangeData &md)
+    {
+        for (const auto &bc : md.new_net_bounds)
+            net_bounds[bc.first] = bc.second;
+        for (const auto &tc : md.new_arc_costs)
+            net_arc_tcost[tc.first.first].at(tc.first.second) = tc.second;
+        curr_wirelen_cost += md.wirelen_delta;
+        curr_timing_cost += md.timing_delta;
+    }
+    // Build the cell port -> user index
+    void build_port_index()
+    {
+        for (auto net : sorted(ctx->nets)) {
+            NetInfo *ni = net.second;
+            for (size_t i = 0; i < ni->users.size(); i++) {
+                auto &usr = ni->users.at(i);
+                fast_port_to_user[&(usr.cell->ports.at(usr.port))] = i;
+            }
+        }
+    }
+
+    // Get the combined wirelen/timing metric
+    inline double curr_metric() { return lambda * curr_timing_cost + (1 - lambda) * curr_wirelen_cost; }
+
+    // Map nets to their bounding box (so we can skip recompute for moves that do not exceed the bounds
+    std::vector<BoundingBox> net_bounds;
+    // Map net arcs to their timing cost (criticality * delay ns)
+    std::vector<std::vector<double>> net_arc_tcost;
+
+    // Fast lookup for cell port to net user index
+    std::unordered_map<const PortInfo *, size_t> fast_port_to_user;
+
+    // Wirelength and timing cost at last and current iteration
+    wirelen_t last_wirelen_cost, curr_wirelen_cost;
+    double last_timing_cost, curr_timing_cost;
+
+    // Criticality data from timing analysis
+    NetCriticalityMap net_crit;
+
     Context *ctx;
-    wirelen_t curr_metric = std::numeric_limits<wirelen_t>::max();
-    float curr_tns = 0;
-    float temp = 1000;
+    float temp = 10;
+    float crit_exp = 8;
+    float lambda = 0.5;
     bool improved = false;
     int n_move, n_accept;
     int diameter = 35, max_x = 1, max_y = 1;
     std::unordered_map<IdString, std::tuple<int, int>> bel_types;
+    std::unordered_map<IdString, BoundingBox> region_bounds;
     std::vector<std::vector<std::vector<std::vector<BelId>>>> fast_bels;
     std::unordered_set<BelId> locked_bels;
+    std::vector<NetInfo *> net_by_udata;
+    std::vector<decltype(NetInfo::udata)> old_udata;
     bool require_legal = true;
-    const float legalise_temp = 1;
-    const float post_legalise_temp = 10;
-    const float post_legalise_dia_scale = 1.5;
+    const int legalise_dia = 4;
     Placer1Cfg cfg;
-
-    struct CostChange
-    {
-        wirelen_t curr_cost;
-        wirelen_t new_cost;
-    };
-    std::vector<CostChange> costs;
-    std::vector<decltype(NetInfo::udata)> old_udata;
 };
 
 Placer1Cfg::Placer1Cfg(Context *ctx) : Settings(ctx)
 {
     constraintWeight = get<float>("placer1/constraintWeight", 10);
     minBelsForGridPick = get<int>("placer1/minBelsForGridPick", 64);
+    budgetBased = get<bool>("placer1/budgetBased", false);
+    startTemp = get<float>("placer1/startTemp", 1);
+    timingFanoutThresh = std::numeric_limits<int>::max();
 }
 
 bool placer1(Context *ctx, Placer1Cfg cfg)
@@ -545,4 +964,24 @@ bool placer1(Context *ctx, Placer1Cfg cfg)
     }
 }
 
+bool placer1_refine(Context *ctx, Placer1Cfg cfg)
+{
+    try {
+        SAPlacer placer(ctx, cfg);
+        placer.place(true);
+        log_info("Checksum: 0x%08x\n", ctx->checksum());
+#ifndef NDEBUG
+        ctx->lock();
+        ctx->check();
+        ctx->unlock();
+#endif
+        return true;
+    } catch (log_execution_error_exception) {
+#ifndef NDEBUG
+        ctx->check();
+#endif
+        return false;
+    }
+}
+
 NEXTPNR_NAMESPACE_END
diff --git a/common/placer1.h b/common/placer1.h
index 7305f4b1..4c7c7339 100644
--- a/common/placer1.h
+++ b/common/placer1.h
@@ -29,9 +29,13 @@ struct Placer1Cfg : public Settings
     Placer1Cfg(Context *ctx);
     float constraintWeight;
     int minBelsForGridPick;
+    bool budgetBased;
+    float startTemp;
+    int timingFanoutThresh;
 };
 
 extern bool placer1(Context *ctx, Placer1Cfg cfg);
+extern bool placer1_refine(Context *ctx, Placer1Cfg cfg);
 
 NEXTPNR_NAMESPACE_END
 
diff --git a/common/placer_heap.cc b/common/placer_heap.cc
new file mode 100644
index 00000000..7caf536d
--- /dev/null
+++ b/common/placer_heap.cc
@@ -0,0 +1,1544 @@
+/*
+ *  nextpnr -- Next Generation Place and Route
+ *
+ *  Copyright (C) 2019  David Shah <david@symbioticeda.com>
+ *
+ *  Permission to use, copy, modify, and/or distribute this software for any
+ *  purpose with or without fee is hereby granted, provided that the above
+ *  copyright notice and this permission notice appear in all copies.
+ *
+ *  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ *  WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ *  MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ *  ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ *  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ *  ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ *  OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ *
+ *  [[cite]] HeAP
+ *  Analytical Placement for Heterogeneous FPGAs, Marcel Gort and Jason H. Anderson
+ *  https://janders.eecg.utoronto.ca/pdfs/marcelfpl12.pdf
+ *
+ *  [[cite]] SimPL
+ *  SimPL: An Effective Placement Algorithm, Myung-Chul Kim, Dong-Jin Lee and Igor L. Markov
+ *  http://www.ece.umich.edu/cse/awards/pdfs/iccad10-simpl.pdf
+ *
+ *  Notable changes from the original algorithm
+ *   - Following the other nextpnr placer, Bels are placed rather than CLBs. This means a strict legalisation pass is
+ *     added in addition to coarse legalisation (referred to as "spreading" to avoid confusion with strict legalisation)
+ *     as described in HeAP to ensure validity. This searches random bels in the vicinity of the position chosen by
+ *     spreading, with diameter increasing over iterations, with a heuristic to prefer lower wirelength choices.
+ *   - To make the placer timing-driven, the bound2bound weights are multiplied by (1 + 10 * crit^2)
+ */
+
+#ifdef WITH_HEAP
+
+#include "placer_heap.h"
+#include <Eigen/Core>
+#include <Eigen/IterativeLinearSolvers>
+#include <boost/optional.hpp>
+#include <chrono>
+#include <deque>
+#include <fstream>
+#include <numeric>
+#include <queue>
+#include <thread>
+#include <tuple>
+#include <unordered_map>
+#include "log.h"
+#include "nextpnr.h"
+#include "place_common.h"
+#include "placer1.h"
+#include "timing.h"
+#include "util.h"
+NEXTPNR_NAMESPACE_BEGIN
+
+namespace {
+// A simple internal representation for a sparse system of equations Ax = rhs
+// This is designed to decouple the functions that build the matrix to the engine that
+// solves it, and the representation that requires
+template <typename T> struct EquationSystem
+{
+
+    EquationSystem(size_t rows, size_t cols)
+    {
+        A.resize(cols);
+        rhs.resize(rows);
+    }
+
+    // Simple sparse format, easy to convert to CCS for solver
+    std::vector<std::vector<std::pair<int, T>>> A; // col -> (row, x[row, col]) sorted by row
+    std::vector<T> rhs;                            // RHS vector
+    void reset()
+    {
+        for (auto &col : A)
+            col.clear();
+        std::fill(rhs.begin(), rhs.end(), T());
+    }
+
+    void add_coeff(int row, int col, T val)
+    {
+        auto &Ac = A.at(col);
+        // Binary search
+        int b = 0, e = int(Ac.size()) - 1;
+        while (b <= e) {
+            int i = (b + e) / 2;
+            if (Ac.at(i).first == row) {
+                Ac.at(i).second += val;
+                return;
+            }
+            if (Ac.at(i).first > row)
+                e = i - 1;
+            else
+                b = i + 1;
+        }
+        Ac.insert(Ac.begin() + b, std::make_pair(row, val));
+    }
+
+    void add_rhs(int row, T val) { rhs[row] += val; }
+
+    void solve(std::vector<T> &x)
+    {
+        using namespace Eigen;
+
+        NPNR_ASSERT(x.size() == A.size());
+
+        VectorXd vx(x.size()), vb(rhs.size());
+        SparseMatrix<T> mat(A.size(), A.size());
+
+        std::vector<int> colnnz;
+        for (auto &Ac : A)
+            colnnz.push_back(int(Ac.size()));
+        mat.reserve(colnnz);
+        for (int col = 0; col < int(A.size()); col++) {
+            auto &Ac = A.at(col);
+            for (auto &el : Ac)
+                mat.insert(el.first, col) = el.second;
+        }
+
+        for (int i = 0; i < int(x.size()); i++)
+            vx[i] = x.at(i);
+        for (int i = 0; i < int(rhs.size()); i++)
+            vb[i] = rhs.at(i);
+
+        ConjugateGradient<SparseMatrix<T>, Lower | Upper> solver;
+        solver.setTolerance(1e-5);
+        VectorXd xr = solver.compute(mat).solveWithGuess(vb, vx);
+        for (int i = 0; i < int(x.size()); i++)
+            x.at(i) = xr[i];
+        // for (int i = 0; i < int(x.size()); i++)
+        //    log_info("x[%d] = %f\n", i, x.at(i));
+    }
+};
+
+} // namespace
+
+class HeAPPlacer
+{
+  public:
+    HeAPPlacer(Context *ctx, PlacerHeapCfg cfg) : ctx(ctx), cfg(cfg) { Eigen::initParallel(); }
+
+    bool place()
+    {
+        auto startt = std::chrono::high_resolution_clock::now();
+
+        ctx->lock();
+        place_constraints();
+        build_fast_bels();
+        seed_placement();
+        update_all_chains();
+        wirelen_t hpwl = total_hpwl();
+        log_info("Creating initial analytic placement for %d cells, random placement wirelen = %d.\n",
+                 int(place_cells.size()), int(hpwl));
+        for (int i = 0; i < 4; i++) {
+            setup_solve_cells();
+            auto solve_startt = std::chrono::high_resolution_clock::now();
+            std::thread xaxis([&]() { build_solve_direction(false, -1); });
+            build_solve_direction(true, -1);
+            xaxis.join();
+            auto solve_endt = std::chrono::high_resolution_clock::now();
+            solve_time += std::chrono::duration<double>(solve_endt - solve_startt).count();
+
+            update_all_chains();
+
+            hpwl = total_hpwl();
+            log_info("    at initial placer iter %d, wirelen = %d\n", i, int(hpwl));
+        }
+
+        wirelen_t solved_hpwl = 0, spread_hpwl = 0, legal_hpwl = 0, best_hpwl = std::numeric_limits<wirelen_t>::max();
+        int iter = 0, stalled = 0;
+
+        std::vector<std::tuple<CellInfo *, BelId, PlaceStrength>> solution;
+
+        std::vector<std::unordered_set<IdString>> heap_runs;
+        std::unordered_set<IdString> all_celltypes;
+        std::unordered_map<IdString, int> ct_count;
+
+        for (auto cell : place_cells) {
+            if (!all_celltypes.count(cell->type)) {
+                heap_runs.push_back(std::unordered_set<IdString>{cell->type});
+                all_celltypes.insert(cell->type);
+            }
+            ct_count[cell->type]++;
+        }
+        // If more than 98% of cells are one cell type, always solve all at once
+        // Otherwise, follow full HeAP strategy of rotate&all
+        for (auto &c : ct_count)
+            if (c.second >= 0.98 * int(place_cells.size())) {
+                heap_runs.clear();
+                break;
+            }
+
+        heap_runs.push_back(all_celltypes);
+        // The main HeAP placer loop
+        log_info("Running main analytical placer.\n");
+        while (stalled < 5 && (solved_hpwl <= legal_hpwl * 0.8)) {
+            // Alternate between particular Bel types and all bels
+            for (auto &run : heap_runs) {
+                auto run_startt = std::chrono::high_resolution_clock::now();
+
+                setup_solve_cells(&run);
+                if (solve_cells.empty())
+                    continue;
+                // Heuristic: don't bother with threading below a certain size
+                auto solve_startt = std::chrono::high_resolution_clock::now();
+
+                if (solve_cells.size() < 500) {
+                    build_solve_direction(false, (iter == 0) ? -1 : iter);
+                    build_solve_direction(true, (iter == 0) ? -1 : iter);
+                } else {
+                    std::thread xaxis([&]() { build_solve_direction(false, (iter == 0) ? -1 : iter); });
+                    build_solve_direction(true, (iter == 0) ? -1 : iter);
+                    xaxis.join();
+                }
+                auto solve_endt = std::chrono::high_resolution_clock::now();
+                solve_time += std::chrono::duration<double>(solve_endt - solve_startt).count();
+                update_all_chains();
+                solved_hpwl = total_hpwl();
+
+                update_all_chains();
+                for (auto type : sorted(run))
+                    CutSpreader(this, type).run();
+
+                update_all_chains();
+                spread_hpwl = total_hpwl();
+                legalise_placement_strict(true);
+                update_all_chains();
+
+                legal_hpwl = total_hpwl();
+                auto run_stopt = std::chrono::high_resolution_clock::now();
+                log_info("    at iteration #%d, type %s: wirelen solved = %d, spread = %d, legal = %d; time = %.02fs\n",
+                         iter + 1, (run.size() > 1 ? "ALL" : run.begin()->c_str(ctx)), int(solved_hpwl),
+                         int(spread_hpwl), int(legal_hpwl),
+                         std::chrono::duration<double>(run_stopt - run_startt).count());
+            }
+
+            if (ctx->timing_driven)
+                get_criticalities(ctx, &net_crit);
+
+            if (legal_hpwl < best_hpwl) {
+                best_hpwl = legal_hpwl;
+                stalled = 0;
+                // Save solution
+                solution.clear();
+                for (auto cell : sorted(ctx->cells)) {
+                    solution.emplace_back(cell.second, cell.second->bel, cell.second->belStrength);
+                }
+            } else {
+                ++stalled;
+            }
+            for (auto &cl : cell_locs) {
+                cl.second.legal_x = cl.second.x;
+                cl.second.legal_y = cl.second.y;
+            }
+            ctx->yield();
+            ++iter;
+        }
+
+        // Apply saved solution
+        for (auto &sc : solution) {
+            CellInfo *cell = std::get<0>(sc);
+            if (cell->bel != BelId())
+                ctx->unbindBel(cell->bel);
+        }
+        for (auto &sc : solution) {
+            CellInfo *cell;
+            BelId bel;
+            PlaceStrength strength;
+            std::tie(cell, bel, strength) = sc;
+            ctx->bindBel(bel, cell, strength);
+        }
+
+        for (auto cell : sorted(ctx->cells)) {
+            if (cell.second->bel == BelId())
+                log_error("Found unbound cell %s\n", cell.first.c_str(ctx));
+            if (ctx->getBoundBelCell(cell.second->bel) != cell.second)
+                log_error("Found cell %s with mismatched binding\n", cell.first.c_str(ctx));
+            if (ctx->debug)
+                log_info("AP soln: %s -> %s\n", cell.first.c_str(ctx), ctx->getBelName(cell.second->bel).c_str(ctx));
+        }
+
+        ctx->unlock();
+        auto endtt = std::chrono::high_resolution_clock::now();
+        log_info("HeAP Placer Time: %.02fs\n", std::chrono::duration<double>(endtt - startt).count());
+        log_info("  of which solving equations: %.02fs\n", solve_time);
+        log_info("  of which spreading cells: %.02fs\n", cl_time);
+        log_info("  of which strict legalisation: %.02fs\n", sl_time);
+
+        ctx->check();
+
+        placer1_refine(ctx, Placer1Cfg(ctx));
+
+        return true;
+    }
+
+  private:
+    Context *ctx;
+    PlacerHeapCfg cfg;
+
+    int max_x = 0, max_y = 0;
+    std::vector<std::vector<std::vector<std::vector<BelId>>>> fast_bels;
+    std::unordered_map<IdString, std::tuple<int, int>> bel_types;
+
+    // For fast handling of heterogeneosity during initial placement without full legalisation,
+    // for each Bel type this goes from x or y to the nearest x or y where a Bel of a given type exists
+    // This is particularly important for the iCE40 architecture, where multipliers and BRAM only exist at the
+    // edges and corners respectively
+    std::vector<std::vector<int>> nearest_row_with_bel;
+    std::vector<std::vector<int>> nearest_col_with_bel;
+
+    // In some cases, we can't use bindBel because we allow overlap in the earlier stages. So we use this custom
+    // structure instead
+    struct CellLocation
+    {
+        int x, y;
+        int legal_x, legal_y;
+        double rawx, rawy;
+        bool locked, global;
+    };
+    std::unordered_map<IdString, CellLocation> cell_locs;
+    // The set of cells that we will actually place. This excludes locked cells and children cells of macros/chains
+    // (only the root of each macro is placed.)
+    std::vector<CellInfo *> place_cells;
+
+    // The cells in the current equation being solved (a subset of place_cells in some cases, where we only place
+    // cells of a certain type)
+    std::vector<CellInfo *> solve_cells;
+
+    // For cells in a chain, this is the ultimate root cell of the chain (sometimes this is not constr_parent
+    // where chains are within chains
+    std::unordered_map<IdString, CellInfo *> chain_root;
+    std::unordered_map<IdString, int> chain_size;
+
+    // The offset from chain_root to a cell in the chain
+    std::unordered_map<IdString, std::pair<int, int>> cell_offsets;
+
+    // Performance counting
+    double solve_time = 0, cl_time = 0, sl_time = 0;
+
+    NetCriticalityMap net_crit;
+
+    // Place cells with the BEL attribute set to constrain them
+    void place_constraints()
+    {
+        size_t placed_cells = 0;
+        // Initial constraints placer
+        for (auto &cell_entry : ctx->cells) {
+            CellInfo *cell = cell_entry.second.get();
+            auto loc = cell->attrs.find(ctx->id("BEL"));
+            if (loc != cell->attrs.end()) {
+                std::string loc_name = loc->second;
+                BelId bel = ctx->getBelByName(ctx->id(loc_name));
+                if (bel == BelId()) {
+                    log_error("No Bel named \'%s\' located for "
+                              "this chip (processing BEL attribute on \'%s\')\n",
+                              loc_name.c_str(), cell->name.c_str(ctx));
+                }
+
+                IdString bel_type = ctx->getBelType(bel);
+                if (bel_type != cell->type) {
+                    log_error("Bel \'%s\' of type \'%s\' does not match cell "
+                              "\'%s\' of type \'%s\'\n",
+                              loc_name.c_str(), bel_type.c_str(ctx), cell->name.c_str(ctx), cell->type.c_str(ctx));
+                }
+                if (!ctx->isValidBelForCell(cell, bel)) {
+                    log_error("Bel \'%s\' of type \'%s\' is not valid for cell "
+                              "\'%s\' of type \'%s\'\n",
+                              loc_name.c_str(), bel_type.c_str(ctx), cell->name.c_str(ctx), cell->type.c_str(ctx));
+                }
+
+                auto bound_cell = ctx->getBoundBelCell(bel);
+                if (bound_cell) {
+                    log_error("Cell \'%s\' cannot be bound to bel \'%s\' since it is already bound to cell \'%s\'\n",
+                              cell->name.c_str(ctx), loc_name.c_str(), bound_cell->name.c_str(ctx));
+                }
+
+                ctx->bindBel(bel, cell, STRENGTH_USER);
+                placed_cells++;
+            }
+        }
+        log_info("Placed %d cells based on constraints.\n", int(placed_cells));
+        ctx->yield();
+    }
+
+    // Construct the fast_bels, nearest_row_with_bel and nearest_col_with_bel
+    void build_fast_bels()
+    {
+
+        int num_bel_types = 0;
+        for (auto bel : ctx->getBels()) {
+            IdString type = ctx->getBelType(bel);
+            if (bel_types.find(type) == bel_types.end()) {
+                bel_types[type] = std::tuple<int, int>(num_bel_types++, 1);
+            } else {
+                std::get<1>(bel_types.at(type))++;
+            }
+        }
+        for (auto bel : ctx->getBels()) {
+            if (!ctx->checkBelAvail(bel))
+                continue;
+            Loc loc = ctx->getBelLocation(bel);
+            IdString type = ctx->getBelType(bel);
+            int type_idx = std::get<0>(bel_types.at(type));
+            if (int(fast_bels.size()) < type_idx + 1)
+                fast_bels.resize(type_idx + 1);
+            if (int(fast_bels.at(type_idx).size()) < (loc.x + 1))
+                fast_bels.at(type_idx).resize(loc.x + 1);
+            if (int(fast_bels.at(type_idx).at(loc.x).size()) < (loc.y + 1))
+                fast_bels.at(type_idx).at(loc.x).resize(loc.y + 1);
+            max_x = std::max(max_x, loc.x);
+            max_y = std::max(max_y, loc.y);
+            fast_bels.at(type_idx).at(loc.x).at(loc.y).push_back(bel);
+        }
+
+        nearest_row_with_bel.resize(num_bel_types, std::vector<int>(max_y + 1, -1));
+        nearest_col_with_bel.resize(num_bel_types, std::vector<int>(max_x + 1, -1));
+        for (auto bel : ctx->getBels()) {
+            if (!ctx->checkBelAvail(bel))
+                continue;
+            Loc loc = ctx->getBelLocation(bel);
+            int type_idx = std::get<0>(bel_types.at(ctx->getBelType(bel)));
+            auto &nr = nearest_row_with_bel.at(type_idx), &nc = nearest_col_with_bel.at(type_idx);
+            // Traverse outwards through nearest_row_with_bel and nearest_col_with_bel, stopping once
+            // another row/col is already recorded as being nearer
+            for (int x = loc.x; x <= max_x; x++) {
+                if (nc.at(x) != -1 && std::abs(loc.x - nc.at(x)) <= (x - loc.x))
+                    break;
+                nc.at(x) = loc.x;
+            }
+            for (int x = loc.x - 1; x >= 0; x--) {
+                if (nc.at(x) != -1 && std::abs(loc.x - nc.at(x)) <= (loc.x - x))
+                    break;
+                nc.at(x) = loc.x;
+            }
+            for (int y = loc.y; y <= max_y; y++) {
+                if (nr.at(y) != -1 && std::abs(loc.y - nr.at(y)) <= (y - loc.y))
+                    break;
+                nr.at(y) = loc.y;
+            }
+            for (int y = loc.y - 1; y >= 0; y--) {
+                if (nr.at(y) != -1 && std::abs(loc.y - nr.at(y)) <= (loc.y - y))
+                    break;
+                nr.at(y) = loc.y;
+            }
+        }
+    }
+
+    // Build and solve in one direction
+    void build_solve_direction(bool yaxis, int iter)
+    {
+        for (int i = 0; i < 5; i++) {
+            EquationSystem<double> esx(solve_cells.size(), solve_cells.size());
+            build_equations(esx, yaxis, iter);
+            solve_equations(esx, yaxis);
+        }
+    }
+
+    // Check if a cell has any meaningful connectivity
+    bool has_connectivity(CellInfo *cell)
+    {
+        for (auto port : cell->ports) {
+            if (port.second.net != nullptr && port.second.net->driver.cell != nullptr &&
+                !port.second.net->users.empty())
+                return true;
+        }
+        return false;
+    }
+
+    // Build up a random initial placement, without regard to legality
+    // FIXME: Are there better approaches to the initial placement (e.g. greedy?)
+    void seed_placement()
+    {
+        std::unordered_map<IdString, std::deque<BelId>> available_bels;
+        for (auto bel : ctx->getBels()) {
+            if (!ctx->checkBelAvail(bel))
+                continue;
+            available_bels[ctx->getBelType(bel)].push_back(bel);
+        }
+        for (auto &t : available_bels) {
+            std::random_shuffle(t.second.begin(), t.second.end(), [&](size_t n) { return ctx->rng(int(n)); });
+        }
+        for (auto cell : sorted(ctx->cells)) {
+            CellInfo *ci = cell.second;
+            if (ci->bel != BelId()) {
+                Loc loc = ctx->getBelLocation(ci->bel);
+                cell_locs[cell.first].x = loc.x;
+                cell_locs[cell.first].y = loc.y;
+                cell_locs[cell.first].locked = true;
+                cell_locs[cell.first].global = ctx->getBelGlobalBuf(ci->bel);
+            } else if (ci->constr_parent == nullptr) {
+                bool placed = false;
+                while (!placed) {
+                    if (!available_bels.count(ci->type) || available_bels.at(ci->type).empty())
+                        log_error("Unable to place cell '%s', no Bels remaining of type '%s'\n", ci->name.c_str(ctx),
+                                  ci->type.c_str(ctx));
+                    BelId bel = available_bels.at(ci->type).back();
+                    available_bels.at(ci->type).pop_back();
+                    Loc loc = ctx->getBelLocation(bel);
+                    cell_locs[cell.first].x = loc.x;
+                    cell_locs[cell.first].y = loc.y;
+                    cell_locs[cell.first].locked = false;
+                    cell_locs[cell.first].global = ctx->getBelGlobalBuf(bel);
+                    // FIXME
+                    if (has_connectivity(cell.second) && !cfg.ioBufTypes.count(ci->type)) {
+                        place_cells.push_back(ci);
+                        placed = true;
+                    } else {
+                        if (ctx->isValidBelForCell(ci, bel)) {
+                            ctx->bindBel(bel, ci, STRENGTH_STRONG);
+                            cell_locs[cell.first].locked = true;
+                            placed = true;
+                        } else {
+                            available_bels.at(ci->type).push_front(bel);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // Setup the cells to be solved, returns the number of rows
+    int setup_solve_cells(std::unordered_set<IdString> *celltypes = nullptr)
+    {
+        int row = 0;
+        solve_cells.clear();
+        // First clear the udata of all cells
+        for (auto cell : sorted(ctx->cells))
+            cell.second->udata = dont_solve;
+        // Then update cells to be placed, which excludes cell children
+        for (auto cell : place_cells) {
+            if (celltypes && !celltypes->count(cell->type))
+                continue;
+            cell->udata = row++;
+            solve_cells.push_back(cell);
+        }
+        // Finally, update the udata of children
+        for (auto chained : chain_root)
+            ctx->cells.at(chained.first)->udata = chained.second->udata;
+        return row;
+    }
+
+    // Update the location of all children of a chain
+    void update_chain(CellInfo *cell, CellInfo *root)
+    {
+        const auto &base = cell_locs[cell->name];
+        for (auto child : cell->constr_children) {
+            chain_size[root->name]++;
+            if (child->constr_x != child->UNCONSTR)
+                cell_locs[child->name].x = std::min(max_x, base.x + child->constr_x);
+            else
+                cell_locs[child->name].x = base.x; // better handling of UNCONSTR?
+            if (child->constr_y != child->UNCONSTR)
+                cell_locs[child->name].y = std::min(max_y, base.y + child->constr_y);
+            else
+                cell_locs[child->name].y = base.y; // better handling of UNCONSTR?
+            chain_root[child->name] = root;
+            if (!child->constr_children.empty())
+                update_chain(child, root);
+        }
+    }
+
+    // Update all chains
+    void update_all_chains()
+    {
+        for (auto cell : place_cells) {
+            chain_size[cell->name] = 1;
+            if (!cell->constr_children.empty())
+                update_chain(cell, cell);
+        }
+    }
+
+    // Run a function on all ports of a net - including the driver and all users
+    template <typename Tf> void foreach_port(NetInfo *net, Tf func)
+    {
+        if (net->driver.cell != nullptr)
+            func(net->driver, -1);
+        for (size_t i = 0; i < net->users.size(); i++)
+            func(net->users.at(i), i);
+    }
+
+    // Build the system of equations for either X or Y
+    void build_equations(EquationSystem<double> &es, bool yaxis, int iter = -1)
+    {
+        // Return the x or y position of a cell, depending on ydir
+        auto cell_pos = [&](CellInfo *cell) { return yaxis ? cell_locs.at(cell->name).y : cell_locs.at(cell->name).x; };
+        auto legal_pos = [&](CellInfo *cell) {
+            return yaxis ? cell_locs.at(cell->name).legal_y : cell_locs.at(cell->name).legal_x;
+        };
+
+        es.reset();
+
+        for (auto net : sorted(ctx->nets)) {
+            NetInfo *ni = net.second;
+            if (ni->driver.cell == nullptr)
+                continue;
+            if (ni->users.empty())
+                continue;
+            if (cell_locs.at(ni->driver.cell->name).global)
+                continue;
+            // Find the bounds of the net in this axis, and the ports that correspond to these bounds
+            PortRef *lbport = nullptr, *ubport = nullptr;
+            int lbpos = std::numeric_limits<int>::max(), ubpos = std::numeric_limits<int>::min();
+            foreach_port(ni, [&](PortRef &port, int user_idx) {
+                int pos = cell_pos(port.cell);
+                if (pos < lbpos) {
+                    lbpos = pos;
+                    lbport = &port;
+                }
+                if (pos > ubpos) {
+                    ubpos = pos;
+                    ubport = &port;
+                }
+            });
+            NPNR_ASSERT(lbport != nullptr);
+            NPNR_ASSERT(ubport != nullptr);
+
+            auto stamp_equation = [&](PortRef &var, PortRef &eqn, double weight) {
+                if (eqn.cell->udata == dont_solve)
+                    return;
+                int row = eqn.cell->udata;
+                int v_pos = cell_pos(var.cell);
+                if (var.cell->udata != dont_solve) {
+                    es.add_coeff(row, var.cell->udata, weight);
+                } else {
+                    es.add_rhs(row, -v_pos * weight);
+                }
+                if (cell_offsets.count(var.cell->name)) {
+                    es.add_rhs(row, -(yaxis ? cell_offsets.at(var.cell->name).second
+                                            : cell_offsets.at(var.cell->name).first) *
+                                            weight);
+                }
+            };
+
+            // Add all relevant connections to the matrix
+            foreach_port(ni, [&](PortRef &port, int user_idx) {
+                int this_pos = cell_pos(port.cell);
+                auto process_arc = [&](PortRef *other) {
+                    if (other == &port)
+                        return;
+                    int o_pos = cell_pos(other->cell);
+                    double weight = 1.0 / (ni->users.size() * std::max<double>(1, std::abs(o_pos - this_pos)));
+
+                    if (user_idx != -1 && net_crit.count(ni->name)) {
+                        auto &nc = net_crit.at(ni->name);
+                        if (user_idx < int(nc.criticality.size()))
+                            weight *= (1.0 + cfg.timingWeight *
+                                                     std::pow(nc.criticality.at(user_idx), cfg.criticalityExponent));
+                    }
+
+                    // If cell 0 is not fixed, it will stamp +w on its equation and -w on the other end's equation,
+                    // if the other end isn't fixed
+                    stamp_equation(port, port, weight);
+                    stamp_equation(port, *other, -weight);
+                    stamp_equation(*other, *other, weight);
+                    stamp_equation(*other, port, -weight);
+                };
+                process_arc(lbport);
+                process_arc(ubport);
+            });
+        }
+        if (iter != -1) {
+            float alpha = cfg.alpha;
+            for (size_t row = 0; row < solve_cells.size(); row++) {
+                int l_pos = legal_pos(solve_cells.at(row));
+                int c_pos = cell_pos(solve_cells.at(row));
+
+                double weight = alpha * iter / std::max<double>(1, std::abs(l_pos - c_pos));
+                // Add an arc from legalised to current position
+                es.add_coeff(row, row, weight);
+                es.add_rhs(row, weight * l_pos);
+            }
+        }
+    }
+
+    // Build the system of equations for either X or Y
+    void solve_equations(EquationSystem<double> &es, bool yaxis)
+    {
+        // Return the x or y position of a cell, depending on ydir
+        auto cell_pos = [&](CellInfo *cell) { return yaxis ? cell_locs.at(cell->name).y : cell_locs.at(cell->name).x; };
+        std::vector<double> vals;
+        std::transform(solve_cells.begin(), solve_cells.end(), std::back_inserter(vals), cell_pos);
+        es.solve(vals);
+        for (size_t i = 0; i < vals.size(); i++)
+            if (yaxis) {
+                cell_locs.at(solve_cells.at(i)->name).rawy = vals.at(i);
+                cell_locs.at(solve_cells.at(i)->name).y = std::min(max_y, std::max(0, int(vals.at(i))));
+            } else {
+                cell_locs.at(solve_cells.at(i)->name).rawx = vals.at(i);
+                cell_locs.at(solve_cells.at(i)->name).x = std::min(max_x, std::max(0, int(vals.at(i))));
+            }
+    }
+
+    // Compute HPWL
+    wirelen_t total_hpwl()
+    {
+        wirelen_t hpwl = 0;
+        for (auto net : sorted(ctx->nets)) {
+            NetInfo *ni = net.second;
+            if (ni->driver.cell == nullptr)
+                continue;
+            CellLocation &drvloc = cell_locs.at(ni->driver.cell->name);
+            if (drvloc.global)
+                continue;
+            int xmin = drvloc.x, xmax = drvloc.x, ymin = drvloc.y, ymax = drvloc.y;
+            for (auto &user : ni->users) {
+                CellLocation &usrloc = cell_locs.at(user.cell->name);
+                xmin = std::min(xmin, usrloc.x);
+                xmax = std::max(xmax, usrloc.x);
+                ymin = std::min(ymin, usrloc.y);
+                ymax = std::max(ymax, usrloc.y);
+            }
+            hpwl += (xmax - xmin) + (ymax - ymin);
+        }
+        return hpwl;
+    }
+
+    // Strict placement legalisation, performed after the initial HeAP spreading
+    void legalise_placement_strict(bool require_validity = false)
+    {
+        auto startt = std::chrono::high_resolution_clock::now();
+
+        // Unbind all cells placed in this solution
+        for (auto cell : sorted(ctx->cells)) {
+            CellInfo *ci = cell.second;
+            if (ci->bel != BelId() && (ci->udata != dont_solve ||
+                                       (chain_root.count(ci->name) && chain_root.at(ci->name)->udata != dont_solve)))
+                ctx->unbindBel(ci->bel);
+        }
+
+        // At the moment we don't follow the full HeAP algorithm using cuts for legalisation, instead using
+        // the simple greedy largest-macro-first approach.
+        std::priority_queue<std::pair<int, IdString>> remaining;
+        for (auto cell : solve_cells) {
+            remaining.emplace(chain_size[cell->name], cell->name);
+        }
+        int ripup_radius = 2;
+        int total_iters = 0;
+        while (!remaining.empty()) {
+            auto top = remaining.top();
+            remaining.pop();
+
+            CellInfo *ci = ctx->cells.at(top.second).get();
+            // Was now placed, ignore
+            if (ci->bel != BelId())
+                continue;
+            // log_info("   Legalising %s (%s)\n", top.second.c_str(ctx), ci->type.c_str(ctx));
+            int bt = std::get<0>(bel_types.at(ci->type));
+            auto &fb = fast_bels.at(bt);
+            int radius = 0;
+            int iter = 0;
+            int iter_at_radius = 0;
+            bool placed = false;
+            BelId bestBel;
+            int best_inp_len = std::numeric_limits<int>::max();
+
+            total_iters++;
+            if (total_iters > int(solve_cells.size())) {
+                total_iters = 0;
+                ripup_radius = std::max(std::max(max_x, max_y), ripup_radius * 2);
+            }
+
+            while (!placed) {
+
+                int nx = ctx->rng(2 * radius + 1) + std::max(cell_locs.at(ci->name).x - radius, 0);
+                int ny = ctx->rng(2 * radius + 1) + std::max(cell_locs.at(ci->name).y - radius, 0);
+
+                iter++;
+                iter_at_radius++;
+                if (iter >= (10 * (radius + 1))) {
+                    radius = std::min(std::max(max_x, max_y), radius + 1);
+                    while (radius < std::max(max_x, max_y)) {
+                        for (int x = std::max(0, cell_locs.at(ci->name).x - radius);
+                             x <= std::min(max_x, cell_locs.at(ci->name).x + radius); x++) {
+                            if (x >= int(fb.size()))
+                                break;
+                            for (int y = std::max(0, cell_locs.at(ci->name).y - radius);
+                                 y <= std::min(max_y, cell_locs.at(ci->name).y + radius); y++) {
+                                if (y >= int(fb.at(x).size()))
+                                    break;
+                                if (fb.at(x).at(y).size() > 0)
+                                    goto notempty;
+                            }
+                        }
+                        radius = std::min(std::max(max_x, max_y), radius + 1);
+                    }
+                notempty:
+                    iter_at_radius = 0;
+                    iter = 0;
+                }
+                if (nx < 0 || nx > max_x)
+                    continue;
+                if (ny < 0 || ny > max_y)
+                    continue;
+
+                // ny = nearest_row_with_bel.at(bt).at(ny);
+                // nx = nearest_col_with_bel.at(bt).at(nx);
+
+                if (nx >= int(fb.size()))
+                    continue;
+                if (ny >= int(fb.at(nx).size()))
+                    continue;
+                if (fb.at(nx).at(ny).empty())
+                    continue;
+
+                int need_to_explore = 2 * radius;
+
+                if (iter_at_radius >= need_to_explore && bestBel != BelId()) {
+                    CellInfo *bound = ctx->getBoundBelCell(bestBel);
+                    if (bound != nullptr) {
+                        ctx->unbindBel(bound->bel);
+                        remaining.emplace(chain_size[bound->name], bound->name);
+                    }
+                    ctx->bindBel(bestBel, ci, STRENGTH_WEAK);
+                    placed = true;
+                    Loc loc = ctx->getBelLocation(bestBel);
+                    cell_locs[ci->name].x = loc.x;
+                    cell_locs[ci->name].y = loc.y;
+                    break;
+                }
+
+                if (ci->constr_children.empty() && !ci->constr_abs_z) {
+                    for (auto sz : fb.at(nx).at(ny)) {
+                        if (ctx->checkBelAvail(sz) || (radius > ripup_radius || ctx->rng(20000) < 10)) {
+                            CellInfo *bound = ctx->getBoundBelCell(sz);
+                            if (bound != nullptr) {
+                                if (bound->constr_parent != nullptr || !bound->constr_children.empty() ||
+                                    bound->constr_abs_z)
+                                    continue;
+                                ctx->unbindBel(bound->bel);
+                            }
+                            ctx->bindBel(sz, ci, STRENGTH_WEAK);
+                            if (require_validity && !ctx->isBelLocationValid(sz)) {
+                                ctx->unbindBel(sz);
+                                if (bound != nullptr)
+                                    ctx->bindBel(sz, bound, STRENGTH_WEAK);
+                            } else if (iter_at_radius < need_to_explore) {
+                                ctx->unbindBel(sz);
+                                if (bound != nullptr)
+                                    ctx->bindBel(sz, bound, STRENGTH_WEAK);
+                                int input_len = 0;
+                                for (auto &port : ci->ports) {
+                                    auto &p = port.second;
+                                    if (p.type != PORT_IN || p.net == nullptr || p.net->driver.cell == nullptr)
+                                        continue;
+                                    CellInfo *drv = p.net->driver.cell;
+                                    auto drv_loc = cell_locs.find(drv->name);
+                                    if (drv_loc == cell_locs.end())
+                                        continue;
+                                    if (drv_loc->second.global)
+                                        continue;
+                                    input_len += std::abs(drv_loc->second.x - nx) + std::abs(drv_loc->second.y - ny);
+                                }
+                                if (input_len < best_inp_len) {
+                                    best_inp_len = input_len;
+                                    bestBel = sz;
+                                }
+                                break;
+                            } else {
+                                if (bound != nullptr)
+                                    remaining.emplace(chain_size[bound->name], bound->name);
+                                Loc loc = ctx->getBelLocation(sz);
+                                cell_locs[ci->name].x = loc.x;
+                                cell_locs[ci->name].y = loc.y;
+                                placed = true;
+                                break;
+                            }
+                        }
+                    }
+                } else {
+                    for (auto sz : fb.at(nx).at(ny)) {
+                        Loc loc = ctx->getBelLocation(sz);
+                        if (ci->constr_abs_z && loc.z != ci->constr_z)
+                            continue;
+                        std::vector<std::pair<CellInfo *, BelId>> targets;
+                        std::vector<std::pair<BelId, CellInfo *>> swaps_made;
+                        std::queue<std::pair<CellInfo *, Loc>> visit;
+                        visit.emplace(ci, loc);
+                        while (!visit.empty()) {
+                            CellInfo *vc = visit.front().first;
+                            NPNR_ASSERT(vc->bel == BelId());
+                            Loc ploc = visit.front().second;
+                            visit.pop();
+                            BelId target = ctx->getBelByLocation(ploc);
+                            CellInfo *bound;
+                            if (target == BelId() || ctx->getBelType(target) != vc->type)
+                                goto fail;
+                            bound = ctx->getBoundBelCell(target);
+                            // Chains cannot overlap
+                            if (bound != nullptr)
+                                if (bound->constr_z != bound->UNCONSTR || bound->constr_parent != nullptr ||
+                                    !bound->constr_children.empty() || bound->belStrength > STRENGTH_WEAK)
+                                    goto fail;
+                            targets.emplace_back(vc, target);
+                            for (auto child : vc->constr_children) {
+                                Loc cloc = ploc;
+                                if (child->constr_x != child->UNCONSTR)
+                                    cloc.x += child->constr_x;
+                                if (child->constr_y != child->UNCONSTR)
+                                    cloc.y += child->constr_y;
+                                if (child->constr_z != child->UNCONSTR)
+                                    cloc.z = child->constr_abs_z ? child->constr_z : (ploc.z + child->constr_z);
+                                visit.emplace(child, cloc);
+                            }
+                        }
+
+                        for (auto &target : targets) {
+                            CellInfo *bound = ctx->getBoundBelCell(target.second);
+                            if (bound != nullptr)
+                                ctx->unbindBel(target.second);
+                            ctx->bindBel(target.second, target.first, STRENGTH_STRONG);
+                            swaps_made.emplace_back(target.second, bound);
+                        }
+
+                        for (auto &sm : swaps_made) {
+                            if (!ctx->isBelLocationValid(sm.first))
+                                goto fail;
+                        }
+
+                        if (false) {
+                        fail:
+                            for (auto &swap : swaps_made) {
+                                ctx->unbindBel(swap.first);
+                                if (swap.second != nullptr)
+                                    ctx->bindBel(swap.first, swap.second, STRENGTH_WEAK);
+                            }
+                            continue;
+                        }
+                        for (auto &target : targets) {
+                            Loc loc = ctx->getBelLocation(target.second);
+                            cell_locs[target.first->name].x = loc.x;
+                            cell_locs[target.first->name].y = loc.y;
+                            // log_info("%s %d %d %d\n", target.first->name.c_str(ctx), loc.x, loc.y, loc.z);
+                        }
+                        for (auto &swap : swaps_made) {
+                            if (swap.second != nullptr)
+                                remaining.emplace(chain_size[swap.second->name], swap.second->name);
+                        }
+
+                        placed = true;
+                        break;
+                    }
+                }
+            }
+        }
+        auto endt = std::chrono::high_resolution_clock::now();
+        sl_time += std::chrono::duration<float>(endt - startt).count();
+    }
+    // Implementation of the cut-based spreading as described in the HeAP/SimPL papers
+    static constexpr float beta = 0.9;
+
+    struct ChainExtent
+    {
+        int x0, y0, x1, y1;
+    };
+
+    struct SpreaderRegion
+    {
+        int id;
+        int x0, y0, x1, y1;
+        int cells, bels;
+        bool overused() const
+        {
+            if (bels < 4)
+                return cells > bels;
+            else
+                return cells > beta * bels;
+        }
+    };
+
+    class CutSpreader
+    {
+      public:
+        CutSpreader(HeAPPlacer *p, IdString beltype)
+                : p(p), ctx(p->ctx), beltype(beltype), fb(p->fast_bels.at(std::get<0>(p->bel_types.at(beltype))))
+        {
+        }
+        static int seq;
+        void run()
+        {
+            auto startt = std::chrono::high_resolution_clock::now();
+            init();
+            find_overused_regions();
+            for (auto &r : regions) {
+                if (merged_regions.count(r.id))
+                    continue;
+#if 0
+                log_info("%s (%d, %d) |_> (%d, %d) %d/%d\n", beltype.c_str(ctx), r.x0, r.y0, r.x1, r.y1, r.cells,
+                         r.bels);
+#endif
+            }
+            expand_regions();
+            std::queue<std::pair<int, bool>> workqueue;
+#if 0
+            std::vector<std::pair<double, double>> orig;
+            if (ctx->debug)
+                for (auto c : p->solve_cells)
+                    orig.emplace_back(p->cell_locs[c->name].rawx, p->cell_locs[c->name].rawy);
+#endif
+            for (auto &r : regions) {
+                if (merged_regions.count(r.id))
+                    continue;
+#if 0
+                log_info("%s (%d, %d) |_> (%d, %d) %d/%d\n", beltype.c_str(ctx), r.x0, r.y0, r.x1, r.y1, r.cells,
+                         r.bels);
+#endif
+                workqueue.emplace(r.id, false);
+                // cut_region(r, false);
+            }
+            while (!workqueue.empty()) {
+                auto front = workqueue.front();
+                workqueue.pop();
+                auto &r = regions.at(front.first);
+                if (r.cells == 0)
+                    continue;
+                auto res = cut_region(r, front.second);
+                if (res) {
+                    workqueue.emplace(res->first, !front.second);
+                    workqueue.emplace(res->second, !front.second);
+                } else {
+                    // Try the other dir, in case stuck in one direction only
+                    auto res2 = cut_region(r, !front.second);
+                    if (res2) {
+                        // log_info("RETRY SUCCESS\n");
+                        workqueue.emplace(res2->first, front.second);
+                        workqueue.emplace(res2->second, front.second);
+                    }
+                }
+            }
+#if 0
+            if (ctx->debug) {
+                std::ofstream sp("spread" + std::to_string(seq) + ".csv");
+                for (size_t i = 0; i < p->solve_cells.size(); i++) {
+                    auto &c = p->solve_cells.at(i);
+                    if (c->type != beltype)
+                        continue;
+                    sp << orig.at(i).first << "," << orig.at(i).second << "," << p->cell_locs[c->name].rawx << "," << p->cell_locs[c->name].rawy << std::endl;
+                }
+                std::ofstream oc("cells" + std::to_string(seq) + ".csv");
+                for (size_t y = 0; y <= p->max_y; y++) {
+                    for (size_t x = 0; x <= p->max_x; x++) {
+                        oc << cells_at_location.at(x).at(y).size() << ", ";
+                    }
+                    oc << std::endl;
+                }
+                ++seq;
+            }
+#endif
+            auto endt = std::chrono::high_resolution_clock::now();
+            p->cl_time += std::chrono::duration<float>(endt - startt).count();
+        }
+
+      private:
+        HeAPPlacer *p;
+        Context *ctx;
+        IdString beltype;
+        std::vector<std::vector<int>> occupancy;
+        std::vector<std::vector<int>> groups;
+        std::vector<std::vector<ChainExtent>> chaines;
+        std::map<IdString, ChainExtent> cell_extents;
+
+        std::vector<std::vector<std::vector<BelId>>> &fb;
+
+        std::vector<SpreaderRegion> regions;
+        std::unordered_set<int> merged_regions;
+        // Cells at a location, sorted by real (not integer) x and y
+        std::vector<std::vector<std::vector<CellInfo *>>> cells_at_location;
+
+        int occ_at(int x, int y) { return occupancy.at(x).at(y); }
+
+        int bels_at(int x, int y)
+        {
+            if (x >= int(fb.size()) || y >= int(fb.at(x).size()))
+                return 0;
+            return int(fb.at(x).at(y).size());
+        }
+
+        void init()
+        {
+            occupancy.resize(p->max_x + 1, std::vector<int>(p->max_y + 1, 0));
+            groups.resize(p->max_x + 1, std::vector<int>(p->max_y + 1, -1));
+            chaines.resize(p->max_x + 1, std::vector<ChainExtent>(p->max_y + 1));
+            cells_at_location.resize(p->max_x + 1, std::vector<std::vector<CellInfo *>>(p->max_y + 1));
+            for (int x = 0; x <= p->max_x; x++)
+                for (int y = 0; y <= p->max_y; y++) {
+                    occupancy.at(x).at(y) = 0;
+                    groups.at(x).at(y) = -1;
+                    chaines.at(x).at(y) = {x, y, x, y};
+                }
+
+            auto set_chain_ext = [&](IdString cell, int x, int y) {
+                if (!cell_extents.count(cell))
+                    cell_extents[cell] = {x, y, x, y};
+                else {
+                    cell_extents[cell].x0 = std::min(cell_extents[cell].x0, x);
+                    cell_extents[cell].y0 = std::min(cell_extents[cell].y0, y);
+                    cell_extents[cell].x1 = std::max(cell_extents[cell].x1, x);
+                    cell_extents[cell].y1 = std::max(cell_extents[cell].y1, y);
+                }
+            };
+
+            for (auto &cell : p->cell_locs) {
+                if (ctx->cells.at(cell.first)->type != beltype)
+                    continue;
+                if (ctx->cells.at(cell.first)->belStrength > STRENGTH_STRONG)
+                    continue;
+                occupancy.at(cell.second.x).at(cell.second.y)++;
+                // Compute ultimate extent of each chain root
+                if (p->chain_root.count(cell.first)) {
+                    set_chain_ext(p->chain_root.at(cell.first)->name, cell.second.x, cell.second.y);
+                } else if (!ctx->cells.at(cell.first)->constr_children.empty()) {
+                    set_chain_ext(cell.first, cell.second.x, cell.second.y);
+                }
+            }
+            for (auto &cell : p->cell_locs) {
+                if (ctx->cells.at(cell.first)->type != beltype)
+                    continue;
+                // Transfer chain extents to the actual chaines structure
+                ChainExtent *ce = nullptr;
+                if (p->chain_root.count(cell.first))
+                    ce = &(cell_extents.at(p->chain_root.at(cell.first)->name));
+                else if (!ctx->cells.at(cell.first)->constr_children.empty())
+                    ce = &(cell_extents.at(cell.first));
+                if (ce) {
+                    auto &lce = chaines.at(cell.second.x).at(cell.second.y);
+                    lce.x0 = std::min(lce.x0, ce->x0);
+                    lce.y0 = std::min(lce.y0, ce->y0);
+                    lce.x1 = std::max(lce.x1, ce->x1);
+                    lce.y1 = std::max(lce.y1, ce->y1);
+                }
+            }
+            for (auto cell : p->solve_cells) {
+                if (cell->type != beltype)
+                    continue;
+                cells_at_location.at(p->cell_locs.at(cell->name).x).at(p->cell_locs.at(cell->name).y).push_back(cell);
+            }
+        }
+        void merge_regions(SpreaderRegion &merged, SpreaderRegion &mergee)
+        {
+            // Prevent grow_region from recursing while doing this
+            for (int x = mergee.x0; x <= mergee.x1; x++)
+                for (int y = mergee.y0; y <= mergee.y1; y++) {
+                    // log_info("%d %d\n", groups.at(x).at(y), mergee.id);
+                    NPNR_ASSERT(groups.at(x).at(y) == mergee.id);
+                    groups.at(x).at(y) = merged.id;
+                    merged.cells += occ_at(x, y);
+                    merged.bels += bels_at(x, y);
+                }
+            merged_regions.insert(mergee.id);
+            grow_region(merged, mergee.x0, mergee.y0, mergee.x1, mergee.y1);
+        }
+
+        void grow_region(SpreaderRegion &r, int x0, int y0, int x1, int y1, bool init = false)
+        {
+            // log_info("growing to (%d, %d) |_> (%d, %d)\n", x0, y0, x1, y1);
+            if ((x0 >= r.x0 && y0 >= r.y0 && x1 <= r.x1 && y1 <= r.y1) || init)
+                return;
+            int old_x0 = r.x0 + (init ? 1 : 0), old_y0 = r.y0, old_x1 = r.x1, old_y1 = r.y1;
+            r.x0 = std::min(r.x0, x0);
+            r.y0 = std::min(r.y0, y0);
+            r.x1 = std::max(r.x1, x1);
+            r.y1 = std::max(r.y1, y1);
+
+            auto process_location = [&](int x, int y) {
+                // Merge with any overlapping regions
+                if (groups.at(x).at(y) == -1) {
+                    r.bels += bels_at(x, y);
+                    r.cells += occ_at(x, y);
+                }
+                if (groups.at(x).at(y) != -1 && groups.at(x).at(y) != r.id)
+                    merge_regions(r, regions.at(groups.at(x).at(y)));
+                groups.at(x).at(y) = r.id;
+                // Grow to cover any chains
+                auto &chaine = chaines.at(x).at(y);
+                grow_region(r, chaine.x0, chaine.y0, chaine.x1, chaine.y1);
+            };
+            for (int x = r.x0; x < old_x0; x++)
+                for (int y = r.y0; y <= r.y1; y++)
+                    process_location(x, y);
+            for (int x = old_x1 + 1; x <= x1; x++)
+                for (int y = r.y0; y <= r.y1; y++)
+                    process_location(x, y);
+            for (int y = r.y0; y < old_y0; y++)
+                for (int x = r.x0; x <= r.x1; x++)
+                    process_location(x, y);
+            for (int y = old_y1 + 1; y <= r.y1; y++)
+                for (int x = r.x0; x <= r.x1; x++)
+                    process_location(x, y);
+        }
+
+        void find_overused_regions()
+        {
+            for (int x = 0; x <= p->max_x; x++)
+                for (int y = 0; y <= p->max_y; y++) {
+                    // Either already in a group, or not overutilised. Ignore
+                    if (groups.at(x).at(y) != -1 || (occ_at(x, y) <= bels_at(x, y)))
+                        continue;
+                    // log_info("%d %d %d\n", x, y, occ_at(x, y));
+                    int id = int(regions.size());
+                    groups.at(x).at(y) = id;
+                    SpreaderRegion reg;
+                    reg.id = id;
+                    reg.x0 = reg.x1 = x;
+                    reg.y0 = reg.y1 = y;
+                    reg.bels = bels_at(x, y);
+                    reg.cells = occ_at(x, y);
+                    // Make sure we cover carries, etc
+                    grow_region(reg, reg.x0, reg.y0, reg.x1, reg.y1, true);
+
+                    bool expanded = true;
+                    while (expanded) {
+                        expanded = false;
+                        // Keep trying expansion in x and y, until we find no over-occupancy cells
+                        // or hit grouped cells
+
+                        // First try expanding in x
+                        if (reg.x1 < p->max_x) {
+                            bool over_occ_x = false;
+                            for (int y1 = reg.y0; y1 <= reg.y1; y1++) {
+                                if (occ_at(reg.x1 + 1, y1) > bels_at(reg.x1 + 1, y1)) {
+                                    // log_info("(%d, %d) occ %d bels %d\n", reg.x1+ 1, y1, occ_at(reg.x1 + 1, y1),
+                                    // bels_at(reg.x1 + 1, y1));
+                                    over_occ_x = true;
+                                    break;
+                                }
+                            }
+                            if (over_occ_x) {
+                                expanded = true;
+                                grow_region(reg, reg.x0, reg.y0, reg.x1 + 1, reg.y1);
+                            }
+                        }
+
+                        if (reg.y1 < p->max_y) {
+                            bool over_occ_y = false;
+                            for (int x1 = reg.x0; x1 <= reg.x1; x1++) {
+                                if (occ_at(x1, reg.y1 + 1) > bels_at(x1, reg.y1 + 1)) {
+                                    // log_info("(%d, %d) occ %d bels %d\n", x1, reg.y1 + 1, occ_at(x1, reg.y1 + 1),
+                                    // bels_at(x1, reg.y1 + 1));
+                                    over_occ_y = true;
+                                    break;
+                                }
+                            }
+                            if (over_occ_y) {
+                                expanded = true;
+                                grow_region(reg, reg.x0, reg.y0, reg.x1, reg.y1 + 1);
+                            }
+                        }
+                    }
+                    regions.push_back(reg);
+                }
+        }
+
+        void expand_regions()
+        {
+            std::queue<int> overu_regions;
+            for (auto &r : regions) {
+                if (!merged_regions.count(r.id) && r.overused())
+                    overu_regions.push(r.id);
+            }
+            while (!overu_regions.empty()) {
+                int rid = overu_regions.front();
+                overu_regions.pop();
+                if (merged_regions.count(rid))
+                    continue;
+                auto &reg = regions.at(rid);
+                while (reg.overused()) {
+                    bool changed = false;
+                    if (reg.x0 > 0) {
+                        grow_region(reg, reg.x0 - 1, reg.y0, reg.x1, reg.y1);
+                        changed = true;
+                        if (!reg.overused())
+                            break;
+                    }
+                    if (reg.x1 < p->max_x) {
+                        grow_region(reg, reg.x0, reg.y0, reg.x1 + 1, reg.y1);
+                        changed = true;
+                        if (!reg.overused())
+                            break;
+                    }
+                    if (reg.y0 > 0) {
+                        grow_region(reg, reg.x0, reg.y0 - 1, reg.x1, reg.y1);
+                        changed = true;
+                        if (!reg.overused())
+                            break;
+                    }
+                    if (reg.y1 < p->max_y) {
+                        grow_region(reg, reg.x0, reg.y0, reg.x1, reg.y1 + 1);
+                        changed = true;
+                        if (!reg.overused())
+                            break;
+                    }
+                    if (!changed) {
+                        if (reg.cells > reg.bels)
+                            log_error("Failed to expand region (%d, %d) |_> (%d, %d) of %d %ss\n", reg.x0, reg.y0,
+                                      reg.x1, reg.y1, reg.cells, beltype.c_str(ctx));
+                        else
+                            break;
+                    }
+                }
+            }
+        }
+
+        // Implementation of the recursive cut-based spreading as described in the HeAP paper
+        // Note we use "left" to mean "-x/-y" depending on dir and "right" to mean "+x/+y" depending on dir
+
+        std::vector<CellInfo *> cut_cells;
+
+        boost::optional<std::pair<int, int>> cut_region(SpreaderRegion &r, bool dir)
+        {
+            cut_cells.clear();
+            auto &cal = cells_at_location;
+            int total_cells = 0, total_bels = 0;
+            for (int x = r.x0; x <= r.x1; x++) {
+                for (int y = r.y0; y <= r.y1; y++) {
+                    std::copy(cal.at(x).at(y).begin(), cal.at(x).at(y).end(), std::back_inserter(cut_cells));
+                    total_bels += bels_at(x, y);
+                }
+            }
+            for (auto &cell : cut_cells) {
+                total_cells += p->chain_size.count(cell->name) ? p->chain_size.at(cell->name) : 1;
+            }
+            std::sort(cut_cells.begin(), cut_cells.end(), [&](const CellInfo *a, const CellInfo *b) {
+                return dir ? (p->cell_locs.at(a->name).rawy < p->cell_locs.at(b->name).rawy)
+                           : (p->cell_locs.at(a->name).rawx < p->cell_locs.at(b->name).rawx);
+            });
+
+            if (cut_cells.size() < 2)
+                return {};
+            // Find the cells midpoint, counting chains in terms of their total size - making the initial source cut
+            int pivot_cells = 0;
+            int pivot = 0;
+            for (auto &cell : cut_cells) {
+                pivot_cells += p->chain_size.count(cell->name) ? p->chain_size.at(cell->name) : 1;
+                if (pivot_cells >= total_cells / 2)
+                    break;
+                pivot++;
+            }
+            if (pivot == int(cut_cells.size()))
+                pivot = int(cut_cells.size()) - 1;
+            // log_info("orig pivot %d lc %d rc %d\n", pivot, pivot_cells, r.cells - pivot_cells);
+
+            // Find the clearance required either side of the pivot
+            int clearance_l = 0, clearance_r = 0;
+            for (size_t i = 0; i < cut_cells.size(); i++) {
+                int size;
+                if (cell_extents.count(cut_cells.at(i)->name)) {
+                    auto &ce = cell_extents.at(cut_cells.at(i)->name);
+                    size = dir ? (ce.y1 - ce.y0 + 1) : (ce.x1 - ce.x0 + 1);
+                } else {
+                    size = 1;
+                }
+                if (int(i) < pivot)
+                    clearance_l = std::max(clearance_l, size);
+                else
+                    clearance_r = std::max(clearance_r, size);
+            }
+            // Find the target cut that minimises difference in utilisation, whilst trying to ensure that all chains
+            // still fit
+
+            // First trim the boundaries of the region in the axis-of-interest, skipping any rows/cols without any
+            // bels of the appropriate type
+            int trimmed_l = dir ? r.y0 : r.x0, trimmed_r = dir ? r.y1 : r.x1;
+            while (trimmed_l < (dir ? r.y1 : r.x1)) {
+                bool have_bels = false;
+                for (int i = dir ? r.x0 : r.y0; i <= (dir ? r.x1 : r.y1); i++)
+                    if (bels_at(dir ? i : trimmed_l, dir ? trimmed_l : i) > 0) {
+                        have_bels = true;
+                        break;
+                    }
+                if (have_bels)
+                    break;
+                trimmed_l++;
+            }
+            while (trimmed_r > (dir ? r.y0 : r.x0)) {
+                bool have_bels = false;
+                for (int i = dir ? r.x0 : r.y0; i <= (dir ? r.x1 : r.y1); i++)
+                    if (bels_at(dir ? i : trimmed_r, dir ? trimmed_r : i) > 0) {
+                        have_bels = true;
+                        break;
+                    }
+                if (have_bels)
+                    break;
+                trimmed_r--;
+            }
+            // log_info("tl %d tr %d cl %d cr %d\n", trimmed_l, trimmed_r, clearance_l, clearance_r);
+            if ((trimmed_r - trimmed_l + 1) <= std::max(clearance_l, clearance_r))
+                return {};
+            // Now find the initial target cut that minimises utilisation imbalance, whilst
+            // meeting the clearance requirements for any large macros
+            int left_cells = pivot_cells, right_cells = total_cells - pivot_cells;
+            int left_bels = 0, right_bels = total_bels;
+            int best_tgt_cut = -1;
+            double best_deltaU = std::numeric_limits<double>::max();
+            std::pair<int, int> target_cut_bels;
+            for (int i = trimmed_l; i <= trimmed_r; i++) {
+                int slither_bels = 0;
+                for (int j = dir ? r.x0 : r.y0; j <= (dir ? r.x1 : r.y1); j++) {
+                    slither_bels += dir ? bels_at(j, i) : bels_at(i, j);
+                }
+                left_bels += slither_bels;
+                right_bels -= slither_bels;
+                if (((i - trimmed_l) + 1) >= clearance_l && ((trimmed_r - i) + 1) >= clearance_r) {
+                    // Solution is potentially valid
+                    double aU =
+                            std::abs(double(left_cells) / double(left_bels) - double(right_cells) / double(right_bels));
+                    if (aU < best_deltaU) {
+                        best_deltaU = aU;
+                        best_tgt_cut = i;
+                        target_cut_bels = std::make_pair(left_bels, right_bels);
+                    }
+                }
+            }
+            if (best_tgt_cut == -1)
+                return {};
+            left_bels = target_cut_bels.first;
+            right_bels = target_cut_bels.second;
+            // log_info("pivot %d target cut %d lc %d lb %d rc %d rb %d\n", pivot, best_tgt_cut, left_cells, left_bels,
+            // right_cells, right_bels);
+
+            // Peturb the source cut to eliminate overutilisation
+            while (pivot > 0 && (double(left_cells) / double(left_bels) > double(right_cells) / double(right_bels))) {
+                auto &move_cell = cut_cells.at(pivot);
+                int size = p->chain_size.count(move_cell->name) ? p->chain_size.at(move_cell->name) : 1;
+                left_cells -= size;
+                right_cells += size;
+                pivot--;
+            }
+            while (pivot < int(cut_cells.size()) - 1 &&
+                   (double(left_cells) / double(left_bels) < double(right_cells) / double(right_bels))) {
+                auto &move_cell = cut_cells.at(pivot + 1);
+                int size = p->chain_size.count(move_cell->name) ? p->chain_size.at(move_cell->name) : 1;
+                left_cells += size;
+                right_cells -= size;
+                pivot++;
+            }
+            // log_info("peturbed pivot %d lc %d lb %d rc %d rb %d\n", pivot, left_cells, left_bels, right_cells,
+            // right_bels);
+            // Split regions into bins, and then spread cells by linear interpolation within those bins
+            auto spread_binlerp = [&](int cells_start, int cells_end, double area_l, double area_r) {
+                int N = cells_end - cells_start;
+                if (N <= 2) {
+                    for (int i = cells_start; i < cells_end; i++) {
+                        auto &pos = dir ? p->cell_locs.at(cut_cells.at(i)->name).rawy
+                                        : p->cell_locs.at(cut_cells.at(i)->name).rawx;
+                        pos = area_l + i * ((area_r - area_l) / N);
+                    }
+                    return;
+                }
+                // Split region into up to 10 (K) bins
+                int K = std::min<int>(N, 10);
+                std::vector<std::pair<int, double>> bin_bounds; // [(cell start, area start)]
+                bin_bounds.emplace_back(cells_start, area_l);
+                for (int i = 1; i < K; i++)
+                    bin_bounds.emplace_back(cells_start + (N * i) / K, area_l + ((area_r - area_l + 0.99) * i) / K);
+                bin_bounds.emplace_back(cells_end, area_r + 0.99);
+                for (int i = 0; i < K; i++) {
+                    auto &bl = bin_bounds.at(i), br = bin_bounds.at(i + 1);
+                    double orig_left = dir ? p->cell_locs.at(cut_cells.at(bl.first)->name).rawy
+                                           : p->cell_locs.at(cut_cells.at(bl.first)->name).rawx;
+                    double orig_right = dir ? p->cell_locs.at(cut_cells.at(br.first - 1)->name).rawy
+                                            : p->cell_locs.at(cut_cells.at(br.first - 1)->name).rawx;
+                    double m = (br.second - bl.second) / std::max(0.00001, orig_right - orig_left);
+                    for (int j = bl.first; j < br.first; j++) {
+                        auto &pos = dir ? p->cell_locs.at(cut_cells.at(j)->name).rawy
+                                        : p->cell_locs.at(cut_cells.at(j)->name).rawx;
+                        NPNR_ASSERT(pos >= orig_left && pos <= orig_right);
+                        pos = bl.second + m * (pos - orig_left);
+                        // log("[%f, %f] -> [%f, %f]: %f -> %f\n", orig_left, orig_right, bl.second, br.second,
+                        // orig_pos, pos);
+                    }
+                }
+            };
+            spread_binlerp(0, pivot + 1, trimmed_l, best_tgt_cut);
+            spread_binlerp(pivot + 1, int(cut_cells.size()), best_tgt_cut + 1, trimmed_r);
+            // Update various data structures
+            for (int x = r.x0; x <= r.x1; x++)
+                for (int y = r.y0; y <= r.y1; y++) {
+                    cells_at_location.at(x).at(y).clear();
+                }
+            for (auto cell : cut_cells) {
+                auto &cl = p->cell_locs.at(cell->name);
+                cl.x = std::min(r.x1, std::max(r.x0, int(cl.rawx)));
+                cl.y = std::min(r.y1, std::max(r.y0, int(cl.rawy)));
+                cells_at_location.at(cl.x).at(cl.y).push_back(cell);
+                // log_info("spread pos %d %d\n", cl.x, cl.y);
+            }
+            SpreaderRegion rl, rr;
+            rl.id = int(regions.size());
+            rl.x0 = r.x0;
+            rl.y0 = r.y0;
+            rl.x1 = dir ? r.x1 : best_tgt_cut;
+            rl.y1 = dir ? best_tgt_cut : r.y1;
+            rl.cells = left_cells;
+            rl.bels = left_bels;
+            rr.id = int(regions.size()) + 1;
+            rr.x0 = dir ? r.x0 : (best_tgt_cut + 1);
+            rr.y0 = dir ? (best_tgt_cut + 1) : r.y0;
+            rr.x1 = r.x1;
+            rr.y1 = r.y1;
+            rr.cells = right_cells;
+            rr.bels = right_bels;
+            regions.push_back(rl);
+            regions.push_back(rr);
+            for (int x = rl.x0; x <= rl.x1; x++)
+                for (int y = rl.y0; y <= rl.y1; y++)
+                    groups.at(x).at(y) = rl.id;
+            for (int x = rr.x0; x <= rr.x1; x++)
+                for (int y = rr.y0; y <= rr.y1; y++)
+                    groups.at(x).at(y) = rr.id;
+            return std::make_pair(rl.id, rr.id);
+        };
+    };
+    typedef decltype(CellInfo::udata) cell_udata_t;
+    cell_udata_t dont_solve = std::numeric_limits<cell_udata_t>::max();
+};
+int HeAPPlacer::CutSpreader::seq = 0;
+
+bool placer_heap(Context *ctx, PlacerHeapCfg cfg) { return HeAPPlacer(ctx, cfg).place(); }
+
+PlacerHeapCfg::PlacerHeapCfg(Context *ctx) : Settings(ctx)
+{
+    alpha = get<float>("placerHeap/alpha", 0.1);
+    criticalityExponent = get<int>("placerHeap/criticalityExponent", 2);
+    timingWeight = get<int>("placerHeap/timingWeight", 10);
+}
+
+NEXTPNR_NAMESPACE_END
+
+#else
+
+#include "log.h"
+#include "nextpnr.h"
+#include "placer_heap.h"
+
+NEXTPNR_NAMESPACE_BEGIN
+bool placer_heap(Context *ctx, PlacerHeapCfg cfg)
+{
+    log_error("nextpnr was built without the HeAP placer\n");
+    return false;
+}
+
+PlacerHeapCfg::PlacerHeapCfg(Context *ctx) : Settings(ctx) {}
+
+NEXTPNR_NAMESPACE_END
+
+#endif
diff --git a/common/placer_heap.h b/common/placer_heap.h
new file mode 100644
index 00000000..f94489c1
--- /dev/null
+++ b/common/placer_heap.h
@@ -0,0 +1,47 @@
+/*
+ *  nextpnr -- Next Generation Place and Route
+ *
+ *  Copyright (C) 2019  David Shah <david@symbioticeda.com>
+ *
+ *  Permission to use, copy, modify, and/or distribute this software for any
+ *  purpose with or without fee is hereby granted, provided that the above
+ *  copyright notice and this permission notice appear in all copies.
+ *
+ *  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ *  WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ *  MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ *  ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ *  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ *  ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ *  OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ *
+ *  [[cite]] HeAP
+ *  Analytical Placement for Heterogeneous FPGAs, Marcel Gort and Jason H. Anderson
+ *  https://janders.eecg.utoronto.ca/pdfs/marcelfpl12.pdf
+ *
+ *  [[cite]] SimPL
+ *  SimPL: An Effective Placement Algorithm, Myung-Chul Kim, Dong-Jin Lee and Igor L. Markov
+ *  http://www.ece.umich.edu/cse/awards/pdfs/iccad10-simpl.pdf
+ */
+
+#ifndef PLACER_HEAP_H
+#define PLACER_HEAP
+#include "nextpnr.h"
+#include "settings.h"
+
+NEXTPNR_NAMESPACE_BEGIN
+
+struct PlacerHeapCfg : public Settings
+{
+    PlacerHeapCfg(Context *ctx);
+
+    float alpha;
+    float criticalityExponent;
+    float timingWeight;
+
+    std::unordered_set<IdString> ioBufTypes;
+};
+
+extern bool placer_heap(Context *ctx, PlacerHeapCfg cfg);
+NEXTPNR_NAMESPACE_END
+#endif
diff --git a/common/pybindings.cc b/common/pybindings.cc
index 6cae889d..eee78b5e 100644
--- a/common/pybindings.cc
+++ b/common/pybindings.cc
@@ -104,6 +104,7 @@ BOOST_PYTHON_MODULE(MODULE_NAME)
     typedef std::unordered_map<IdString, std::string> AttrMap;
     typedef std::unordered_map<IdString, PortInfo> PortMap;
     typedef std::unordered_map<IdString, IdString> PinMap;
+    typedef std::unordered_map<IdString, std::unique_ptr<Region>> RegionMap;
 
     class_<BaseCtx, BaseCtx *, boost::noncopyable>("BaseCtx", no_init);
 
@@ -135,6 +136,8 @@ BOOST_PYTHON_MODULE(MODULE_NAME)
 
     typedef std::vector<PortRef> PortRefVector;
     typedef std::unordered_map<WireId, PipMap> WireMap;
+    typedef std::unordered_set<BelId> BelSet;
+    typedef std::unordered_set<WireId> WireSet;
 
     auto ni_cls = class_<ContextualWrapper<NetInfo &>>("NetInfo", no_init);
     readwrite_wrapper<NetInfo &, decltype(&NetInfo::name), &NetInfo::name, conv_to_str<IdString>,
@@ -163,10 +166,25 @@ BOOST_PYTHON_MODULE(MODULE_NAME)
     def("parse_json", parse_json_shim);
     def("load_design", load_design_shim, return_value_policy<manage_new_object>());
 
+    auto region_cls = class_<ContextualWrapper<Region &>>("Region", no_init);
+    readwrite_wrapper<Region &, decltype(&Region::name), &Region::name, conv_to_str<IdString>,
+                      conv_from_str<IdString>>::def_wrap(region_cls, "name");
+    readwrite_wrapper<Region &, decltype(&Region::constr_bels), &Region::constr_bels, pass_through<bool>,
+                      pass_through<bool>>::def_wrap(region_cls, "constr_bels");
+    readwrite_wrapper<Region &, decltype(&Region::constr_wires), &Region::constr_wires, pass_through<bool>,
+                      pass_through<bool>>::def_wrap(region_cls, "constr_bels");
+    readwrite_wrapper<Region &, decltype(&Region::constr_pips), &Region::constr_pips, pass_through<bool>,
+                      pass_through<bool>>::def_wrap(region_cls, "constr_pips");
+    readonly_wrapper<Region &, decltype(&Region::bels), &Region::bels, wrap_context<BelSet &>>::def_wrap(region_cls,
+                                                                                                         "bels");
+    readonly_wrapper<Region &, decltype(&Region::wires), &Region::wires, wrap_context<WireSet &>>::def_wrap(region_cls,
+                                                                                                            "wires");
+
     WRAP_MAP(AttrMap, pass_through<std::string>, "AttrMap");
     WRAP_MAP(PortMap, wrap_context<PortInfo &>, "PortMap");
     WRAP_MAP(PinMap, conv_to_str<IdString>, "PinMap");
     WRAP_MAP(WireMap, wrap_context<PipMap &>, "WireMap");
+    WRAP_MAP_UPTR(RegionMap, "RegionMap");
 
     WRAP_VECTOR(PortRefVector, wrap_context<PortRef &>);
 
diff --git a/common/pycontainers.h b/common/pycontainers.h
index 70f69c51..5de2f6d2 100644
--- a/common/pycontainers.h
+++ b/common/pycontainers.h
@@ -345,6 +345,12 @@ template <typename T, typename value_conv> struct map_wrapper
         std::terminate();
     }
 
+    static bool contains(wrapped_map &x, std::string const &i)
+    {
+        K k = PythonConversion::string_converter<K>().from_str(x.ctx, i);
+        return x.base.count(k);
+    }
+
     static void wrap(const char *map_name, const char *kv_name, const char *kv_iter_name, const char *iter_name)
     {
         map_pair_wrapper<typename KV::first_type, typename KV::second_type, value_conv>::wrap(kv_name, kv_iter_name);
@@ -353,6 +359,7 @@ template <typename T, typename value_conv> struct map_wrapper
         class_<wrapped_map>(map_name, no_init)
                 .def("__iter__", rw::iter)
                 .def("__len__", len)
+                .def("__contains__", contains)
                 .def("__getitem__", get)
                 .def("__setitem__", set, with_custodian_and_ward<1, 2>());
     }
@@ -465,6 +472,12 @@ template <typename T> struct map_wrapper_uptr
         std::terminate();
     }
 
+    static bool contains(wrapped_map &x, std::string const &i)
+    {
+        K k = PythonConversion::string_converter<K>().from_str(x.ctx, i);
+        return x.base.count(k);
+    }
+
     static void wrap(const char *map_name, const char *kv_name, const char *kv_iter_name, const char *iter_name)
     {
         map_pair_wrapper_uptr<typename KV::first_type, typename KV::second_type>::wrap(kv_name, kv_iter_name);
@@ -473,6 +486,7 @@ template <typename T> struct map_wrapper_uptr
         class_<wrapped_map>(map_name, no_init)
                 .def("__iter__", rw::iter)
                 .def("__len__", len)
+                .def("__contains__", contains)
                 .def("__getitem__", get)
                 .def("__setitem__", set, with_custodian_and_ward<1, 2>());
     }
diff --git a/common/pywrappers.h b/common/pywrappers.h
index 4e463afd..427c3623 100644
--- a/common/pywrappers.h
+++ b/common/pywrappers.h
@@ -269,7 +269,7 @@ template <typename Class, typename FuncT, FuncT fn, typename arg1_conv, typename
     template <typename WrapCls> static void def_wrap(WrapCls cls_, const char *name) { cls_.def(name, wrapped_fn); }
 };
 
-// Three parameters, one return
+// Three parameters, no return
 template <typename Class, typename FuncT, FuncT fn, typename arg1_conv, typename arg2_conv, typename arg3_conv>
 struct fn_wrapper_3a_v
 {
@@ -288,6 +288,30 @@ struct fn_wrapper_3a_v
     template <typename WrapCls> static void def_wrap(WrapCls cls_, const char *name) { cls_.def(name, wrapped_fn); }
 };
 
+// Five parameters, no return
+template <typename Class, typename FuncT, FuncT fn, typename arg1_conv, typename arg2_conv, typename arg3_conv,
+          typename arg4_conv, typename arg5_conv>
+struct fn_wrapper_5a_v
+{
+    using class_type = typename WrapIfNotContext<Class>::maybe_wrapped_t;
+    using conv_arg1_type = typename arg1_conv::arg_type;
+    using conv_arg2_type = typename arg2_conv::arg_type;
+    using conv_arg3_type = typename arg3_conv::arg_type;
+    using conv_arg4_type = typename arg4_conv::arg_type;
+    using conv_arg5_type = typename arg5_conv::arg_type;
+
+    static void wrapped_fn(class_type &cls, conv_arg1_type arg1, conv_arg2_type arg2, conv_arg3_type arg3,
+                           conv_arg4_type arg4, conv_arg5_type arg5)
+    {
+        Context *ctx = get_ctx<Class>(cls);
+        Class &base = get_base<Class>(cls);
+        return (base.*fn)(arg1_conv()(ctx, arg1), arg2_conv()(ctx, arg2), arg3_conv()(ctx, arg3),
+                          arg4_conv()(ctx, arg4), arg5_conv()(ctx, arg5));
+    }
+
+    template <typename WrapCls> static void def_wrap(WrapCls cls_, const char *name) { cls_.def(name, wrapped_fn); }
+};
+
 // Wrapped getter
 template <typename Class, typename MemT, MemT mem, typename v_conv> struct readonly_wrapper
 {
diff --git a/common/settings.h b/common/settings.h
index 0c4a67db..b57947c9 100644
--- a/common/settings.h
+++ b/common/settings.h
@@ -45,19 +45,30 @@ class Settings
         return defaultValue;
     }
 
-    template <typename T> void set(const char *name, T value)
-    {
-        IdString id = ctx->id(name);
-        auto pair = ctx->settings.emplace(id, std::to_string(value));
-        if (!pair.second) {
-            ctx->settings[pair.first->first] = value;
-        }
-    }
+    template <typename T> void set(const char *name, T value);
 
   private:
     Context *ctx;
 };
 
+template <typename T> inline void Settings::set(const char *name, T value)
+{
+    IdString id = ctx->id(name);
+    auto pair = ctx->settings.emplace(id, std::to_string(value));
+    if (!pair.second) {
+        ctx->settings[pair.first->first] = value;
+    }
+}
+
+template <> inline void Settings::set<std::string>(const char *name, std::string value)
+{
+    IdString id = ctx->id(name);
+    auto pair = ctx->settings.emplace(id, value);
+    if (!pair.second) {
+        ctx->settings[pair.first->first] = value;
+    }
+}
+
 NEXTPNR_NAMESPACE_END
 
 #endif // SETTINGS_H
diff --git a/common/timing.cc b/common/timing.cc
index 2a0af874..2ce9eea3 100644
--- a/common/timing.cc
+++ b/common/timing.cc
@@ -904,10 +904,9 @@ void timing_analysis(Context *ctx, bool print_histogram, bool print_fmax, bool p
             if (!warn_on_failure || passed)
                 log_info("Max frequency for clock %*s'%s': %.02f MHz (%s at %.02f MHz)\n", width, "",
                          clock_name.c_str(), clock_fmax[clock.first], passed ? "PASS" : "FAIL", target);
-            else
-            if (bool_or_default(ctx->settings, ctx->id("timing/allowFail"), false))
+            else if (bool_or_default(ctx->settings, ctx->id("timing/allowFail"), false))
                 log_warning("Max frequency for clock %*s'%s': %.02f MHz (%s at %.02f MHz)\n", width, "",
-                         clock_name.c_str(), clock_fmax[clock.first], passed ? "PASS" : "FAIL", target);
+                            clock_name.c_str(), clock_fmax[clock.first], passed ? "PASS" : "FAIL", target);
             else
                 log_nonfatal_error("Max frequency for clock %*s'%s': %.02f MHz (%s at %.02f MHz)\n", width, "",
                                    clock_name.c_str(), clock_fmax[clock.first], passed ? "PASS" : "FAIL", target);
diff --git a/docs/archapi.md b/docs/archapi.md
index 3c938865..6e59ecdb 100644
--- a/docs/archapi.md
+++ b/docs/archapi.md
@@ -490,3 +490,14 @@ a certain number of different clock signals allowed for a group of bels.
 
 Returns true if a bell in the current configuration is valid, i.e. if
 `isValidBelForCell()` would return true for the current mapping.
+
+
+### static const std::string defaultPlacer
+
+Name of the default placement algorithm for the architecture, if
+`--placer` isn't specified on the command line.
+
+### static const std::vector\<std::string\> availablePlacers
+
+Name of available placer algorithms for the architecture, used
+to provide help for and validate `--placer`.
\ No newline at end of file
diff --git a/ecp5/arch.cc b/ecp5/arch.cc
index da0f7b1a..9da8abdf 100644
--- a/ecp5/arch.cc
+++ b/ecp5/arch.cc
@@ -28,6 +28,7 @@
 #include "log.h"
 #include "nextpnr.h"
 #include "placer1.h"
+#include "placer_heap.h"
 #include "router1.h"
 #include "timing.h"
 #include "util.h"
@@ -456,6 +457,7 @@ delay_t Arch::estimateDelay(WireId src, WireId dst) const
     auto src_loc = est_location(src), dst_loc = est_location(dst);
 
     int dx = abs(src_loc.first - dst_loc.first), dy = abs(src_loc.second - dst_loc.second);
+
     return (130 - 25 * args.speed) *
            (6 + std::max(dx - 5, 0) + std::max(dy - 5, 0) + 2 * (std::min(dx, 5) + std::min(dy, 5)));
 }
@@ -467,7 +469,6 @@ delay_t Arch::predictDelay(const NetInfo *net_info, const PortRef &sink) const
         return 0;
     auto driver_loc = getBelLocation(driver.cell->bel);
     auto sink_loc = getBelLocation(sink.cell->bel);
-
     // Encourage use of direct interconnect
     if (driver_loc.x == sink_loc.x && driver_loc.y == sink_loc.y) {
         if ((sink.port == id_A0 || sink.port == id_A1) && (driver.port == id_F1) &&
@@ -485,6 +486,7 @@ delay_t Arch::predictDelay(const NetInfo *net_info, const PortRef &sink) const
     }
 
     int dx = abs(driver_loc.x - sink_loc.x), dy = abs(driver_loc.y - sink_loc.y);
+
     return (130 - 25 * args.speed) *
            (6 + std::max(dx - 5, 0) + std::max(dy - 5, 0) + 2 * (std::min(dx, 5) + std::min(dy, 5)));
 }
@@ -506,10 +508,23 @@ bool Arch::getBudgetOverride(const NetInfo *net_info, const PortRef &sink, delay
 
 bool Arch::place()
 {
-    bool result = placer1(getCtx(), Placer1Cfg(getCtx()));
-    if (result)
-        permute_luts();
-    return result;
+    std::string placer = str_or_default(settings, id("placer"), defaultPlacer);
+
+    if (placer == "heap") {
+        PlacerHeapCfg cfg(getCtx());
+        cfg.criticalityExponent = 7;
+        cfg.ioBufTypes.insert(id_TRELLIS_IO);
+        if (!placer_heap(getCtx(), cfg))
+            return false;
+    } else if (placer == "sa") {
+        if (!placer1(getCtx(), Placer1Cfg(getCtx())))
+            return false;
+    } else {
+        log_error("ECP5 architecture does not support placer '%s'\n", placer.c_str());
+    }
+
+    permute_luts();
+    return true;
 }
 
 bool Arch::route()
@@ -605,6 +620,11 @@ DecalXY Arch::getGroupDecal(GroupId pip) const { return {}; };
 
 bool Arch::getDelayFromTimingDatabase(IdString tctype, IdString from, IdString to, DelayInfo &delay) const
 {
+    auto fnd_dk = celldelay_cache.find({tctype, from, to});
+    if (fnd_dk != celldelay_cache.end()) {
+        delay = fnd_dk->second.second;
+        return fnd_dk->second.first;
+    }
     for (int i = 0; i < speed_grade->num_cell_timings; i++) {
         const auto &tc = speed_grade->cell_timings[i];
         if (tc.cell_type == tctype.index) {
@@ -613,9 +633,11 @@ bool Arch::getDelayFromTimingDatabase(IdString tctype, IdString from, IdString t
                 if (dly.from_port == from.index && dly.to_port == to.index) {
                     delay.max_delay = dly.max_delay;
                     delay.min_delay = dly.min_delay;
+                    celldelay_cache[{tctype, from, to}] = std::make_pair(true, delay);
                     return true;
                 }
             }
+            celldelay_cache[{tctype, from, to}] = std::make_pair(false, DelayInfo());
             return false;
         }
     }
@@ -645,7 +667,6 @@ void Arch::getSetupHoldFromTimingDatabase(IdString tctype, IdString clock, IdStr
 
 bool Arch::getCellDelay(const CellInfo *cell, IdString fromPort, IdString toPort, DelayInfo &delay) const
 {
-
     // Data for -8 grade
     if (cell->type == id_TRELLIS_SLICE) {
         bool has_carry = cell->sliceInfo.is_carry;
@@ -965,4 +986,16 @@ WireId Arch::getBankECLK(int bank, int eclk)
     return getWireByLocAndBasename(Location(0, 0), "G_BANK" + std::to_string(bank) + "ECLK" + std::to_string(eclk));
 }
 
+#ifdef WITH_HEAP
+const std::string Arch::defaultPlacer = "heap";
+#else
+const std::string Arch::defaultPlacer = "sa";
+#endif
+
+const std::vector<std::string> Arch::availablePlacers = {"sa",
+#ifdef WITH_HEAP
+                                                         "heap"
+#endif
+};
+
 NEXTPNR_NAMESPACE_END
diff --git a/ecp5/arch.h b/ecp5/arch.h
index ab4a4e00..3de06a42 100644
--- a/ecp5/arch.h
+++ b/ecp5/arch.h
@@ -448,6 +448,30 @@ struct ArchArgs
     } speed = SPEED_6;
 };
 
+struct DelayKey
+{
+    IdString celltype, from, to;
+    inline bool operator==(const DelayKey &other) const
+    {
+        return celltype == other.celltype && from == other.from && to == other.to;
+    }
+};
+
+NEXTPNR_NAMESPACE_END
+namespace std {
+template <> struct hash<NEXTPNR_NAMESPACE_PREFIX DelayKey>
+{
+    std::size_t operator()(const NEXTPNR_NAMESPACE_PREFIX DelayKey &dk) const noexcept
+    {
+        std::size_t seed = std::hash<NEXTPNR_NAMESPACE_PREFIX IdString>()(dk.celltype);
+        seed ^= std::hash<NEXTPNR_NAMESPACE_PREFIX IdString>()(dk.from) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+        seed ^= std::hash<NEXTPNR_NAMESPACE_PREFIX IdString>()(dk.to) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+        return seed;
+    }
+};
+} // namespace std
+NEXTPNR_NAMESPACE_BEGIN
+
 struct Arch : BaseCtx
 {
     const ChipInfoPOD *chip_info;
@@ -1019,6 +1043,11 @@ struct Arch : BaseCtx
     IdString id_clk, id_lsr;
     IdString id_clkmux, id_lsrmux;
     IdString id_srmode, id_mode;
+
+    mutable std::unordered_map<DelayKey, std::pair<bool, DelayInfo>> celldelay_cache;
+
+    static const std::string defaultPlacer;
+    static const std::vector<std::string> availablePlacers;
 };
 
 NEXTPNR_NAMESPACE_END
diff --git a/ecp5/arch_pybindings.cc b/ecp5/arch_pybindings.cc
index 5e73a673..18d21112 100644
--- a/ecp5/arch_pybindings.cc
+++ b/ecp5/arch_pybindings.cc
@@ -133,6 +133,13 @@ void arch_wrap_python()
 
     fn_wrapper_2a_v<Context, decltype(&Context::addClock), &Context::addClock, conv_from_str<IdString>,
                     pass_through<float>>::def_wrap(ctx_cls, "addClock");
+    fn_wrapper_5a_v<Context, decltype(&Context::createRectangularRegion), &Context::createRectangularRegion,
+                    conv_from_str<IdString>, pass_through<int>, pass_through<int>, pass_through<int>,
+                    pass_through<int>>::def_wrap(ctx_cls, "createRectangularRegion");
+    fn_wrapper_2a_v<Context, decltype(&Context::addBelToRegion), &Context::addBelToRegion, conv_from_str<IdString>,
+                    conv_from_str<BelId>>::def_wrap(ctx_cls, "addBelToRegion");
+    fn_wrapper_2a_v<Context, decltype(&Context::constrainCellToRegion), &Context::constrainCellToRegion,
+                    conv_from_str<IdString>, conv_from_str<IdString>>::def_wrap(ctx_cls, "constrainCellToRegion");
 
     WRAP_RANGE(Bel, conv_to_str<BelId>);
     WRAP_RANGE(Wire, conv_to_str<WireId>);
diff --git a/ecp5/main.cc b/ecp5/main.cc
index 15027a5a..bb18aa58 100644
--- a/ecp5/main.cc
+++ b/ecp5/main.cc
@@ -149,8 +149,8 @@ std::unique_ptr<Context> ECP5CommandHandler::createContext()
             chipArgs.speed = ArchArgs::SPEED_6;
         }
     }
-
-    return std::unique_ptr<Context>(new Context(chipArgs));
+    auto ctx = std::unique_ptr<Context>(new Context(chipArgs));
+    return ctx;
 }
 
 void ECP5CommandHandler::customAfterLoad(Context *ctx)
diff --git a/generic/arch.cc b/generic/arch.cc
index 77417d27..aca81559 100644
--- a/generic/arch.cc
+++ b/generic/arch.cc
@@ -21,6 +21,7 @@
 #include "nextpnr.h"
 #include "placer1.h"
 #include "router1.h"
+#include "util.h"
 
 NEXTPNR_NAMESPACE_BEGIN
 
@@ -439,7 +440,16 @@ bool Arch::getBudgetOverride(const NetInfo *net_info, const PortRef &sink, delay
 
 // ---------------------------------------------------------------
 
-bool Arch::place() { return placer1(getCtx(), Placer1Cfg(getCtx())); }
+bool Arch::place()
+{
+    std::string placer = str_or_default(settings, id("placer"), defaultPlacer);
+    // FIXME: No HeAP because it needs a list of IO buffers
+    if (placer == "sa") {
+        return placer1(getCtx(), Placer1Cfg(getCtx()));
+    } else {
+        log_error("Generic architecture does not support placer '%s'\n", placer.c_str());
+    }
+}
 
 bool Arch::route() { return router1(getCtx(), Router1Cfg(getCtx())); }
 
@@ -476,4 +486,7 @@ TimingClockingInfo Arch::getPortClockingInfo(const CellInfo *cell, IdString port
 bool Arch::isValidBelForCell(CellInfo *cell, BelId bel) const { return true; }
 bool Arch::isBelLocationValid(BelId bel) const { return true; }
 
+const std::string Arch::defaultPlacer = "sa";
+const std::vector<std::string> Arch::availablePlacers = {"sa"};
+
 NEXTPNR_NAMESPACE_END
diff --git a/generic/arch.h b/generic/arch.h
index dc4258cc..5b5d8c55 100644
--- a/generic/arch.h
+++ b/generic/arch.h
@@ -240,6 +240,9 @@ struct Arch : BaseCtx
 
     bool isValidBelForCell(CellInfo *cell, BelId bel) const;
     bool isBelLocationValid(BelId bel) const;
+
+    static const std::string defaultPlacer;
+    static const std::vector<std::string> availablePlacers;
 };
 
 NEXTPNR_NAMESPACE_END
diff --git a/gui/fpgaviewwidget.cc b/gui/fpgaviewwidget.cc
index c932c3e7..5eab20ed 100644
--- a/gui/fpgaviewwidget.cc
+++ b/gui/fpgaviewwidget.cc
@@ -645,10 +645,10 @@ void FPGAViewWidget::mousePressEvent(QMouseEvent *event)
         return;
 
     bool shift = QApplication::keyboardModifiers().testFlag(Qt::ShiftModifier);
-    bool ctrl  = QApplication::keyboardModifiers().testFlag(Qt::ControlModifier);
+    bool ctrl = QApplication::keyboardModifiers().testFlag(Qt::ControlModifier);
     bool btn_right = event->buttons() & Qt::RightButton;
-    bool btn_mid   = event->buttons() & Qt::MidButton;
-    bool btn_left  = event->buttons() & Qt::LeftButton;
+    bool btn_mid = event->buttons() & Qt::MidButton;
+    bool btn_left = event->buttons() & Qt::LeftButton;
 
     if (btn_right || btn_mid || (btn_left && shift)) {
         lastDragPos_ = event->pos();
@@ -687,8 +687,8 @@ void FPGAViewWidget::mouseMoveEvent(QMouseEvent *event)
 
     bool shift = QApplication::keyboardModifiers().testFlag(Qt::ShiftModifier);
     bool btn_right = event->buttons() & Qt::RightButton;
-    bool btn_mid   = event->buttons() & Qt::MidButton;
-    bool btn_left  = event->buttons() & Qt::LeftButton;
+    bool btn_mid = event->buttons() & Qt::MidButton;
+    bool btn_left = event->buttons() & Qt::LeftButton;
 
     if (btn_right || btn_mid || (btn_left && shift)) {
         const int dx = event->x() - lastDragPos_.x();
diff --git a/ice40/arch.cc b/ice40/arch.cc
index fbe882fc..b0839fa5 100644
--- a/ice40/arch.cc
+++ b/ice40/arch.cc
@@ -26,10 +26,10 @@
 #include "log.h"
 #include "nextpnr.h"
 #include "placer1.h"
+#include "placer_heap.h"
 #include "router1.h"
 #include "timing_opt.h"
 #include "util.h"
-
 NEXTPNR_NAMESPACE_BEGIN
 
 // -----------------------------------------------------------------------
@@ -671,8 +671,18 @@ bool Arch::getBudgetOverride(const NetInfo *net_info, const PortRef &sink, delay
 
 bool Arch::place()
 {
-    if (!placer1(getCtx(), Placer1Cfg(getCtx())))
-        return false;
+    std::string placer = str_or_default(settings, id("placer"), defaultPlacer);
+    if (placer == "heap") {
+        PlacerHeapCfg cfg(getCtx());
+        cfg.ioBufTypes.insert(id_SB_IO);
+        if (!placer_heap(getCtx(), cfg))
+            return false;
+    } else if (placer == "sa") {
+        if (!placer1(getCtx(), Placer1Cfg(getCtx())))
+            return false;
+    } else {
+        log_error("iCE40 architecture does not support placer '%s'\n", placer.c_str());
+    }
     if (bool_or_default(settings, id("opt_timing"), false)) {
         TimingOptCfg tocfg(getCtx());
         tocfg.cellTypes.insert(id_ICESTORM_LC);
@@ -1198,4 +1208,12 @@ void Arch::assignCellInfo(CellInfo *cell)
     }
 }
 
+const std::string Arch::defaultPlacer = "sa";
+
+const std::vector<std::string> Arch::availablePlacers = {"sa",
+#ifdef WITH_HEAP
+                                                         "heap"
+#endif
+};
+
 NEXTPNR_NAMESPACE_END
diff --git a/ice40/arch.h b/ice40/arch.h
index 706043b2..ea29f4f1 100644
--- a/ice40/arch.h
+++ b/ice40/arch.h
@@ -897,6 +897,9 @@ struct Arch : BaseCtx
         IdString glb_net = getWireName(getBelPinWire(bel, id_GLOBAL_BUFFER_OUTPUT));
         return std::stoi(std::string("") + glb_net.str(this).back());
     }
+
+    static const std::string defaultPlacer;
+    static const std::vector<std::string> availablePlacers;
 };
 
 void ice40DelayFuzzerMain(Context *ctx);
diff --git a/ice40/arch_pybindings.cc b/ice40/arch_pybindings.cc
index f0ca584b..bc0bfb84 100644
--- a/ice40/arch_pybindings.cc
+++ b/ice40/arch_pybindings.cc
@@ -144,6 +144,13 @@ void arch_wrap_python()
 
     fn_wrapper_2a_v<Context, decltype(&Context::addClock), &Context::addClock, conv_from_str<IdString>,
                     pass_through<float>>::def_wrap(ctx_cls, "addClock");
+    fn_wrapper_5a_v<Context, decltype(&Context::createRectangularRegion), &Context::createRectangularRegion,
+                    conv_from_str<IdString>, pass_through<int>, pass_through<int>, pass_through<int>,
+                    pass_through<int>>::def_wrap(ctx_cls, "createRectangularRegion");
+    fn_wrapper_2a_v<Context, decltype(&Context::addBelToRegion), &Context::addBelToRegion, conv_from_str<IdString>,
+                    conv_from_str<BelId>>::def_wrap(ctx_cls, "addBelToRegion");
+    fn_wrapper_2a_v<Context, decltype(&Context::constrainCellToRegion), &Context::constrainCellToRegion,
+                    conv_from_str<IdString>, conv_from_str<IdString>>::def_wrap(ctx_cls, "constrainCellToRegion");
 
     WRAP_RANGE(Bel, conv_to_str<BelId>);
     WRAP_RANGE(Wire, conv_to_str<WireId>);
diff --git a/ice40/blinky.pcf b/ice40/examples/blinky/blinky.pcf
index 141dfcc8..141dfcc8 100644
--- a/ice40/blinky.pcf
+++ b/ice40/examples/blinky/blinky.pcf
diff --git a/ice40/blinky.proj b/ice40/examples/blinky/blinky.proj
index f5bb9f88..f5bb9f88 100644
--- a/ice40/blinky.proj
+++ b/ice40/examples/blinky/blinky.proj
diff --git a/ice40/blinky.sh b/ice40/examples/blinky/blinky.sh
index a2326fc3..a2326fc3 100755
--- a/ice40/blinky.sh
+++ b/ice40/examples/blinky/blinky.sh
diff --git a/ice40/blinky.v b/ice40/examples/blinky/blinky.v
index 36eaee86..36eaee86 100644
--- a/ice40/blinky.v
+++ b/ice40/examples/blinky/blinky.v
diff --git a/ice40/blinky.ys b/ice40/examples/blinky/blinky.ys
index a5dd2c85..a5dd2c85 100644
--- a/ice40/blinky.ys
+++ b/ice40/examples/blinky/blinky.ys
diff --git a/ice40/blinky_tb.v b/ice40/examples/blinky/blinky_tb.v
index f80b5e64..f80b5e64 100644
--- a/ice40/blinky_tb.v
+++ b/ice40/examples/blinky/blinky_tb.v
diff --git a/ice40/examples/floorplan/.gitignore b/ice40/examples/floorplan/.gitignore
new file mode 100644
index 00000000..d93659be
--- /dev/null
+++ b/ice40/examples/floorplan/.gitignore
@@ -0,0 +1,4 @@
+*.json
+*.asc
+*.bin
+__pycache__
\ No newline at end of file
diff --git a/ice40/examples/floorplan/floorplan.py b/ice40/examples/floorplan/floorplan.py
new file mode 100644
index 00000000..85c53ccd
--- /dev/null
+++ b/ice40/examples/floorplan/floorplan.py
@@ -0,0 +1,5 @@
+ctx.createRectangularRegion("osc", 1, 1, 1, 4)
+for cell, cellinfo in ctx.cells:
+    if "ringosc" in cellinfo.attrs:
+        print("Floorplanned cell %s" % cell)
+        ctx.constrainCellToRegion(cell, "osc")
diff --git a/ice40/examples/floorplan/floorplan.sh b/ice40/examples/floorplan/floorplan.sh
new file mode 100755
index 00000000..e0ed7a64
--- /dev/null
+++ b/ice40/examples/floorplan/floorplan.sh
@@ -0,0 +1,6 @@
+#!/usr/bin/env bash
+set -ex
+yosys -p "synth_ice40 -top top -json floorplan.json" floorplan.v
+../../../nextpnr-ice40 --up5k --json floorplan.json --pcf icebreaker.pcf --asc floorplan.asc --ignore-loops --pre-place floorplan.py
+icepack floorplan.asc floorplan.bin
+iceprog floorplan.bin
diff --git a/ice40/examples/floorplan/floorplan.v b/ice40/examples/floorplan/floorplan.v
new file mode 100644
index 00000000..8f99ed4e
--- /dev/null
+++ b/ice40/examples/floorplan/floorplan.v
@@ -0,0 +1,22 @@
+module top(output LED1, LED2, LED3, LED4, LED5);
+    localparam N = 31;
+    wire [N:0] x;
+    assign x[0] = x[N];
+
+    genvar ii;
+    generate
+
+    for (ii = 0; ii < N; ii = ii + 1) begin
+        (* ringosc *)
+        SB_LUT4 #(.LUT_INIT(1)) lut_i(.I0(x[ii]), .I1(), .I2(), .I3(), .O(x[ii+1]));
+    end
+    endgenerate
+
+    assign clk = x[N];
+
+
+    reg [19:0] ctr;
+    always @(posedge clk)
+        ctr <= ctr + 1'b1;
+    assign {LED5, LED4, LED3, LED2, LED1} = ctr[19:15];
+endmodule
diff --git a/ice40/examples/floorplan/icebreaker.pcf b/ice40/examples/floorplan/icebreaker.pcf
new file mode 100644
index 00000000..ac7ebf9e
--- /dev/null
+++ b/ice40/examples/floorplan/icebreaker.pcf
@@ -0,0 +1,5 @@
+set_io -nowarn LED1       26
+set_io -nowarn LED2       27
+set_io -nowarn LED3       25
+set_io -nowarn LED4       23
+set_io -nowarn LED5       21
diff --git a/ice40/main.cc b/ice40/main.cc
index 2313c2ae..9b79a08c 100644
--- a/ice40/main.cc
+++ b/ice40/main.cc
@@ -176,7 +176,6 @@ std::unique_ptr<Context> Ice40CommandHandler::createContext()
         ctx->settings[ctx->id("opt_timing")] = "1";
     if (vm.count("pcf-allow-unconstrained"))
         ctx->settings[ctx->id("pcf_allow_unconstrained")] = "1";
-
     return ctx;
 }
 
diff --git a/tests b/tests
-Subproject f29dcbe187b517d01964b1074eb7ff0b90849ee
+Subproject 32a683071758ee59d47e2c5cb29c87882993fac