From abe66efc3c484bf0a9af39e946816d9bcc7f4fdd Mon Sep 17 00:00:00 2001 From: Andrew Eliseev <55635871+jointpoints@users.noreply.github.com> Date: Tue, 12 Mar 2024 17:25:06 +0100 Subject: [PATCH] Version 1.1.0 --- .gitignore | 3 + Cargo.lock | 447 +++++++++++++++++++ Cargo.toml | 16 + LICENSE.md | 9 + README.md | 157 +++++++ build.b64 | 3 + build.py | 166 +++++++ src/main.rs | 236 ++++++++++ src/solver/base_solver.rs | 16 + src/solver/benchmark.rs | 211 +++++++++ src/solver/cplex_solver.rs | 570 ++++++++++++++++++++++++ src/solver/errors.rs | 63 +++ src/solver/mod.rs | 5 + src/solver/tree_decomposition_solver.rs | 404 +++++++++++++++++ src/switch_selection_instance.rs | 226 ++++++++++ src/tree_decomposition.rs | 102 +++++ 16 files changed, 2634 insertions(+) create mode 100644 .gitignore create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 LICENSE.md create mode 100644 README.md create mode 100644 build.b64 create mode 100644 build.py create mode 100644 src/main.rs create mode 100644 src/solver/base_solver.rs create mode 100644 src/solver/benchmark.rs create mode 100644 src/solver/cplex_solver.rs create mode 100644 src/solver/errors.rs create mode 100644 src/solver/mod.rs create mode 100644 src/solver/tree_decomposition_solver.rs create mode 100644 src/switch_selection_instance.rs create mode 100644 src/tree_decomposition.rs diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a83006e --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +/target +/Switch selection +/.vscode diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..05de89d --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,447 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "aho-corasick" +version = "1.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b2969dcb958b36655471fc61f7e416fa76033bdd4bfed0678d8fee1e2d07a1f0" +dependencies = [ + "memchr", +] + +[[package]] +name = "arboretum-td" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3a3ee3383163f08ebb78652f1705349c87b5226dcd1d9b0c9f8adbf2ccc82139" +dependencies = [ + "bitvec", + "fxhash", + "num", + "rand", +] + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "bitvec" +version = "0.19.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "55f93d0ef3363c364d5976646a38f04cf67cfe1d4c8d160cdea02cab2c116b33" +dependencies = [ + "funty", + "radium", + "tap", + "wyz", +] + +[[package]] +name = "byteorder" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fd0f2584146f6f2ef48085050886acf353beff7305ebd1ae69500e27c67f64b" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "const-cstr" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed3d0b5ff30645a68f35ece8cea4556ca14ef8a1651455f789a099a0513532a6" + +[[package]] +name = "cplex_dynamic" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f90c0efc6e0e3205f39858b0c48ee6c8dac5885985053d3ef743f9787ac4b79b" +dependencies = [ + "const-cstr", + "dlopen", + "dlopen_derive", + "glob", + "lazy_static", + "libc", +] + +[[package]] +name = "crabnets" +version = "0.1.0" +source = "git+https://github.com/jointpoints/CrabNets.git?rev=7c316c0#7c316c0f1ca6cd67665ac85e0df01cc501c7c50f" +dependencies = [ + "dyn-clone", + "itertools 0.12.1", + "regex", +] + +[[package]] +name = "dlopen" +version = "0.1.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "71e80ad39f814a9abe68583cd50a2d45c8a67561c3361ab8da240587dda80937" +dependencies = [ + "dlopen_derive", + "lazy_static", + "libc", + "winapi", +] + +[[package]] +name = "dlopen_derive" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f236d9e1b1fbd81cea0f9cbdc8dcc7e8ebcd80e6659cd7cb2ad5f6c05946c581" +dependencies = [ + "libc", + "quote", + "syn", +] + +[[package]] +name = "dyn-clone" +version = "1.0.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0d6ef0072f8a535281e4876be788938b528e9a1d43900b82c2569af7da799125" + +[[package]] +name = "either" +version = "1.10.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "11157ac094ffbdde99aa67b23417ebdd801842852b500e395a45a9c0aac03e4a" + +[[package]] +name = "funty" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fed34cd105917e91daa4da6b3728c47b068749d6a62c59811f06ed2ac71d9da7" + +[[package]] +name = "fxhash" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c31b6d751ae2c7f11320402d34e41349dd1016f8d5d45e48c4312bc8625af50c" +dependencies = [ + "byteorder", +] + +[[package]] +name = "getrandom" +version = "0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "190092ea657667030ac6a35e305e62fc4dd69fd98ac98631e5d3a2b1575a12b5" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + +[[package]] +name = "glob" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2fabcfbdc87f4758337ca535fb41a6d701b65693ce38287d856d1674551ec9b" + +[[package]] +name = "hermit-abi" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d231dfb89cfffdbc30e7fc41579ed6066ad03abda9e567ccafae602b97ec5024" + +[[package]] +name = "itertools" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b1c173a5686ce8bfa551b3563d0c2170bf24ca44da99c7ca4bfdab5418c3fe57" +dependencies = [ + "either", +] + +[[package]] +name = "itertools" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba291022dbbd398a455acf126c1e341954079855bc60dfdda641363bd6922569" +dependencies = [ + "either", +] + +[[package]] +name = "lazy_static" +version = "1.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" + +[[package]] +name = "libc" +version = "0.2.153" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9c198f91728a82281a64e1f4f9eeb25d82cb32a5de251c6bd1b5154d63a8e7bd" + +[[package]] +name = "memchr" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "523dc4f511e55ab87b694dc30d0f820d60906ef06413f93d4d7a1385599cc149" + +[[package]] +name = "num" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b05180d69e3da0e530ba2a1dae5110317e49e3b7f3d41be227dc5f92e49ee7af" +dependencies = [ + "num-bigint", + "num-complex", + "num-integer", + "num-iter", + "num-rational", + "num-traits", +] + +[[package]] +name = "num-bigint" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "608e7659b5c3d7cba262d894801b9ec9d00de989e8a82bd4bef91d08da45cdc0" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-complex" +version = "0.4.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23c6602fda94a57c990fe0df199a035d83576b496aa29f4e634a8ac6004e68a6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.46" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-iter" +version = "0.1.44" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d869c01cc0c455284163fd0092f1f93835385ccab5a98a0dcc497b2f8bf055a9" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +dependencies = [ + "autocfg", + "num-bigint", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "da0df0e5185db44f69b44f26786fe401b6c293d1907744beaa7fa62b2e5a517a" +dependencies = [ + "autocfg", +] + +[[package]] +name = "num_cpus" +version = "1.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4161fcb6d602d4d2081af7c3a45852d875a03dd337a6bfdd6e06407b61342a43" +dependencies = [ + "hermit-abi", + "libc", +] + +[[package]] +name = "ppv-lite86" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de" + +[[package]] +name = "proc-macro2" +version = "0.4.30" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cf3d2011ab5c909338f7887f4fc896d35932e29146c12c8d01da6b22a80ba759" +dependencies = [ + "unicode-xid", +] + +[[package]] +name = "quote" +version = "0.6.13" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ce23b6b870e8f94f81fb0a363d65d86675884b34a09043c81e5562f11c1f8e1" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "radium" +version = "0.5.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "941ba9d78d8e2f7ce474c015eea4d9c6d25b6a3327f9832ee29a4de27f91bbb8" + +[[package]] +name = "rand" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404" +dependencies = [ + "libc", + "rand_chacha", + "rand_core", +] + +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" +dependencies = [ + "getrandom", +] + +[[package]] +name = "rand_xoshiro" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6f97cdb2a36ed4183de61b2f824cc45c9f1037f28afe0a322e9fff4c108b5aaa" +dependencies = [ + "rand_core", +] + +[[package]] +name = "regex" +version = "1.10.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b62dbe01f0b06f9d8dc7d49e05a0785f153b00b2c227856282f671e0318c9b15" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "86b83b8b9847f9bf95ef68afb0b8e6cdb80f498442f5179a29fad448fcc1eaea" +dependencies = [ + "aho-corasick", + "memchr", + "regex-syntax", +] + +[[package]] +name = "regex-syntax" +version = "0.8.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c08c74e62047bb2de4ff487b251e4a92e24f48745648451635cec7d591162d9f" + +[[package]] +name = "switch-selection" +version = "1.1.0" +dependencies = [ + "arboretum-td", + "cplex_dynamic", + "crabnets", + "dyn-clone", + "itertools 0.11.0", + "num_cpus", + "rand", + "rand_xoshiro", +] + +[[package]] +name = "syn" +version = "0.15.44" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ca4b3b69a77cbe1ffc9e198781b7acb0c7365a883670e8f1c1bc66fba79a5c5" +dependencies = [ + "proc-macro2", + "quote", + "unicode-xid", +] + +[[package]] +name = "tap" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "55937e1799185b12863d447f42597ed69d9928686b8d88a1df17376a097d8369" + +[[package]] +name = "unicode-xid" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc72304796d0818e357ead4e000d19c9c174ab23dc11093ac919054d20a6a7fc" + +[[package]] +name = "wasi" +version = "0.11.0+wasi-snapshot-preview1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423" + +[[package]] +name = "winapi" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419" +dependencies = [ + "winapi-i686-pc-windows-gnu", + "winapi-x86_64-pc-windows-gnu", +] + +[[package]] +name = "winapi-i686-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" + +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" + +[[package]] +name = "wyz" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85e60b0d1b5f99db2556934e21937020776a5d31520bf169e851ac44e6420214" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..d0cfe93 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,16 @@ +[package] +name = "switch-selection" +version = "1.1.0" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +arboretum-td = "0.1.0" +cplex_dynamic = "0.1.1" +crabnets = { git = "https://github.com/jointpoints/CrabNets.git", rev = "7c316c0" } +dyn-clone = "1.0.12" +itertools = "0.11.0" +num_cpus = "1.16.0" +rand = "0.8.5" +rand_xoshiro = "0.6.0" diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..4aea7fe --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,9 @@ +# License + +Copyright © Andrew Eliseev (JointPoints), 2023 + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the 'Software') to use, copy, modify, merge, distribute copies of the Software under the same conditions, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..03ce79e --- /dev/null +++ b/README.md @@ -0,0 +1,157 @@ +# Switch selection + + + + + +
+ Version + + Licence + + PSCC paper +
+ + + + + +## Contents + +* [Brief](#brief) +* [Input format](#input-format) +* [Output format](#output-format) +* [What is GNBS?](#what-is-gnbs) +* [Available solvers](#available-solvers) +* [Compilation](#compilation) +* [Usage](#usage) +* [Benchmarking](#benchmarking) + * [Reproduction of the results](#reproduction-of-the-results) + * [Interpretation of the output](#interpretation-of-the-output) + + + +## Brief + +A tool in this repository solves the following problem. + +**Input** + +> * $G$ — a MV distribution grid that has a DG-kernel $D$ such that: +> * $\text{V}(D)$ are all primary substations of $G$. +> * There is an edge replacement sequence that constructs $G$ from $D$ with all replaced edges coming from $\text{E}(D)$. +> * All radial subnetworks of $G$ are degenerate. +> * $p : \text{V}(G) \setminus \text{V}(D) \rightarrow \mathbb{Q}$ — active power at each secondary substation of $G$. +> * $q : \text{V}(G) \setminus \text{V}(D) \rightarrow \mathbb{Q}$ — reactive power at each secondary substation of $G$. +> * $r : \text{E}(G) \rightarrow \mathbb{Q}$ — resistance of each edge. +> * $x : \text{E}(G) \rightarrow \mathbb{Q}$ — reactance of each edge. + +**Output** + +> * $v : \text{V}(D) \rightarrow \{ -10, \dots, 10 \}$ — an optimal tap position for each primary substation. +> * $S \subseteq \text{E}(G)$ — set of edges where switches should be opened. + +For further detail and definitions see our paper (link coming soon). + + + +## Input format + +The program expects a GNBS file as an input. The GNBS file must describe grid $G$ with following attributes defined: + +| Vertex attribute name | Type | Meaning | Possible values | +|:-----------------------:|:----:|:--------|:----------------| +| `is primary substation` | `B` | Flag of primary substations | `T` for primary substations, `F` for secondary substations | +| `p` | `F8` | Active power | A 64-bit float if `is primary substation == F`, `X` otherwise | +| `q` | `F8` | Reactive power | A 64-bit float if `is primary substation == F`, `X` otherwise | + +| Edge attribute name | Type | Meaning | Possible values | +|:-------------------:|:----:|:--------|:----------------| +| `r` | `F8` | Resistance | A 64-bit float | +| `x` | `F8` | Reactance | A 64-bit float | + + + +## Output format + +The program produces a GNBS file as an output. This GNBS file is a full copy of the input file with the following additional attributes defined: + +| Vertex attribute name | Type | Meaning | Possible values | +|:---------------------:|:----:|:--------|:----------------| +| `tap position` | `I1` | Tap position $v$ that defines the base voltage $1 + 0.01 v$ | An integer from $\{ -10, \dots, 10 \}$ if `is primary substation == T`, `X` otherwise | + +| Edge attribute name | Type | Meaning | Possible values | +|:-------------------:|:----:|:--------|:----------------| +| `opened switch` | `B` | Flag of an opened switch | `T` or `F` | + + + +## What is GNBS? + +You can find a full specification of GNBS format [here](https://github.com/jointpoints/GNBSFormat/blob/main/Specification.md). + + + +## Available solvers + +The following solvers are implemented in this tool: + +* `TreeDecompositionSolver` — a solver that solves the problem with dynamic programming using tree decompositions. +* `CPLEXSolver` — a solver that solves the problem formulated as a MILP with the help of CPLEX. + +To use the `CPLEXSolver` or to run the benchmark, you must have a copy of [CPLEX](https://www.ibm.com/products/ilog-cplex-optimization-studio/cplex-optimizer) installed on your computer. CPLEX is proprietary software owned by IBM. If you don't own a licence of CPLEX, you can still use our `TreeDecompositionSolver` without any problems or restrictions. + + + +## Compilation + +Clone this repository and run `build.py`, which is located in the root folder. Note that in order for compilation to be successful, the following must be present on your system: + +* `python` — a Python interpreter, version 3, distributed either with pip or conda. +* `rustc` — a Rust compiler, version 1.76 or newer. +* `cargo` — a Rust build system and package manager, version 1.76 or newer. +* Internet connection. + +The compiled program will be saved in the `Switch selection` folder, which is created automatically. + + + +## Usage + +The tool is supplied with a CLI. Run the following command to learn how to use it: + +``` +.\switch-selection.exe -h +``` +on Windows or +``` +switch-selection -h +``` +on Linux. + + + +## Benchmarking + +#### Reproduction of the results + +To reproduce the results from our PSCC paper, run + +``` +.\switch-selection.exe -b +``` +on Windows or +``` +switch-selection -b +``` +on Linux. You must have a copy of CPLEX installed on your computer to run the benchmark. + +#### Interpretation of the output + +When launched in the benchmark mode, the program will print messages to the console. For each sample, the output will contain a line that looks something like `PDXPDXPDXPD` or, more generally, in the language of regular expressions, `(PDX)*PD`. Here is what these symbols mean: + +* `P` — primary substations have been successfully generated. +* `D` — the corresponding distribution grid has been successfully generated. +* `X` — the generated sample turned out infeasible. + +If the sample is found to be infeasible, the whole process starts again and is repeated until a feasible sample is produced. The time taken to identify infeasibility is not recorded and doesn't affect the metrics. diff --git a/build.b64 b/build.b64 new file mode 100644 index 0000000..14ff937 --- /dev/null +++ b/build.b64 @@ -0,0 +1,3 @@ +IyBFeGFtcGxlIDEKIwojIFRoaXMgaXMgYSBzaW1wbGUgZGlzdHJpYnV0aW9uIGdyaWQgd2l0aCAyIHByaW1hcnkgc3Vic3RhdGlvbnMgdGhhdCBjYW4gYmUgdXNlZCBhcyBhbgojIGlucHV0IGluc3RhbmNlIGZvciB0aGUgc3dpdGNoLXNlbGVjdGlvbiB0b29sLiBUaGlzIGluc3RhbmNlIGNhbiBhbHNvIGJlIGZvdW5kIGluCiMgRmlnLiA3IG9mIG91ciBQU0NDIHBhcGVyLgojCiMgVGhpcyBncmlkIGhhcyB0aGUgZm9sbG93aW5nIHN0cnVjdHVyZToKIwojICAgICAgIDLigJTigJTigJTigJTigJTigJQzCiMgICAgICAvICAgICAgICBcCiMgIFsgMCBdICAgICAgICBbIDEgXQojICAgICAgXCAgICAgICAgLwojICAgICAgIDTigJTigJTigJTigJTigJTigJQ1CiMKIyBOdW1iZXJzIGluIHNxdWFyZSBiYWNrZXRzIHJlcHJlc2VudCBwcmltYXJ5IHN1YnN0YXRpb25zCgoKCiMgVkVSVEVYIEFUVFJJQlVURVMKCkFWIEIgCWlzIHByaW1hcnkgc3Vic3RhdGlvbgpBViBGOAlwCkFWIEY4CXEKCgoKIyBFREdFIEFUVFJJQlVURVMKCkFFIEY4CXIKQUUgRjgJeAoKCgojIFZFUlRJQ0VTCgojIElEICAgICAgICBpcyBwcmltYXJ5IHN1YnN0YXRpb24gICAgICAgIHAgICAgICAgIHEKClYgMCAgICAgICAgIFQgICAgICAgICAgICAgICAgICAgICAgICAgICAgWCAgICAgICAgWApWIDEgICAgICAgICBUICAgICAgICAgICAgICAgICAgICAgICAgICAgIFggICAgICAgIFgKViAyICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAwLjI4ICAgICAwLjAKViAzICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAwLjA1ICAgICAwLjAKViA0ICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMC4xNyAgICAwLjAKViA1ICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMC4wNSAgICAwLjAKCgoKIyBFREdFUwoKIyBJRDEgICAgICAgIElEMiAgICAgICAgciAgICAgICAgeAoKRSAwICAgICAgICAgIDIgICAgICAgICAgMS4wICAgICAgMC4wCkUgMiAgICAgICAgICAzICAgICAgICAgIDEuMCAgICAgIDAuMApFIDMgICAgICAgICAgMSAgICAgICAgICAxLjAgICAgICAwLjAKRSAwICAgICAgICAgIDQgICAgICAgICAgMS4wICAgICAgMC4wCkUgNCAgICAgICAgICA1ICAgICAgICAgIDEuMCAgICAgIDAuMApFIDUgICAgICAgICAgMSAgICAgICAgICAxLjAgICAgICAwLjAK +IyBFeGFtcGxlIDIKIwojIFRoaXMgaXMgYSBzaW1wbGUgZGlzdHJpYnV0aW9uIGdyaWQgd2l0aCAzIHByaW1hcnkgc3Vic3RhdGlvbnMgdGhhdCBjYW4gYmUgdXNlZCBhcyBhbgojIGlucHV0IGluc3RhbmNlIGZvciB0aGUgc3dpdGNoLXNlbGVjdGlvbiB0b29sLgojCiMgVGhpcyBncmlkIGhhcyB0aGUgZm9sbG93aW5nIHN0cnVjdHVyZToKIwojICBbIDAgXeKAlOKAlOKAlAojICAgIOKUgiAgICAgXAojICAgIDUgICAgICAzCiMgICAg4pSCICAgICAgIFwKIyAgICA2ICAgICAgICA0CiMgICAg4pSCICAgICAgICAgXAojICBbIDEgXeKAlOKAlDfigJTigJQ44oCU4oCUWyAyIF0KIwojIE51bWJlcnMgaW4gc3F1YXJlIGJhY2tldHMgcmVwcmVzZW50IHByaW1hcnkgc3Vic3RhdGlvbnMKCgoKIyBWRVJURVggQVRUUklCVVRFUwoKQVYgQglpcyBwcmltYXJ5IHN1YnN0YXRpb24KQVYgRjgJcApBViBGOAlxCgoKCiMgRURHRSBBVFRSSUJVVEVTCgpBRSBGOAlyCkFFIEY4CXgKCgoKIyBWRVJUSUNFUwoKIyBJRCAgICAgICAgaXMgcHJpbWFyeSBzdWJzdGF0aW9uICAgICAgICBwICAgICAgICBxCgpWIDAgICAgICAgICBUICAgICAgICAgICAgICAgICAgICAgICAgICAgIFggICAgICAgIFgKViAxICAgICAgICAgVCAgICAgICAgICAgICAgICAgICAgICAgICAgICBYICAgICAgICBYClYgMiAgICAgICAgIFQgICAgICAgICAgICAgICAgICAgICAgICAgICAgWCAgICAgICAgWApWIDMgICAgICAgICBGICAgICAgICAgICAgICAgICAgICAgICAgICAgIDIuMiAgICAgIDAuMApWIDQgICAgICAgICBGICAgICAgICAgICAgICAgICAgICAgICAgICAgIDEuMzYgICAgIDAuMApWIDUgICAgICAgICBGICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0yLjIgICAgIDAuMApWIDYgICAgICAgICBGICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0xLjM2ICAgIDAuMApWIDcgICAgICAgICBGICAgICAgICAgICAgICAgICAgICAgICAgICAgIC0yLjIgICAgIDAuMApWIDggICAgICAgICBGICAgICAgICAgICAgICAgICAgICAgICAgICAgIDIuMiAgICAgIDAuMAoKCgojIEVER0VTCgojIElEMSAgICAgICAgSUQyICAgICAgICByICAgICAgICB4CgpFIDAgICAgICAgICAgMyAgICAgICAgICAwLjA1ICAgICAwLjAxCkUgMyAgICAgICAgICA0ICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSA0ICAgICAgICAgIDIgICAgICAgICAgMC4wNSAgICAgMC4wMQpFIDAgICAgICAgICAgNSAgICAgICAgICAwLjA1ICAgICAwLjAxCkUgNSAgICAgICAgICA2ICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSA2ICAgICAgICAgIDEgICAgICAgICAgMC4wNSAgICAgMC4wMQpFIDEgICAgICAgICAgNyAgICAgICAgICAwLjA1ICAgICAwLjAxCkUgNyAgICAgICAgICA4ICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSA4ICAgICAgICAgIDIgICAgICAgICAgMC4wNSAgICAgMC4wMQo= +IyBFeGFtcGxlIDMKIwojIFRoaXMgaXMgYSBzaW1wbGUgZGlzdHJpYnV0aW9uIGdyaWQgd2l0aCA0IHByaW1hcnkgc3Vic3RhdGlvbnMgdGhhdCBjYW4gYmUgdXNlZCBhcyBhbgojIGlucHV0IGluc3RhbmNlIGZvciB0aGUgc3dpdGNoLXNlbGVjdGlvbiB0b29sLgojCiMgVGhpcyBncmlkIGhhcyB0aGUgZm9sbG93aW5nIHN0cnVjdHVyZToKIwojICAgICAgNOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlOKAlDUKIyAgICAgLyAgICAgICAgICAgICAgICAgIFwKIyAgWyAwIF3igJTigJTigJQ24oCU4oCU4oCU4oCU4oCU4oCU4oCU4oCUN+KAlOKAlOKAlFsgMSBdCiMgICAvIFwgICAgICAgICAgICAgICAgLy8vIFwKIyAgOCAgMTAgICAgICAgICAgICAgIC8vMTYgMTgKIyAg4pSCICAg4pSCICAgMTLigJTigJTigJTigJTigJTigJTigJQxMy8g4pSCICAg4pSCCiMgIOKUgiAgIOKUgiAgLyAgICAgICAgICAvICDilIIgICDilIIKIyAg4pSCICAg4pSCIC8xNOKAlOKAlOKAlOKAlOKAlOKAlOKAlDE1ICAg4pSCICAg4pSCCiMgIDkgIDExLy8gICAgICAgICAgICAgIDE3IDE5CiMgICBcIC8vLyAgICAgICAgICAgICAgICBcIC8KIyAgWyAyIF3igJTigJTigJQyMOKAlOKAlOKAlOKAlOKAlOKAlDIx4oCU4oCU4oCUWyAzIF0KIyAgICAgXCAgICAgICAgICAgICAgICAgIC8KIyAgICAgIDIy4oCU4oCU4oCU4oCU4oCU4oCU4oCU4oCU4oCU4oCU4oCU4oCU4oCU4oCUMjMKIwojIE51bWJlcnMgaW4gc3F1YXJlIGJhY2tldHMgcmVwcmVzZW50IHByaW1hcnkgc3Vic3RhdGlvbnMKCgoKIyBWRVJURVggQVRUUklCVVRFUwoKQVYgQglpcyBwcmltYXJ5IHN1YnN0YXRpb24KQVYgRjgJcApBViBGOAlxCgoKCiMgRURHRSBBVFRSSUJVVEVTCgpBRSBGOAlyCkFFIEY4CXgKCgoKIyBWRVJUSUNFUwoKIyBJRCAgICAgICAgaXMgcHJpbWFyeSBzdWJzdGF0aW9uICAgICAgICBwICAgICAgICBxCgpWIDAgICAgICAgICBUICAgICAgICAgICAgICAgICAgICAgICAgICAgIFggICAgICAgIFgKViAxICAgICAgICAgVCAgICAgICAgICAgICAgICAgICAgICAgICAgICBYICAgICAgICBYClYgMiAgICAgICAgIFQgICAgICAgICAgICAgICAgICAgICAgICAgICAgWCAgICAgICAgWApWIDMgICAgICAgICBUICAgICAgICAgICAgICAgICAgICAgICAgICAgIFggICAgICAgIFgKViA0ICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICA4LjAgICAgICAwLjAKViA1ICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAxMS4wICAgICAwLjAKViA2ICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMS4wICAgICAwLjAKViA3ICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMS4wICAgICAwLjAKViA4ICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMTUuMCAgICAwLjAKViA5ICAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICA3LjAgICAgICAwLjAKViAxMCAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICA3LjAgICAgICAwLjAKViAxMSAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMTUuMCAgICAwLjAKViAxMiAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtNS4wICAgICAwLjAKViAxMyAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMi4wICAgICAwLjAKViAxNCAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMTAuMCAgICAwLjAKViAxNSAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtMS4wICAgICAwLjAKViAxNiAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAtOC4wICAgICAwLjAKViAxNyAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICA4LjAgICAgICAwLjAKViAxOCAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICA1LjAgICAgICAwLjAKViAxOSAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICA2LjAgICAgICAwLjAKViAyMCAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAxLjAgICAgICAwLjAKViAyMSAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAxLjAgICAgICAwLjAKViAyMiAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAxLjAgICAgICAwLjAKViAyMyAgICAgICAgRiAgICAgICAgICAgICAgICAgICAgICAgICAgICAxLjAgICAgICAwLjAKCgoKIyBFREdFUwoKIyBJRDEgICAgICAgIElEMiAgICAgICAgciAgICAgICAgeAoKRSAwICAgICAgICAgIDQgICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDQgICAgICAgICAgNSAgICAgICAgICAwLjAxICAgICAwLjAxCkUgNSAgICAgICAgICAxICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAwICAgICAgICAgIDYgICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDYgICAgICAgICAgNyAgICAgICAgICAwLjAxICAgICAwLjAxCkUgNyAgICAgICAgICAxICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAwICAgICAgICAgIDggICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDggICAgICAgICAgOSAgICAgICAgICAwLjAxICAgICAwLjAxCkUgOSAgICAgICAgICAyICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAwICAgICAgICAgIDEwICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDEwICAgICAgICAgMTEgICAgICAgICAwLjAxICAgICAwLjAxCkUgMTEgICAgICAgICAyICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAyICAgICAgICAgIDEyICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDEyICAgICAgICAgMTMgICAgICAgICAwLjAxICAgICAwLjAxCkUgMTMgICAgICAgICAxICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAyICAgICAgICAgIDE0ICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDE0ICAgICAgICAgMTUgICAgICAgICAwLjAxICAgICAwLjAxCkUgMTUgICAgICAgICAxICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAxICAgICAgICAgIDE2ICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDE2ICAgICAgICAgMTcgICAgICAgICAwLjAxICAgICAwLjAxCkUgMTcgICAgICAgICAzICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAxICAgICAgICAgIDE4ICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDE4ICAgICAgICAgMTkgICAgICAgICAwLjAxICAgICAwLjAxCkUgMTkgICAgICAgICAzICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAyICAgICAgICAgIDIwICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDIwICAgICAgICAgMjEgICAgICAgICAwLjAxICAgICAwLjAxCkUgMjEgICAgICAgICAzICAgICAgICAgIDAuMDEgICAgIDAuMDEKRSAyICAgICAgICAgIDIyICAgICAgICAgMC4wMSAgICAgMC4wMQpFIDIyICAgICAgICAgMjMgICAgICAgICAwLjAxICAgICAwLjAxCkUgMjMgICAgICAgICAzICAgICAgICAgIDAuMDEgICAgIDAuMDEK \ No newline at end of file diff --git a/build.py b/build.py new file mode 100644 index 0000000..1b85f36 --- /dev/null +++ b/build.py @@ -0,0 +1,166 @@ +# A script to manually compile the tool +# Change this value to False to make a debug build +BUILD_RELEASE = True +REMOVE_TARGET = True + + + + + +error = False +warning = False + +import os +import subprocess +import shutil +import re +import pathlib +import base64 +try: + from termcolor import cprint +except ModuleNotFoundError: + print('''It seems like your Python environment doesn't have a termcolor module installed. +This module isn't neccessary but it's recommended for a better experience with this build script. + +Shall it be installed now? +\t- yes, let's install it with (p)ip using 'pip install termcolor'. +\t- yes, let's install it with (c)onda using 'conda install termcolor'. +\t- (n)o, let's continue without it. +\t- (a)bort the execution of this script.''') + command = input('Enter a symbol in brackets to choose from one of the options: ') + if command == 'p': + return_value = subprocess.run(['pip', 'install', 'termcolor'], capture_output=True, shell=True) + if return_value.returncode == 0: + from termcolor import cprint + else: + print('ERROR! Command \'pip install termcolor\' wasn\'t executed successfully. The exectution of this script will be aborted.') + error = True + elif command == 'c': + return_value = subprocess.run(['conda', 'install', 'termcolor'], capture_output=True, shell=True) + if return_value.returncode == 0: + from termcolor import cprint + else: + print('ERROR! Command \'conda install termcolor\' wasn\'t executed successfully. The exectution of this script will be aborted.') + error = True + elif command == 'n': + global cprint + def cprint(message, _colour, end='\n'): + print(message, end=end) + elif command == 'a': + error = True + else: + print(f'ERROR! Unknown command {command}. The execution of this script will be aborted.') + error = True + + + + + +width = os.get_terminal_size().columns +return_value = None + +if not error: + os.system('cls' if os.name == 'nt' else 'clear') + print(' '.join(['*'] * (width // 2 + (width % 2)))) + print('* SWITCH SELECTION BUILD SCRIPT' + ' ' * (width - 32 - (width % 2 == 0)) + '*') + print('*' + ' ' * (width - 2 - (width % 2 == 0)) + '*') + print('*' + ' ' * (width - 34 - (width % 2 == 0)) + 'by Andrew Eliseev (JointPoints) *') + print(' '.join(['*'] * (width // 2 + (width % 2)))) + print() + + if not BUILD_RELEASE: + cprint('You are about to make a DEBUG build! If this is not your intention, open the source code', 'red') + cprint('of this script and set BUILD_RELEASE to True!', 'red') + print() + + # Check the presence of the Rust compiler + print('1. Checking the availability of necessary programs') + print('\trustc\t: ', end='') + return_value = subprocess.run(['rustc'], capture_output=True, shell=True) + if return_value.returncode == 0: + cprint('OK', 'green') + else: + cprint('NOT FOUND', 'red') + error = True + print('\tcargo\t: ', end='') + return_value = subprocess.run(['cargo'], capture_output=True, shell=True) + if return_value.returncode == 0: + cprint('OK', 'green') + else: + cprint('NOT FOUND', 'red') + error = True + print() + if error: + cprint('ERROR! ', 'red', end='') + print('Rust compiler or Cargo build system not found.') + print() + +# Try to build the tool with Cargo +if not error: + print('2. Building the tool (might take some time)') + return_value = subprocess.run(['cargo', 'build', '-r'] if BUILD_RELEASE else ['cargo', 'build'], capture_output=True, shell=True) + if return_value.returncode != 0: + print() + cprint('ERROR! ', 'red', end='') + print('Compilation failed. The raw output of the Rust compiler follows.') + print(return_value.stderr.decode()) + print() + error = True + else: + cprint('\tOK', 'green') + print() + +# Replace the tool to our custom-made folder +if not error: + print('3. Final arrangements') + os.makedirs('./Switch selection/Graphs', exist_ok=True) + try: + subfolder_name = 'release' if BUILD_RELEASE else 'debug' + shutil.copy(f'./target/{subfolder_name}/switch-selection' + '.exe' if os.name == 'nt' else '', './Switch selection') + except shutil.SameFileError: + pass + if BUILD_RELEASE and REMOVE_TARGET: + shutil.rmtree('./target') + with open('./build.b64') as b64file: + # Decode build.b64 into different files + files = (('./Switch selection/Graphs/example1.gnbs', 'w', -1, 'utf-8'), ('./Switch selection/Graphs/example2.gnbs', 'w', -1, 'utf-8'), ('./Switch selection/Graphs/example3.gnbs', 'w', -1, 'utf-8')) + for file in files: + print(f'\t{file[0][19:]:<30}: ', end='') + line = b64file.readline() + with open(*file) as f: + decoded = base64.b64decode(line) + f.write(decoded.decode() if file[1] == 'w' else decoded) + cprint('OK', 'green') + # Try to find a CPLEX dynamic library + cplex_lib_name = 'cplex*.' + 'dll' if os.name == 'nt' else 'so' + print(f'\t{cplex_lib_name:<30}: ', end='') + matching_env_vars = [] + for env_var_name in os.environ.keys(): + if re.match(r'CPLEX_STUDIO_DIR[0-9]+', env_var_name) != None: + matching_env_vars.append(env_var_name) + matching_env_vars = sorted(matching_env_vars)[::-1] + for env_var_name in matching_env_vars: + for cplex_dll_path in pathlib.Path(os.environ[env_var_name]).rglob(cplex_lib_name): + shutil.copy(cplex_dll_path.resolve(), './Switch selection') + cprint('OK', 'green') + break + else: + continue + break + else: + cprint('NOT FOUND', 'red') + warning = True + print() + +# End +if not error: + cprint('SUCCESS! ', 'green', end='') + print('Build finished. You can find your compiled tool in the \'Switch selection\' folder.') + print() + if warning: + cprint('WARNING! ', 'light_yellow', end='') + print('CPLEX dynamic library file (cplex*.dll or cplex*.so) couldn\'t be automatically located on your system. Please, copy it manually into the \'Switch selection\' folder. Otherwise, you won\'t be able to use CPLEXSolver.') + print() + +print('Press to finish...', end='') +input() diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..9ee5521 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,236 @@ +mod switch_selection_instance; +mod tree_decomposition; +mod solver; + +use std::{env, process::exit, time::Instant}; +use crabnets::{io::IO, Graph}; +use switch_selection_instance::{SwitchSelectionInstance, SwitchSelectionGraph}; +use solver::{base_solver::BaseSolver, cplex_solver::CPLEXSolver, tree_decomposition_solver::TreeDecompositionSolver, benchmark::start_benchmark}; +use crate::solver::base_solver::TapValue; + + + + + +macro_rules! pretty_panic { + ($message: expr) => { + { + println!("{}", $message); + exit(1); + } + }; +} + +macro_rules! pretty_unwrap { + ($result: expr) => { + match $result { + Ok(value) => value, + Err(error) => pretty_panic!(error), + } + }; +} + + + + + +fn main() { + const HELP_STRING: &str = +"switch-selection + +A tool to solve a variation of the SwitchSelection problem described in our PSCC paper. See our +repository https://github.com/EINS-TUDa/PSCC2024-SwitchSelection for a detailed user guide. + +USAGE + OPTION 1: switch-selection [...] (-h|--help) [...] + OPTION 2: switch-selection (-b|--benchmark) + OPTION 3: switch-selection [] + + OPTION 1 prints this message. OPTION 2 launches benchmark: use it to reproduce the results from + our paper. Use OPTION 3 to solve a specific instance of the SwitchSelection problem. + +ARGUMENTS + (-i|--input) PATH Set the path to the input file in GNBS format. + Default value if this argument is omitted: -i input.gnbs + (-o|--output) PATH Set the path to the output file in GNBS format. If the file doesn't + exist, it'll be created automatically. If the file exists, its + contents will be rewritten. + Default value if this argument is omitted: -o output.gnbs + (-s|--solver) SOLVER Set a solver to solve the problem instance with. Possible values for + SOLVER: + o TreeDecompositionSolver - solve the problem using the + dynamic programming approach described in our paper. + o CPLEXSolver - solve the problem in its MILP formulation + using CPLEX (requires CPLEX to be installed). + Default value if this argument is omitted: -s TreeDecompositionSolver + +OPTIONS + --dgkernel [PATH] Save a DG-kernel of the input graph into a GNBS file. If PATH is not + given, value 'dgkernel.gnbs' is assumed. + +EXAMPLES + switch-selection + Equivalent to 'switch-selection -i input.gnbs -o output.gnbs -s TreeDecompositionSolver'. + switch-selection -i ./Graphs/example1.gnbs -s CPLEXSolver + Solve the SwitchSelection instance given by ./Graphs/example1.gnbs with CPLEX, save the + optimal solution into output.gnbs. + switch-selection -o 123.gnbs --dgkernel dgk.gnbs + Solve the SwitchSelection instance given by input.gnbs with TreeDecompositionSolver, save + the optimal solution into 123.gnbs and save the DG-kernel into dgk.gnbs."; + + + + let mut benchmark_mode: Option = None; + let mut input_path: String = "input.gnbs".to_string(); + let mut output_path: String = "output.gnbs".to_string(); + let mut solver_name: String = "TreeDecompositionSolver".to_string(); + let mut dg_kernel_path: Option = None; + let mut timeit: Option<(usize, usize)> = None; + let mut target_path: &mut String = &mut input_path; + // Command-line parser states + enum CLParserState { + ExpectParameter, + ExpectPath, + ExpectPathOrParameter, + ExpectSolver, + ExpectNumber1OrParameter, + ExpectNumber2, + } + + + + // Read command-line arguments + let mut state: CLParserState = CLParserState::ExpectParameter; + let mut arguments: env::Args = env::args(); + arguments.next(); + for argument in arguments { + match argument.as_str() { + "-b" | "--benchmark" => if benchmark_mode == Some(false) { + pretty_panic!(format!("You can't use {} together with any other parameters.", argument)); + } else { + benchmark_mode = Some(true); + }, + "-h" | "--help" => { + println!("{}", HELP_STRING); + exit(0); + }, + "-i" | "--input" => match state { + CLParserState::ExpectParameter | CLParserState::ExpectPathOrParameter | CLParserState::ExpectNumber1OrParameter => { + if benchmark_mode == Some(true) { + pretty_panic!(format!("You can't use {} in the benchmark mode.", argument)); + } + benchmark_mode = Some(false); + target_path = &mut input_path; + state = CLParserState::ExpectPath; + }, + _ => pretty_panic!(format!("Unexpected command-line value {}.", argument)), + }, + "-o" | "--output" => match state { + CLParserState::ExpectParameter | CLParserState::ExpectPathOrParameter | CLParserState::ExpectNumber1OrParameter => { + if benchmark_mode == Some(true) { + pretty_panic!(format!("You can't use {} in the benchmark mode.", argument)); + } + benchmark_mode = Some(false); + target_path = &mut output_path; + state = CLParserState::ExpectPath; + }, + _ => pretty_panic!(format!("Unexpected command-line value {}.", argument)), + }, + "-s" | "--solver" => match state { + CLParserState::ExpectParameter | CLParserState::ExpectPathOrParameter | CLParserState::ExpectNumber1OrParameter => { + if benchmark_mode == Some(true) { + pretty_panic!(format!("You can't use {} in the benchmark mode.", argument)); + } + benchmark_mode = Some(false); + state = CLParserState::ExpectSolver; + }, + _ => pretty_panic!(format!("Unexpected command-line value {}.", argument)), + }, + "--dgkernel" => match state { + CLParserState::ExpectParameter | CLParserState::ExpectNumber1OrParameter => { + if benchmark_mode == Some(true) { + pretty_panic!(format!("You can't use {} in the benchmark mode.", argument)); + } + benchmark_mode = Some(false); + dg_kernel_path = Some("dgkernel.gnbs".to_string()); + state = CLParserState::ExpectPathOrParameter; + }, + _ => pretty_panic!(format!("Unexpected command-line value {}.", argument)), + }, + "--timeit" => match state { + CLParserState::ExpectParameter | CLParserState::ExpectPathOrParameter => { + if benchmark_mode == Some(true) { + pretty_panic!(format!("You can't use {} in the benchmark mode.", argument)); + } + benchmark_mode = Some(false); + timeit = Some((100, 10)); + state = CLParserState::ExpectNumber1OrParameter; + }, + _ => pretty_panic!(format!("Unexpected command-line value {}.", argument)), + }, + a => match state { + CLParserState::ExpectParameter => pretty_panic!(format!("Unknown command-line parameter {}.", a)), + CLParserState::ExpectPath => { + *target_path = argument; + state = CLParserState::ExpectParameter; + }, + CLParserState::ExpectPathOrParameter => { + dg_kernel_path = Some(a.to_string()); + state = CLParserState::ExpectParameter; + } + CLParserState::ExpectSolver => { + match a { + "TreeDecompositionSolver" | "CPLEXSolver" => solver_name = a.to_string(), + _ => pretty_panic!(format!("Unknown solver {}.", a)), + } + state = CLParserState::ExpectParameter; + }, + CLParserState::ExpectNumber1OrParameter => { + timeit = Some((pretty_unwrap!(a.parse()), 0)); + state = CLParserState::ExpectNumber2; + }, + CLParserState::ExpectNumber2 => { + timeit = Some((timeit.unwrap().0, pretty_unwrap!(a.parse()))); + state = CLParserState::ExpectParameter; + }, + }, + } + } + match state { + CLParserState::ExpectParameter | CLParserState::ExpectPathOrParameter | CLParserState::ExpectNumber1OrParameter => (), + _ => pretty_panic!(format!("Unexpected end of command line.")), + } + + + + // Act according to the parameters set + if benchmark_mode == Some(true) { + pretty_unwrap!(start_benchmark(2, 10, 30, 20, 1, 0)); + } else { + let input: SwitchSelectionGraph = pretty_unwrap!(Graph::from_file(&input_path)); + let problem_instance: SwitchSelectionInstance = pretty_unwrap!(SwitchSelectionInstance::new(input)); + if let Some(value) = dg_kernel_path { + pretty_unwrap!(problem_instance.dg_kernel_for_switch_selection().into_file(&value)); + } + let solver_begin_time: Instant = Instant::now(); + match solver_name.as_str() { + "TreeDecompositionSolver" => { + let mut solver: TreeDecompositionSolver = pretty_unwrap!(TreeDecompositionSolver::with_input(problem_instance)); + pretty_unwrap!(solver.solve()); + println!("{} solved the problem instance in {} s.", solver_name, solver_begin_time.elapsed().as_secs_f64()); + let solution: (SwitchSelectionGraph, TapValue) = solver.get_solution().unwrap(); + println!("Objective value = {}.", solution.1); + pretty_unwrap!(solution.0.into_file(&output_path)); + }, + "CPLEXSolver" => { + let mut solver: CPLEXSolver = pretty_unwrap!(CPLEXSolver::with_input(problem_instance)); + pretty_unwrap!(solver.solve()); + println!("{} solved the problem instance in {} s.", solver_name, solver_begin_time.elapsed().as_secs_f64()); + let solution: (SwitchSelectionGraph, TapValue) = solver.get_solution().unwrap(); + println!("Objective value = {}.", solution.1); + pretty_unwrap!(solution.0.into_file(&output_path)); + }, + _ => (), + } + } +} diff --git a/src/solver/base_solver.rs b/src/solver/base_solver.rs new file mode 100644 index 0000000..6132d7d --- /dev/null +++ b/src/solver/base_solver.rs @@ -0,0 +1,16 @@ +use crate::switch_selection_instance::{SwitchSelectionInstance, SwitchSelectionGraph}; +use super::errors::SolverError; + + + + + +pub type TapValue = i8; + + + +pub trait BaseSolver: Sized { + fn with_input(input: SwitchSelectionInstance) -> Result; + fn get_solution(&self) -> Option<(SwitchSelectionGraph, TapValue)>; + fn solve(&mut self) -> Result<(), SolverError>; +} diff --git a/src/solver/benchmark.rs b/src/solver/benchmark.rs new file mode 100644 index 0000000..61cf029 --- /dev/null +++ b/src/solver/benchmark.rs @@ -0,0 +1,211 @@ +use std::{cmp::Ordering, io::{stdout, Write}, time::Instant}; +use crabnets::{topology_tests::TopologyTests, BasicImmutableGraph, BasicMutableGraph}; +use itertools::Itertools; +use rand::{Rng, distributions::Uniform, prelude::Distribution, seq::IteratorRandom}; +use rand_xoshiro::{Xoroshiro128PlusPlus, rand_core::SeedableRng}; +use crate::switch_selection_instance::{SwitchSelectionGraph, SwitchSelectionInstance}; +use super::{cplex_solver::CPLEXSolver, base_solver::BaseSolver, errors::SolverError, tree_decomposition_solver::TreeDecompositionSolver}; + + + + + +type PRNG = Xoroshiro128PlusPlus; + + + +fn random_partial_k_tree(treewidth: usize, vertex_count: u8, prng: &mut PRNG) -> SwitchSelectionGraph { + loop { + let mut answer = SwitchSelectionGraph::new(); + // Populate with a (treewidth + 1)-clique + for vertex_id1 in 0..=treewidth { + answer.add_v(None); + for vertex_id2 in 0..vertex_id1 { + answer.add_e(&vertex_id1, &vertex_id2, false, None).unwrap(); + } + } + // Add other vertices one by one + for new_vertex_id in (treewidth + 1)..(vertex_count as usize) { + answer.add_v(None); + // Sample a random clique + let mut clique_candidate; + 'clique_check: + loop { + clique_candidate = (0..new_vertex_id).choose_multiple(prng, treewidth); + for vertex_i1 in 0..treewidth { + for vertex_i2 in 0..vertex_i1 { + if answer.contains_e(&clique_candidate[vertex_i1], &clique_candidate[vertex_i2], &0).is_none() { + continue 'clique_check; + } + } + } + break; + } + // Attach the new vertex to the found clique + for clique_vertex_id in clique_candidate { + answer.add_e(&new_vertex_id, &clique_vertex_id, false, None).unwrap(); + } + } + // Randomly remove some edges + let edges_to_remove_count: usize = prng.sample(Uniform::new(treewidth, answer.count_e() / 2)); + for _ in 0..edges_to_remove_count { + let edge = answer.iter_e().sorted_by(|x, y| + match x.id1.cmp(&y.id1) { + Ordering::Less => Ordering::Less, + Ordering::Equal => x.id2.cmp(&y.id2), + Ordering::Greater => Ordering::Greater, + } + ).choose_stable(prng).unwrap(); + if answer.v_degree(&edge.id1).unwrap() > 1 && answer.v_degree(&edge.id2).unwrap() > 1 + && (edge.id1 > treewidth || edge.id2 > treewidth) { + answer.remove_e(&edge.id1, &edge.id2, &0).unwrap(); + } + } + if answer.is_connected() { + return answer; + } + } +} + +pub fn start_benchmark( + treewidth: usize, + min_primary_substation_count: u8, + max_primary_substation_count: u8, + sample_count: u8, + sample_repeat: u8, + sample_ignore: u8 +) -> Result<(), SolverError> { + let mut prng = PRNG::seed_from_u64(13374); + // Random distributions + let feeder_count_distribution: Uniform = Uniform::new(2, 6); + let substation_count_distribution: Uniform = Uniform::new(5, 11); + let rx_distribution: Uniform = Uniform::new(0.01, 0.05); + // Benchmarking results + let mut cplex_max_times: Vec = Vec::with_capacity((max_primary_substation_count - min_primary_substation_count + 1) as usize); + let mut cplex_avg_times: Vec = Vec::with_capacity((max_primary_substation_count - min_primary_substation_count + 1) as usize); + let mut cplex_min_times: Vec = Vec::with_capacity((max_primary_substation_count - min_primary_substation_count + 1) as usize); + let mut td_max_times: Vec = Vec::with_capacity((max_primary_substation_count - min_primary_substation_count + 1) as usize); + let mut td_avg_times: Vec = Vec::with_capacity((max_primary_substation_count - min_primary_substation_count + 1) as usize); + let mut td_min_times: Vec = Vec::with_capacity((max_primary_substation_count - min_primary_substation_count + 1) as usize); + // Test random distribution grids + for primary_substation_count in min_primary_substation_count..=max_primary_substation_count { + println!("Treewidth = {}, # primary substations = {}", treewidth, primary_substation_count); + let mut cplex_times: Vec = Vec::with_capacity(sample_count as usize); + let mut td_times: Vec = Vec::with_capacity(sample_count as usize); + let mut successful_samples: u8 = 0; + let mut restart = false; + loop { + if restart { + print!("X"); + stdout().flush().unwrap(); + } else { + print!("\tSample {}\n\t\t", successful_samples + 1); + } + // Generate a random partial k-tree + let mut graph: SwitchSelectionGraph = random_partial_k_tree(treewidth, primary_substation_count, &mut prng); + print!("P"); + stdout().flush().unwrap(); + // Replace each edge of the graph with a bunch of feeders + let old_edges = graph.iter_e().sorted_by(|x, y| + match x.id1.cmp(&y.id1) { + Ordering::Less => Ordering::Less, + Ordering::Equal => x.id2.cmp(&y.id2), + Ordering::Greater => Ordering::Greater, + } + ); + for old_edge in old_edges { + let feeder_count: u8 = feeder_count_distribution.sample(&mut prng); + graph.v_attrs_mut(&old_edge.id1).unwrap().tap_position = Some(0); + graph.v_attrs_mut(&old_edge.id2).unwrap().tap_position = Some(0); + graph.remove_e(&old_edge.id1, &old_edge.id2, &0).unwrap(); + for _ in 0..feeder_count { + let substation_count: u8 = substation_count_distribution.sample(&mut prng); + let mut last_substation_id = old_edge.id1; + let pq_distribution: Uniform = Uniform::new(-0.5, 0.5); + for _ in 0..substation_count { + let new_substation_id: usize = graph.add_v(None); + { + let attributes = graph.v_attrs_mut(&new_substation_id).unwrap(); + attributes.p = pq_distribution.sample(&mut prng); + attributes.q = pq_distribution.sample(&mut prng); + } + graph.add_e(&last_substation_id, &new_substation_id, false, None).unwrap(); + { + let attributes = graph.e_attrs_mut(&last_substation_id, &new_substation_id, &0).unwrap(); + attributes.r = rx_distribution.sample(&mut prng); + attributes.x = rx_distribution.sample(&mut prng); + } + last_substation_id = new_substation_id; + } + graph.add_e(&last_substation_id, &old_edge.id2, false, None).unwrap(); + { + let attributes = graph.e_attrs_mut(&last_substation_id, &old_edge.id2, &0).unwrap(); + attributes.r = rx_distribution.sample(&mut prng); + attributes.x = rx_distribution.sample(&mut prng); + } + } + } + print!("D"); + stdout().flush().unwrap(); + // Create a problem instance instance out of graph + let instance = SwitchSelectionInstance::new(graph).unwrap(); + // Time the TreeDecompositionSolver + let mut solver: TreeDecompositionSolver = TreeDecompositionSolver::with_input(instance.clone())?; + match timeit(&mut solver, sample_repeat, sample_ignore) { + Ok(value) => td_times.push(value), + Err(_) => { + restart = true; + continue; + }, + } + println!("\n\t\tTreeDecompositionSolver finished ({} s). Optimal value = {}.", td_times.last().unwrap(), solver.get_solution().unwrap().1); + // Time the CPLEXSolver + let mut solver: CPLEXSolver = CPLEXSolver::with_input(instance.clone())?; + match timeit(&mut solver, sample_repeat, sample_ignore) { + Ok(value) => cplex_times.push(value), + Err(_) => { + restart = true; + continue; + }, + } + println!("\t\tCPLEXSolver finished ({} s). Optimal value = {}.", cplex_times.last().unwrap(), solver.get_solution().unwrap().1); + successful_samples += 1; + restart = false; + if successful_samples == sample_count { + break; + } + } + // Print the results for this number of primary substations + println!("\tResults for this # primary substations"); + println!("\tTreeDecompositionSolver: {:?}", td_times); + println!("\tCPLEXSolver: {:?}", cplex_times); + // Remember the results for this number of primary substations + cplex_max_times.push(cplex_times.iter().map(|x: &f64| *x).reduce(f64::max).unwrap()); + cplex_avg_times.push(cplex_times.iter().sum::() / sample_count as f64); + cplex_min_times.push(cplex_times.iter().map(|x: &f64| *x).reduce(f64::min).unwrap()); + td_max_times.push(td_times.iter().map(|x: &f64| *x).reduce(f64::max).unwrap()); + td_avg_times.push(td_times.iter().sum::() / sample_count as f64); + td_min_times.push(td_times.iter().map(|x: &f64| *x).reduce(f64::min).unwrap()); + } + println!("Results"); + println!("TreeDecompositionSolver max: {:?}", td_max_times); + println!("TreeDecompositionSolver avg: {:?}", td_avg_times); + println!("TreeDecompositionSolver min: {:?}", td_min_times); + println!("CPLEXSolver max: {:?}", cplex_max_times); + println!("CPLEXSolver avg: {:?}", cplex_avg_times); + println!("CPLEXSolver min: {:?}", cplex_min_times); + Ok(()) +} + +pub fn timeit(solver: &mut S, repeat: u8, ignore: u8) -> Result { + for _ in 0..ignore { + solver.solve()?; + } + let mut cumulative_time: f64 = 0.0; + for _ in 0..(repeat - ignore) { + let solver_begin_time: Instant = Instant::now(); + solver.solve()?; + cumulative_time += solver_begin_time.elapsed().as_secs_f64(); + } + Ok(cumulative_time / repeat as f64) +} diff --git a/src/solver/cplex_solver.rs b/src/solver/cplex_solver.rs new file mode 100644 index 0000000..2671708 --- /dev/null +++ b/src/solver/cplex_solver.rs @@ -0,0 +1,570 @@ +use std::{collections::HashMap, pin::Pin, ptr::NonNull}; +use cplex_dynamic::{Constraint, ConstraintType, Env, Problem, ProblemType, Solution, Variable, VariableType, VariableValue, WeightedVariable}; +use crabnets::{BasicImmutableGraph, BasicMutableGraph, ImmutableGraphContainer}; +use crate::switch_selection_instance::{SwitchSelectionGraph, SwitchSelectionInstance}; +use super::{base_solver::*, errors::SolverError}; + + + + + +macro_rules! cplex_unwrap { + ($expr: expr) => { + match $expr { + Ok(value) => value, + Err(error) => return Err(SolverError::from_string(format!("CPLEXSolver. {}", error))), + } + }; +} + + + +pub struct CPLEXSolverCore<'a> { + input: SwitchSelectionInstance, + solution: Option, + variables: HashMap, + problem: Option>, + env: Env, +} + + + +pub type CPLEXSolver<'a> = Pin>>; + + + +trait CPLEXSolverTools<'a> { + fn get_problem_mut(&mut self) -> &mut Problem<'a>; +} + + + +// CPLEXSolver::CPLEXSolverTools +impl<'a> CPLEXSolverTools<'a> for CPLEXSolver<'a> { + fn get_problem_mut(&mut self) -> &mut Problem<'a> { + unsafe { + self.as_mut().get_unchecked_mut().problem.as_mut().unwrap() + } + } +} + +// CPLEXSolver::BaseSolver +impl<'a> BaseSolver for CPLEXSolver<'a> { + fn with_input(input: SwitchSelectionInstance) -> Result { + let mut solver: CPLEXSolver = Box::pin(CPLEXSolverCore { input: input.clone(), env: cplex_unwrap!(Env::new()), variables: HashMap::new(), problem: None, solution: None }); + unsafe { + // Create a problem instance + solver.as_mut().get_unchecked_mut().problem = Some(cplex_unwrap!(Problem::new(NonNull::from(&solver.as_ref().env).as_ref(), "name"))); + } + // Populate the problem with variables + // * max_tap_abs : {0, ..., 10} -- max |tap(s)| = max tap_abs(s) over all s in P(input) + // * u(s) for s in S(input) : [0.81, 1.21] -- square voltage at substation s + // * tap(i, s) for s in P(input) : {0, 1} -- whether tap at a primary substation s is in position i in {-10, ..., 10} + // * tap_abs(s) for s in P(input) : {0, ..., 10} -- |tap(s)| = 1 * (tap(-1, s) + tap(1, s)) + ... + 10 * (tap(-10, s) + tap(10, s)) + // * part(s) for s in S(input) \ P(input) : {0, 1} -- shows to which primary substation s should be attributed to: left (0) or right (1) + // * u_right(s) for s in S(input) \ P(input) : [0.0, 1.21] -- an alias for part(s) * u(s') where s' is the substation to the right of s + // * u_left(s) for s in S(input) \ P(input) : [0.0, 1.21] -- an alias for (1 - part(s)) * u(s') where s' is the substation to the left of s + let mut variables: HashMap = HashMap::new(); + variables.insert( + "max_tap_abs".to_string(), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Integer, 1.0, 0.0, 10.0, "max_tap_abs"))) + ); + for substation_id in input.iter_v() { + variables.insert( + format!("u({})", substation_id), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Continuous, 0.0, 0.81, 1.21, format!("u({})", substation_id)))) + ); + if input.v_attrs(&substation_id).unwrap().tap_position.is_some() { + for tap_position in -10..=10 { + variables.insert( + format!("tap({},{})", tap_position, substation_id), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Integer, 0.0, 0.0, 1.0, format!("tap({},{})", tap_position, substation_id)))) + ); + } + variables.insert( + format!("tap_abs({})", substation_id), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Integer, 0.0, 0.0, 10.0, format!("tap_abs({})", substation_id)))) + ); + } else { + variables.insert( + format!("part({})", substation_id), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Integer, 0.0, 0.0, 1.0, format!("part({})", substation_id)))) + ); + variables.insert( + format!("u_right({})", substation_id), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Continuous, 0.0, 0.0, 1.21, format!("u_right({})", substation_id)))) + ); + variables.insert( + format!("u_left({})", substation_id), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Continuous, 0.0, 0.0, 1.21, format!("u_left({})", substation_id)))) + ); + } + } + // Traverse each line, add the remaining variables... + // * right_part(s1, s2) for s1, s2 in S(input) \ P(input) : {0, 1} -- an alias for part(s1) * part(s2) + // * left_part(s1, s2) for s1, s2 in S(input) \ P(input) : {0, 1} -- an alias for (1 - part(s1)) * (1 - part(s2)) + // ... and add all necessary constraints + let mut constraint: Constraint; + for primary_substation_id in input.dg_kernel_for_switch_selection().iter_v() { + // * u(s) = + { + // Possible values of the square voltage at a primary substation. + // Tap positions are: T = {-10, ..., 10}. + // Base voltage: B = {1 + 0.1 * t | t \in T}. + // Square base voltage: {u² | u \in B}; for each t \in T the squared base + // voltage corresponding to t is BASE_VOLTAGE_SQ[t + 10]. + const BASE_VOLTAGE_SQ: [f64; 21] = [0.81, 0.8281, 0.8464, 0.8649, 0.8836, 0.9025, 0.9216, 0.9409, 0.9604, 0.9801, 1.0, + 1.0201, 1.0404, 1.0609, 1.0816, 1.1025, 1.1236, 1.1449, 1.1664, 1.1881, 1.21]; + constraint = Constraint::new(ConstraintType::Eq, 0.0, format!("u({})", primary_substation_id)); + constraint.add_wvar(WeightedVariable::new_idx(variables[&format!("u({})", primary_substation_id)], 1.0)); + for tap_position in -10..=10i8 { + constraint.add_wvar(WeightedVariable::new_idx(variables[&format!("tap({},{})", tap_position, primary_substation_id)], -BASE_VOLTAGE_SQ[(tap_position + 10) as usize])); + } + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // * tap(-10, s) + ... + tap(10, s) = 1 + { + constraint = Constraint::new(ConstraintType::Eq, 1.0, format!("sum_tap({})", primary_substation_id)); + for tap_position in -10..=10i8 { + constraint.add_wvar(WeightedVariable::new_idx(variables[&format!("tap({},{})", tap_position, primary_substation_id)], 1.0)); + } + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // * tap_abs(s) = 1 * (tap(-1, s) + tap(1, s)) + ... + 10 * (tap(-10, s) + tap(10, s)) + { + constraint = Constraint::new(ConstraintType::Eq, 0.0, format!("tap_abs({})", primary_substation_id)); + constraint.add_wvar(WeightedVariable::new_idx(variables[&format!("tap_abs({})", primary_substation_id)], 1.0)); + for tap_position in -10..=10i8 { + constraint.add_wvar(WeightedVariable::new_idx(variables[&format!("tap({},{})", tap_position, primary_substation_id)], -(tap_position.abs()) as f64)); + } + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // * max_tap_abs >= tap_abs(s) + { + constraint = Constraint::new(ConstraintType::GreaterThanEq, 0.0, format!("max_tap_abs({})", primary_substation_id)); + constraint.add_wvar(WeightedVariable::new_idx(variables[&"max_tap_abs".to_string()], 1.0)); + constraint.add_wvar(WeightedVariable::new_idx(variables[&format!("tap_abs({})", primary_substation_id)], -1.0)); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + for adjacent_id in input.iter_adjacent(&primary_substation_id).unwrap() { + let endpoints = input.e_attrs(&primary_substation_id, &adjacent_id, &0).unwrap().line_endpoints.unwrap(); + if input.v_attrs(&adjacent_id).unwrap().line_endpoints.is_some() && endpoints.1 != primary_substation_id { + // Collect the line + // We can guarantee here that line.len() >= 3 + let mut line: Vec = vec![primary_substation_id, adjacent_id]; + while *line.last().unwrap() != endpoints.1 { + for line_neighbour_id in input.iter_adjacent(line.last().unwrap()).unwrap() { + if line_neighbour_id == endpoints.1 + || line[line.len() - 2] != line_neighbour_id { + line.push(line_neighbour_id); + break; + } + } + } + // Process all secondary substations on the line + for substation_i1 in 1..=(line.len() - 2) { + // * u(s_j) = u_right(s_j) + sum_{k = 1}^j [ ( - r(s_j, s_j+1) p(s_k) + x(s_j, s_j+1) q(s_k) ) right_part(s_k, s_j) ] + // + u_left(s_j) + sum_{k = j}^m [ ( - r(s_j-1, s_j) p(s_k) + x(s_j-1, s_j) q(s_k) ) left_part(s_j, s_k) ] + { + constraint = Constraint::new( + ConstraintType::Eq, + 0, + format!("powerbalance({})", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u({})", line[substation_i1])], + -1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_right({})", line[substation_i1])], + 1.0 + )); + // sum for u_right(s_j) + for substation_i2 in 1..=substation_i1 { + variables.insert( + format!("right_part({},{})", line[substation_i2], line[substation_i1]), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Integer, 0.0, 0.0, 1.0, format!("right_part({},{})", line[substation_i2], line[substation_i1])))) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("right_part({},{})", line[substation_i2], line[substation_i1])], + - input.e_attrs(&line[substation_i1], &line[substation_i1 + 1], &0).unwrap().r + * input.v_attrs(&line[substation_i2]).unwrap().p + + input.e_attrs(&line[substation_i1], &line[substation_i1 + 1], &0).unwrap().x + * input.v_attrs(&line[substation_i2]).unwrap().q + )); + } + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_left({})", line[substation_i1])], + 1.0 + )); + // sum for u_left(s_j) + for substation_i2 in substation_i1..=(line.len() - 2) { + variables.insert( + format!("left_part({},{})", line[substation_i1], line[substation_i2]), + cplex_unwrap!(solver.get_problem_mut().add_variable(Variable::new(VariableType::Integer, 0.0, 0.0, 1.0, format!("left_part({},{})", line[substation_i2], line[substation_i1])))) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("left_part({},{})", line[substation_i1], line[substation_i2])], + - input.e_attrs(&line[substation_i1 - 1], &line[substation_i1], &0).unwrap().r + * input.v_attrs(&line[substation_i2]).unwrap().p + + input.e_attrs(&line[substation_i1 - 1], &line[substation_i1], &0).unwrap().x + * input.v_attrs(&line[substation_i2]).unwrap().q + )); + } + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // * u_right(s_j) = part(s_j) * u(s_j+1), which is linearised as + // ----* u_right(s_j) >= 0.81 * part(s_j) + { + constraint = Constraint::new( + ConstraintType::GreaterThanEq, + 0.0, + format!("u_right({})_lin1", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_right({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + -0.81 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* u_right(s_j) <= 1.21 * part(s_j) + { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 0.0, + format!("u_right({})_lin2", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_right({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + -1.21 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* u_right(s_j) <= u(s_j+1) - 2 part(s_j) + 2 + { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 2.0, + format!("u_right({})_lin3", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_right({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u({})", line[substation_i1 + 1])], + -1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + 2.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* u_right(s_j) >= u(s_j+1) + 2 part(s_j) - 2 + { + constraint = Constraint::new( + ConstraintType::GreaterThanEq, + -2.0, + format!("u_right({})_lin4", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_right({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u({})", line[substation_i1 + 1])], + -1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + -2.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // * u_left(s_j) = (1 - part(s_j)) * u(s_j-1), which is linearised as + // ----* u_left(s_j) >= -0.81 * part(s_j) + 0.81 + { + constraint = Constraint::new( + ConstraintType::GreaterThanEq, + 0.81, + format!("u_left({})_lin1", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_left({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + 0.81 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* u_left(s_j) <= -1.21 * part(s_j) + 1.21 + { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 1.21, + format!("u_left({})_lin2", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_left({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + 1.21 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* u_left(s_j) <= u(s_j-1) + 2 * part(s_j) + { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 0.0, + format!("u_left({})_lin3", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_left({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u({})", line[substation_i1 - 1])], + -1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + -2.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* u_left(s_j) >= u(s_j-1) - 2 * part(s_j) + { + constraint = Constraint::new( + ConstraintType::GreaterThanEq, + 0.0, + format!("u_left({})_lin4", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u_left({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("u({})", line[substation_i1 - 1])], + -1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + 2.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // * part(s_j-1) <= part(s_j) + if substation_i1 > 1 { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 0.0, + format!("part({})", line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1 - 1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + -1.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // Constraints for right_part(s_k, s_j) + for substation_i2 in 1..=substation_i1 { + // * right_part(s_k, s_j) = part(s_k) * part(s_j), which is liearised as + // ----* right_part(s_k, s_j) <= part(s_k) + { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 0.0, + format!("right_part({},{})_lin1", line[substation_i2], line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("right_part({},{})", line[substation_i2], line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i2])], + -1.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* right_part(s_k, s_j) <= part(s_j) + { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 0.0, + format!("right_part({},{})_lin2", line[substation_i2], line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("right_part({},{})", line[substation_i2], line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + -1.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* right_part(s_k, s_j) >= part(s_k) + part(s_j) - 1 + { + constraint = Constraint::new( + ConstraintType::GreaterThanEq, + -1.0, + format!("right_part({},{})_lin3", line[substation_i2], line[substation_i1]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("right_part({},{})", line[substation_i2], line[substation_i1])], + 1.0 + )); + if substation_i1 == substation_i2 { + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + -2.0 + )); + } else { + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i2])], + -1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + -1.0 + )); + } + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + } + // Constraints for left_part(s_j, s_k) + for substation_i2 in substation_i1..=(line.len() - 2) { + // * left_part(s_j, s_k) = (1 - part(s_j)) * (1 - part(s_k)), which is linearised as + // ----* left_part(s_j, s_k) <= 1 - part(s_j) + { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 1.0, + format!("left_part({},{})_lin1", line[substation_i1], line[substation_i2]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("left_part({},{})", line[substation_i1], line[substation_i2])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + 1.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* left_part(s_j, s_k) <= 1 - part(s_k) + { + constraint = Constraint::new( + ConstraintType::LessThanEq, + 1.0, + format!("left_part({},{})_lin2", line[substation_i1], line[substation_i2]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("left_part({},{})", line[substation_i1], line[substation_i2])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i2])], + 1.0 + )); + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + // ----* left_part(s_j, s_k) >= - part(s_j) - part(s_k) + 1 + { + constraint = Constraint::new( + ConstraintType::GreaterThanEq, + 1.0, + format!("left_part({},{})_lin3", line[substation_i1], line[substation_i2]) + ); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("left_part({},{})", line[substation_i1], line[substation_i2])], + 1.0 + )); + if substation_i1 == substation_i2 { + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + 2.0 + )); + } else { + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i1])], + 1.0 + )); + constraint.add_wvar(WeightedVariable::new_idx( + variables[&format!("part({})", line[substation_i2])], + 1.0 + )); + } + cplex_unwrap!(solver.get_problem_mut().add_constraint(constraint)); + } + } + } + } + } + } + solver.variables = variables; + Ok(solver) + } + + fn get_solution(&self) -> Option<(SwitchSelectionGraph, TapValue)> { + if self.solution.is_none() { + return None; + } + let mut answer = self.input.unwrap().clone(); + for primary_substation_id in self.input.dg_kernel_for_switch_selection().iter_v() { + for tap_position in -10..10i8 { + if let VariableValue::Integer(1) = self.solution.as_ref().unwrap().variables[self.variables[&format!("tap({},{})", tap_position, primary_substation_id)]] { + answer.v_attrs_mut(&primary_substation_id).unwrap().tap_position = Some(tap_position); + } + } + } + for edge in self.input.iter_e() { + let endpoints = self.input.e_attrs(&edge.id1, &edge.id2, &0).unwrap().line_endpoints.unwrap(); + if self.input.v_attrs(&edge.id1).unwrap().tap_position.is_some() { + if self.input.v_attrs(&edge.id2).unwrap().tap_position.is_some() { + answer.e_attrs_mut(&edge.id1, &edge.id2, &0).unwrap().switch = true; + continue; + } + if let VariableValue::Integer(value) = self.solution.as_ref().unwrap().variables[self.variables[&format!("part({})", edge.id2)]] { + answer.e_attrs_mut(&edge.id1, &edge.id2, &0).unwrap().switch = endpoints.0 == edge.id1 && value == 1 || endpoints.1 == edge.id1 && value == 0; + } + continue; + } + if self.input.v_attrs(&edge.id2).unwrap().tap_position.is_some() { + if let VariableValue::Integer(value) = self.solution.as_ref().unwrap().variables[self.variables[&format!("part({})", edge.id1)]] { + answer.e_attrs_mut(&edge.id1, &edge.id2, &0).unwrap().switch = endpoints.0 == edge.id2 && value == 1 || endpoints.1 == edge.id2 && value == 0; + } + continue; + } + if let VariableValue::Integer(value1) = self.solution.as_ref().unwrap().variables[self.variables[&format!("part({})", edge.id1)]] { + if let VariableValue::Integer(value2) = self.solution.as_ref().unwrap().variables[self.variables[&format!("part({})", edge.id2)]] { + answer.e_attrs_mut(&edge.id1, &edge.id2, &0).unwrap().switch = value1 != value2; + } + } + } + Some((answer, self.solution.as_ref().unwrap().objective as TapValue)) + } + + fn solve(&mut self) -> Result<(), SolverError> { + self.solution = Some(cplex_unwrap!(self.get_problem_mut().solve(ProblemType::MixedInteger))); + Ok(()) + } +} diff --git a/src/solver/errors.rs b/src/solver/errors.rs new file mode 100644 index 0000000..11f93e3 --- /dev/null +++ b/src/solver/errors.rs @@ -0,0 +1,63 @@ +use std::{error::Error, fmt::Display}; + + + + + +#[derive(Debug)] +pub struct GraphError { + description: String, +} + +// GraphError::GraphError +impl GraphError { + #[inline] + pub fn from_str(description: &str) -> Self { + GraphError { description: description.to_string() } + } + + #[inline] + pub fn from_string(description: String) -> Self { + GraphError { description } + } +} + +// GraphError::Error +impl Error for GraphError {} + +// GraphError::Display +impl Display for GraphError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.description) + } +} + + + +#[derive(Debug)] +pub struct SolverError { + description: String, +} + +// SolverError::SolverError +impl SolverError { + #[inline] + pub fn from_str(description: &str) -> Self { + SolverError { description: description.to_string() } + } + + #[inline] + pub fn from_string(description: String) -> Self { + SolverError { description } + } +} + +// SolverError::Error +impl Error for SolverError {} + +// SolverError::Display +impl Display for SolverError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.description) + } +} diff --git a/src/solver/mod.rs b/src/solver/mod.rs new file mode 100644 index 0000000..a8ea7e2 --- /dev/null +++ b/src/solver/mod.rs @@ -0,0 +1,5 @@ +pub mod base_solver; +pub mod benchmark; +pub mod cplex_solver; +pub mod tree_decomposition_solver; +pub mod errors; diff --git a/src/solver/tree_decomposition_solver.rs b/src/solver/tree_decomposition_solver.rs new file mode 100644 index 0000000..b1ef57d --- /dev/null +++ b/src/solver/tree_decomposition_solver.rs @@ -0,0 +1,404 @@ +use std::{collections::{HashMap, HashSet, VecDeque}, sync::{mpsc::{self, Receiver, Sender}, Arc, Mutex}, thread::{self, JoinHandle}}; +use crabnets::{BasicImmutableGraph, BasicMutableGraph, ImmutableGraphContainer}; +use itertools::Itertools; +use crate::{switch_selection_instance::{SwitchSelectionInstance, SwitchSelectionGraph}, tree_decomposition::TreeDecomposition}; +use super::{base_solver::*, errors::SolverError}; + + + + + +#[derive(Clone)] +struct TapsMemo { + pub primary_substations: Vec, + table: HashMap, TapValue>, +} + +// TapsMemo::TapsMemo +impl TapsMemo { + pub fn complete(primary_substations: Vec) -> TapsMemo { + TapsMemo { + primary_substations: primary_substations.clone(), + table: HashMap::from_iter( + primary_substations + .iter() + .map(|_| (-10..=10i8)) + .multi_cartesian_product() + .map(|x| + (x.clone(), x.iter().map(|y| y.abs()).max().unwrap()) + ) + ) + } + } + + #[inline] + pub fn empty(primary_substations: Vec) -> TapsMemo { + TapsMemo { primary_substations: primary_substations.clone(), table: HashMap::with_capacity(21usize.pow(primary_substations.len() as u32)) } + } + + pub fn intersect(&mut self, other: &TapsMemo) { + let common_primary_substations_self_indices = self.primary_substations + .iter() + .enumerate() + .filter(|&(_, x)| other.primary_substations.binary_search(x).is_ok()) + .map(|(x, _)| x) + .collect_vec(); + let common_primary_substations_other_indices = common_primary_substations_self_indices + .iter() + .map(|&x| other.primary_substations.binary_search(&self.primary_substations[x])) + .filter(|x| x.is_ok()) + .map(|x| x.unwrap()) + .collect_vec(); + // Decide what to keep and what to remove + // We only keep entries that have at least one corresponding entry in + // other.table. + // Among all the corresponding entries, we choose one with the lowest + // value of the objective function and then choose the maximum between + // the value stored in that entry and the value stored in the entry + // from self.table. + let mut optimal_corresponding_entries = HashMap::new(); + for (other_taps_positions, &other_obj_value) in other.table.iter() { + let other_taps_positions_for_common_primary_substations = common_primary_substations_other_indices + .iter() + .map(|&x| other_taps_positions[x]) + .collect_vec(); + match optimal_corresponding_entries.get_mut(&other_taps_positions_for_common_primary_substations) { + Some(value) => if other_obj_value < *value { + *value = other_obj_value; + }, + None => { optimal_corresponding_entries.insert(other_taps_positions_for_common_primary_substations, other_obj_value); }, + } + } + let mut entries_to_be_removed = Vec::new(); + for (taps_positions, obj_value) in self.table.iter_mut() { + let taps_positions_for_common_primary_substations = common_primary_substations_self_indices + .iter() + .map(|&x| taps_positions[x]) + .collect_vec(); + match optimal_corresponding_entries.get(&taps_positions_for_common_primary_substations) { + Some(&value) => if *obj_value < value { + *obj_value = value; + }, + None => entries_to_be_removed.push(taps_positions.clone()), + } + } + for taps_position in entries_to_be_removed.drain(..) { + self.table.remove(&taps_position); + } + } +} + + + +fn locally_feasible_taps_positions(input: Arc, bag: &Vec) -> TapsMemo { + // Possible values of the square voltage at a primary substation. + // Tap positions are: T = {-10, ..., 10}. + // Base voltage: B = {1 + 0.1 * t | t \in T}. + // Square base voltage: {u² | u \in B}; for each t \in T the squared base + // voltage corresponding to t is BASE_VOLTAGE_SQ[t + 10]. + const BASE_VOLTAGE_SQ: [f64; 21] = [0.81, 0.8281, 0.8464, 0.8649, 0.8836, 0.9025, 0.9216, 0.9409, 0.9604, 0.9801, 1.0, + 1.0201, 1.0404, 1.0609, 1.0816, 1.1025, 1.1236, 1.1449, 1.1664, 1.1881, 1.21]; + let mut answer = TapsMemo::complete(bag.clone()); + // Consider all possible pairs of primary substations from the bag. If + // there're lines between a pair of the primary substations, try cutting + // each line in different places and see which tap positions are feasible. + for left_primary_substation_i in 0..bag.len() { + let left_primary_substation_id = bag[left_primary_substation_i]; + for right_primary_substation_i in left_primary_substation_i..bag.len() { + let right_primary_substation_id = bag[right_primary_substation_i]; + for adjacent_id in input.iter_adjacent(&left_primary_substation_id).unwrap() { + // If the adjacent vertex is a secondary substation lying on a line between + // left_primary_substation_id and right_primary_substation_id, reconstruct the + // entire line. + match input.v_attrs(&adjacent_id).unwrap().line_endpoints { + Some(value) => if value != (left_primary_substation_id, right_primary_substation_id) { + continue; + }, + None => continue, + } + let mut line = Vec::from([left_primary_substation_id, adjacent_id]); + loop { + let last_discovered_substation_id = line.last().unwrap(); + if input.v_attrs(last_discovered_substation_id).unwrap().tap_position.is_some() { + break; + } + for adjacent_to_last_id in input.iter_adjacent(last_discovered_substation_id).unwrap() { + if adjacent_to_last_id != line[line.len() - 2] { + line.push(adjacent_to_last_id); + break; + } + } + } + // Now that we have a full line, cut all possible edges on it one after another + // and compute feasible tap positions for every cut. + let mut line_memo = TapsMemo::empty(bag.clone()); + for (last_left_substation_i, first_right_substation_i) in (0..(line.len() - 1)).zip(1..line.len()) { + let left_line = &line[..=last_left_substation_i]; + let right_line = &line[first_right_substation_i..]; + let mut voltage_sq = 1.0; + let mut voltage_sq_peak: f64 = 1.0; + let mut voltage_sq_gorge: f64 = 1.0; + for left_substation_i in 1..left_line.len() { + voltage_sq += input.e_attrs(&left_line[left_substation_i - 1], &left_line[left_substation_i], &0).unwrap().x + * left_line[left_substation_i..].iter().map(|x| input.v_attrs(x).unwrap().q).sum::() + - input.e_attrs(&left_line[left_substation_i - 1], &left_line[left_substation_i], &0).unwrap().r + * left_line[left_substation_i..].iter().map(|x| input.v_attrs(x).unwrap().p).sum::(); + voltage_sq_peak = voltage_sq_peak.max(voltage_sq); + voltage_sq_gorge = voltage_sq_gorge.min(voltage_sq); + } + if voltage_sq_peak - voltage_sq_gorge > 1.21 - 0.81 { + continue; + } + let left_tap_position_min = BASE_VOLTAGE_SQ.iter().enumerate().filter(|&(_, &x)| x >= 1.81 - voltage_sq_gorge).next().unwrap().0 as TapValue - 10; + let left_tap_position_max = BASE_VOLTAGE_SQ.iter().enumerate().rev().filter(|&(_, &x)| x <= 2.21 - voltage_sq_peak).next().unwrap().0 as TapValue - 10; + voltage_sq = 1.0; + voltage_sq_peak = 1.0; + voltage_sq_gorge = 1.0; + for right_substation_i in (0..(right_line.len() - 1)).rev() { + voltage_sq += input.e_attrs(&right_line[right_substation_i + 1], &right_line[right_substation_i], &0).unwrap().x + * right_line[..=right_substation_i].iter().map(|x| input.v_attrs(x).unwrap().q).sum::() + - input.e_attrs(&right_line[right_substation_i + 1], &right_line[right_substation_i], &0).unwrap().r + * right_line[..=right_substation_i].iter().map(|x| input.v_attrs(x).unwrap().p).sum::(); + voltage_sq_peak = voltage_sq_peak.max(voltage_sq); + voltage_sq_gorge = voltage_sq_gorge.min(voltage_sq); + } + if voltage_sq_peak - voltage_sq_gorge > 1.21 - 0.81 { + continue; + } + let right_tap_position_min = BASE_VOLTAGE_SQ.iter().enumerate().filter(|&(_, &x)| x >= 1.81 - voltage_sq_gorge).next().unwrap().0 as TapValue - 10; + let right_tap_position_max = BASE_VOLTAGE_SQ.iter().enumerate().rev().filter(|&(_, &x)| x <= 2.21 - voltage_sq_peak).next().unwrap().0 as TapValue - 10; + // The following operation can take almost 90% of all computation time!!! Can + // be optimised by, e.g. replacing these memos with DataFrames and using joins + // instead of extending every line_memo with dozens of rows. + line_memo.table.extend(answer.table.iter().filter(|&(x, _)| + x[left_primary_substation_i] >= left_tap_position_min && + x[left_primary_substation_i] <= left_tap_position_max && + x[right_primary_substation_i] >= right_tap_position_min && + x[right_primary_substation_i] <= right_tap_position_max + ).map(|(k, v)| (k.clone(), v.clone()))); + } + answer = line_memo; + } + } + } + answer +} + +fn thread_workload(input: Arc, memos: Arc>>, td: Arc, bag_id: usize, rx: Receiver) -> Result<(), SolverError> { + // Create a memo for this bag + let bag = td.v_attrs(&bag_id).unwrap().vertices.clone(); + let mut memo = locally_feasible_taps_positions(input, &bag); + // Intersect memo with the memos of the children + let mut remaining_children: HashSet = td.iter_adjacent_out(&bag_id).unwrap().collect(); + while !remaining_children.is_empty() { + let received_bag_id = match rx.recv() { + Ok(value) => value, + Err(_) => return Ok(()), + }; + if remaining_children.contains(&received_bag_id) { + let child_memo = memos.lock().unwrap()[&received_bag_id].clone(); + memo.intersect(&child_memo); + remaining_children.remove(&received_bag_id); + } + } + // If memo is empty, the instance is infeasible + if memo.table.is_empty() { + return Err(SolverError::from_str("TreeDecompositionSolver. The problem instance is infeasible.")); + } + // Otherwise, save in the memos collection and send it to the parent + memos.lock().unwrap().insert(bag_id, memo); + Ok(()) +} + +fn solution_graph_setup(solution: &mut SwitchSelectionGraph, dg_kernel: &SwitchSelectionGraph, taps_positions: &HashMap) { + // Possible values of the square voltage at a primary substation. + // Tap positions are: T = {-10, ..., 10}. + // Base voltage: B = {1 + 0.1 * t | t \in T}. + // Square base voltage: {u² | u \in B}; for each t \in T the squared base + // voltage corresponding to t is BASE_VOLTAGE_SQ[t + 10]. + const BASE_VOLTAGE_SQ: [f64; 21] = [0.81, 0.8281, 0.8464, 0.8649, 0.8836, 0.9025, 0.9216, 0.9409, 0.9604, 0.9801, 1.0, + 1.0201, 1.0404, 1.0609, 1.0816, 1.1025, 1.1236, 1.1449, 1.1664, 1.1881, 1.21]; + for (left_primary_substation_id, right_primary_substation_id) in dg_kernel.iter_e().map(|x| (x.id1, x.id2)) { + for adjacent_id in solution.iter_adjacent(&left_primary_substation_id).unwrap().collect_vec() { + // If the adjacent vertex is a secondary substation lying on a line between + // left_primary_substation_id and right_primary_substation_id, reconstruct the + // entire line. + match solution.v_attrs(&adjacent_id).unwrap().line_endpoints { + Some(value) => if value != (left_primary_substation_id, right_primary_substation_id) { + continue; + }, + None => continue, + } + let mut line = Vec::from([left_primary_substation_id, adjacent_id]); + loop { + let last_discovered_substation_id = line.last().unwrap(); + if solution.v_attrs(last_discovered_substation_id).unwrap().tap_position.is_some() { + break; + } + for adjacent_to_last_id in solution.iter_adjacent(last_discovered_substation_id).unwrap() { + if adjacent_to_last_id != line[line.len() - 2] { + line.push(adjacent_to_last_id); + break; + } + } + } + // Now that we have a full line, cut all possible edges on it one after another + // until a feasible cut is found. + for (last_left_substation_i, first_right_substation_i) in (0..(line.len() - 1)).zip(1..line.len()) { + let left_line = &line[..=last_left_substation_i]; + let right_line = &line[first_right_substation_i..]; + let mut voltage_sq = BASE_VOLTAGE_SQ[(taps_positions[&left_primary_substation_id] + 10) as usize]; + let mut voltage_sq_peak = voltage_sq; + let mut voltage_sq_gorge = voltage_sq; + for left_substation_i in 1..left_line.len() { + voltage_sq += solution.e_attrs(&left_line[left_substation_i - 1], &left_line[left_substation_i], &0).unwrap().x + * left_line[left_substation_i..].iter().map(|x| solution.v_attrs(x).unwrap().q).sum::() + - solution.e_attrs(&left_line[left_substation_i - 1], &left_line[left_substation_i], &0).unwrap().r + * left_line[left_substation_i..].iter().map(|x| solution.v_attrs(x).unwrap().p).sum::(); + voltage_sq_peak = voltage_sq_peak.max(voltage_sq); + voltage_sq_gorge = voltage_sq_gorge.min(voltage_sq); + } + if voltage_sq_peak > 1.21 || voltage_sq_gorge < 0.81 { + continue; + } + voltage_sq = BASE_VOLTAGE_SQ[(taps_positions[&right_primary_substation_id] + 10) as usize]; + voltage_sq_peak = voltage_sq; + voltage_sq_gorge = voltage_sq; + for right_substation_i in (0..(right_line.len() - 1)).rev() { + voltage_sq += solution.e_attrs(&right_line[right_substation_i + 1], &right_line[right_substation_i], &0).unwrap().x + * right_line[..=right_substation_i].iter().map(|x| solution.v_attrs(x).unwrap().q).sum::() + - solution.e_attrs(&right_line[right_substation_i + 1], &right_line[right_substation_i], &0).unwrap().r + * right_line[..=right_substation_i].iter().map(|x| solution.v_attrs(x).unwrap().p).sum::(); + voltage_sq_peak = voltage_sq_peak.max(voltage_sq); + voltage_sq_gorge = voltage_sq_gorge.min(voltage_sq); + } + if voltage_sq_peak > 1.21 || voltage_sq_gorge < 0.81 { + continue; + } + // If we reach this point, we've found a feasible cut. Record it in the graph. + solution.e_attrs_mut(&line[last_left_substation_i], &line[first_right_substation_i], &0).unwrap().switch = true; + break; + } + } + } + for primary_substation_id in dg_kernel.iter_v() { + solution.v_attrs_mut(&primary_substation_id).unwrap().tap_position = Some(taps_positions[&primary_substation_id]); + } +} + + + +struct ThreadMetadata { + bag_id: usize, + join_handle: Option>>, + tx: Option>, +} + + + +pub struct TreeDecompositionSolver { + dg_kernel: SwitchSelectionGraph, + input: Arc, + td: Arc, + memos: Option>, + thread_count: usize, +} + +// TreeDecompositionSolver::BaseSolver +impl BaseSolver for TreeDecompositionSolver { + fn with_input(input: SwitchSelectionInstance) -> Result { + let dg_kernel = input.dg_kernel_for_switch_selection(); + let td = match TreeDecomposition::for_switch_selection_graph(&dg_kernel) { + Ok(value) => value, + Err(value) => return Err(SolverError::from_string(value.to_string())), + }; + Ok(TreeDecompositionSolver { + dg_kernel, + input: Arc::new(input), + td: Arc::new(td), + memos: None, + thread_count: num_cpus::get() - 1, + }) + } + + fn get_solution(&self) -> Option<(SwitchSelectionGraph, TapValue)> { + if self.memos.is_none() { + return None; + } + let memos = self.memos.as_ref().unwrap(); + let mut answer = self.input.as_ref().unwrap().clone(); + // Tap positions will be gradually collected in taps_positions + let mut taps_positions: HashMap = HashMap::new(); + // Traverse the tree decomposition top-down in breadth-first search + // pre-ordering to construct a solution + let mut bag_queue: VecDeque = VecDeque::from([self.td.root_id]); + while !bag_queue.is_empty() { + // Get the ID of the current bag. + let curr_bag_id = bag_queue.pop_front().unwrap(); + // Find the best locally feasible solution that agrees with what + // was already built in taps_positions. + let curr_entry = memos[&curr_bag_id].table + .iter() + .sorted_by_key(|&(_, &x)| x) + .find(|&(k, _)| + k.iter().enumerate().all(|(i, &x)| + match taps_positions.get(&memos[&curr_bag_id].primary_substations[i]) { + Some(&value) => x == value, + None => true, + } + ) + ) + .unwrap() + .0 + .iter() + .enumerate() + .map(|(i, &x)| (memos[&curr_bag_id].primary_substations[i], x)); + taps_positions.extend(curr_entry); + bag_queue.extend(self.td.iter_adjacent_out(&curr_bag_id).unwrap()); + } + solution_graph_setup(&mut answer, &self.dg_kernel, &taps_positions); + Some((answer, memos[&self.td.root_id].table.iter().sorted_by_key(|&(_, &x)| x).next().unwrap().1.clone())) + } + + fn solve(&mut self) -> Result<(), SolverError> { + // Memos + let memos: Arc>> = Arc::new(Mutex::new(HashMap::new())); + // Find out the depth-first search postordering for the bags of self.td + let mut thread_data = self.td.dfs_postordering().into_iter().map(|id: usize| ThreadMetadata { bag_id: id, join_handle: None, tx: None }).collect_vec(); + // Launch threads with the sliding window + let mut left_bound: usize = 0; + let mut right_bound = self.thread_count.min(thread_data.len()) - 1; + for bag_i in left_bound..=right_bound { + let input_clone: Arc = self.input.clone(); + let memos_clone: Arc>> = memos.clone(); + let dtd_clone: Arc = self.td.clone(); + let bag_id_clone: usize = thread_data[bag_i].bag_id; + let (tx, rx) = mpsc::channel(); + thread_data[bag_i].join_handle = Some(thread::spawn(move || thread_workload(input_clone, memos_clone, dtd_clone, bag_id_clone, rx))); + thread_data[bag_i].tx = Some(tx); + } + while left_bound <= right_bound { + thread_data[left_bound].join_handle.take().unwrap().join().unwrap()?; + left_bound += 1; + for bag_i in left_bound..=right_bound { + thread_data[bag_i].tx.as_ref().unwrap().send(thread_data[left_bound - 1].bag_id).unwrap_or(()); + } + if right_bound + 1 < thread_data.len() { + right_bound += 1; + let input_clone: Arc = self.input.clone(); + let memos_clone: Arc>> = memos.clone(); + let dtd_clone: Arc = self.td.clone(); + let bag_id_clone: usize = thread_data[right_bound].bag_id; + let (tx, rx) = mpsc::channel(); + thread_data[right_bound].join_handle = Some(thread::spawn(move || thread_workload(input_clone, memos_clone, dtd_clone, bag_id_clone, rx))); + thread_data[right_bound].tx = Some(tx); + for bag_i in 0..left_bound { + thread_data[right_bound].tx.as_ref().unwrap().send(thread_data[bag_i].bag_id).unwrap(); + } + } + } + // Save the memos + self.memos = Some(Arc::into_inner(memos).unwrap().into_inner().unwrap()); + Ok(()) + } +} diff --git a/src/switch_selection_instance.rs b/src/switch_selection_instance.rs new file mode 100644 index 0000000..62d397e --- /dev/null +++ b/src/switch_selection_instance.rs @@ -0,0 +1,226 @@ +use std::{cmp::Ordering, collections::VecDeque, iter::once}; +use crabnets::{attributes::*, io::{AttributeCollectionIO, AttributeToken}, locales::*, topology_tests::TopologyTests, *}; +use itertools::Itertools; +use crate::solver::errors::GraphError; + + + + + + +#[derive(Clone, Default)] +pub struct DGVertexAttributes { + pub line_endpoints: Option<(usize, usize)>, + pub p: f64, + pub q: f64, + pub tap_position: Option, +} + +// DGVertexAttributes::AttributeCollection +impl AttributeCollection for DGVertexAttributes { + fn new() -> Self { + DGVertexAttributes { line_endpoints: None, p: 0.0, q: 0.0, tap_position: None } + } +} + +// DGVertexAttributes::AttributeCollectionIO +impl AttributeCollectionIO for DGVertexAttributes { + fn io_iter_contents<'a>(&'a self) -> Box> + 'a> { + let tap_position_data = self.tap_position.is_some().then( + || once(AttributeToken { name: "tap position", value: StaticDispatchAttributeValue::Int8(self.tap_position.unwrap()) }) + ).into_iter().flatten(); + Box::new( + once(AttributeToken { name: "p", value: StaticDispatchAttributeValue::Float64(self.p) }) + .chain(once(AttributeToken { name: "q", value: StaticDispatchAttributeValue::Float64(self.q) })) + .chain(tap_position_data) + ) + } + + fn io_query_contents(&self, attribute_name: &str) -> Option { + match attribute_name { + "p" => Some(StaticDispatchAttributeValue::Float64(self.p)), + "q" => Some(StaticDispatchAttributeValue::Float64(self.q)), + "tap position" => match self.tap_position { + Some(value) => Some(StaticDispatchAttributeValue::Int8(value)), + None => None, + }, + _ => None, + } + } + + fn io_reader_callback<'a, EdgeIdType, VertexIdType>(&mut self, token: AttributeToken<'a>) + where + EdgeIdType: Id, + VertexIdType: Id + { + match token.name { + "p" => if let StaticDispatchAttributeValue::Float64(value) = token.value { + self.p = value; + }, + "q" => if let StaticDispatchAttributeValue::Float64(value) = token.value { + self.q = value; + }, + "is primary substation" => if let StaticDispatchAttributeValue::Bool(value) = token.value { + self.tap_position = if value { + Some(0) + } else { + None + } + }, + _ => (), + } + } +} + + + +#[derive(Clone, Default)] +pub struct DGEdgeAttributes { + pub line_endpoints: Option<(usize, usize)>, + pub r: f64, + pub switch: bool, + pub x: f64, +} + +// DGEdgeAttributes::AttributeCollection +impl AttributeCollection for DGEdgeAttributes { + fn new() -> Self { + DGEdgeAttributes { line_endpoints: None, r: 0.0, switch: false, x: 0.0 } + } +} + +// DGEdgeAttributes::AttributeCollectionIO +impl AttributeCollectionIO for DGEdgeAttributes { + fn io_iter_contents<'a>(&'a self) -> Box> + 'a> { + Box::new( + once(AttributeToken { name: "r", value: StaticDispatchAttributeValue::Float64(self.r) }) + .chain(once(AttributeToken { name: "opened switch", value: StaticDispatchAttributeValue::Bool(self.switch) })) + .chain(once(AttributeToken { name: "x", value: StaticDispatchAttributeValue::Float64(self.x) })) + ) + } + + fn io_query_contents(&self, attribute_name: &str) -> Option { + match attribute_name { + "r" => Some(StaticDispatchAttributeValue::Float64(self.r)), + "opened switch" => Some(StaticDispatchAttributeValue::Bool(self.switch)), + "x" => Some(StaticDispatchAttributeValue::Float64(self.x)), + _ => None, + } + } + + fn io_reader_callback<'a, EdgeIdType, VertexIdType>(&mut self, token: AttributeToken<'a>) + where + EdgeIdType: Id, + VertexIdType: Id, + { + match token.name { + "r" => if let StaticDispatchAttributeValue::Float64(value) = token.value { + self.r = value; + }, + "x" => if let StaticDispatchAttributeValue::Float64(value) = token.value { + self.x = value; + }, + _ => (), + } + } +} + + + +pub type SwitchSelectionGraph = graph!(A ---A--- A with VertexAttributeCollectionType = DGVertexAttributes, EdgeAttributeCollectionType = DGEdgeAttributes); + + + +#[derive(Clone, Default)] +pub struct SwitchSelectionInstance { + graph: SwitchSelectionGraph, +} + +// SwitchSelectionInstance::SwitchSelectionInstance +impl SwitchSelectionInstance { + pub fn new(mut graph: SwitchSelectionGraph) -> Result { + if !graph.is_connected() { + return Err(GraphError::from_str("The given ditribution grid is not connected.")); + } + let mut unvisited_primary_substations: VecDeque = VecDeque::from_iter(graph.iter_v().filter(|x| graph.v_attrs(x).unwrap().tap_position.is_some())); + if unvisited_primary_substations.is_empty() { + return Err(GraphError::from_str("The given distribution grid doesn't contain any primary substations.")); + } + // Launch depth-first search from each primary substation to determine + // which secondary substation belongs to a line between which primary + // substations. + while !unvisited_primary_substations.is_empty() { + let primary_substation_id = unvisited_primary_substations.pop_front().unwrap(); + let mut unvisited_vertices_stack = VecDeque::from_iter( + graph.iter_adjacent(&primary_substation_id).unwrap().filter(|x| + graph.v_attrs(x).unwrap().line_endpoints.is_none() + ) + ); + let mut curr_line = Vec::from([primary_substation_id]); + while !unvisited_vertices_stack.is_empty() { + let curr_substation_id = unvisited_vertices_stack.pop_front().unwrap(); + curr_line.push(curr_substation_id); + // If the current substation is a primary substation, we've reached the end of the + // line. Backtrack and set line_endpoints to primary_substation_id and + // curr_substation_id. + if graph.v_attrs(&curr_substation_id).unwrap().tap_position.is_some() { + let line_endpoints = match primary_substation_id.cmp(&curr_substation_id) { + Ordering::Less => (primary_substation_id, curr_substation_id), + Ordering::Equal => return Err(GraphError::from_string(format!("Primary substation {} has a feeder that begins and ends in it.", curr_substation_id))), + Ordering::Greater => (curr_substation_id, primary_substation_id), + }; + graph.e_attrs_mut(curr_line.last().unwrap(), &curr_line[curr_line.len() - 2], &0).unwrap().line_endpoints = Some(line_endpoints.clone()); + for substation_i in (1..(curr_line.len() - 1)).rev() { + graph.v_attrs_mut(&curr_line[substation_i]).unwrap().line_endpoints = Some(line_endpoints.clone()); + graph.e_attrs_mut(&curr_line[substation_i], &curr_line[substation_i - 1], &0).unwrap().line_endpoints = Some(line_endpoints.clone()); + } + curr_line.resize(1, 0); + continue; + } + // If the current substation is a secondary substation, check that it has exactly 2 + // neighbours and add the unvisited one on top of the stack. + let adjacent_substations = graph.iter_adjacent(&curr_substation_id).unwrap().collect_vec(); + if adjacent_substations.len() != 2 { + return Err(GraphError::from_string(format!("Secondary substation {} must have exactly 2 adjacent substations.", curr_substation_id))); + } + unvisited_vertices_stack.push_front(if adjacent_substations[0] != curr_line[curr_line.len() - 2] { adjacent_substations[0] } else { adjacent_substations[1] }); + } + } + Ok(SwitchSelectionInstance { graph }) + } + + pub fn dg_kernel_for_switch_selection(&self) -> SwitchSelectionGraph { + let mut answer = SwitchSelectionGraph::new(); + let lines = self.graph.iter_e().map(|x| self.graph.e_attrs(&x.id1, &x.id2, &x.edge_id).unwrap().line_endpoints.unwrap()).unique(); + for line in lines { + if !answer.contains_v(&line.0) { + answer.add_v(Some(line.0)); + } + if !answer.contains_v(&line.1) { + answer.add_v(Some(line.1)); + } + answer.add_e(&line.0, &line.1, false, None).unwrap(); + } + answer + } +} + +// SwitchSelectionInstance::ImmutableGraphContainer +impl ImmutableGraphContainer for SwitchSelectionInstance { + type EdgeAttributeCollectionType = DGEdgeAttributes; + type EdgeIdType = u8; + type LocaleType = SimpleUndirectedLocale; + type VertexAttributeCollectionType = DGVertexAttributes; + type VertexIdType = usize; + + fn unwrap(&self) -> &Graph { + &self.graph + } +} + +// SwitchSelectionInstance::MutableGraphContainer +impl MutableGraphContainer for SwitchSelectionInstance { + fn unwrap(&mut self) -> &mut Graph { + &mut self.graph + } +} diff --git a/src/tree_decomposition.rs b/src/tree_decomposition.rs new file mode 100644 index 0000000..2bf70db --- /dev/null +++ b/src/tree_decomposition.rs @@ -0,0 +1,102 @@ +use std::collections::VecDeque; +use arboretum_td::{exact::TamakiPid, graph::{HashMapGraph, MutableGraph as ArboretumMutableGraph}, solver::{AtomSolver, ComputationResult}}; +use crabnets::{*, attributes::*, locales::*}; +use itertools::Itertools; +use crate::{switch_selection_instance::SwitchSelectionGraph, solver::errors::GraphError}; + + + + + +#[derive(Clone, Default)] +pub struct Bag { + pub vertices: Vec, +} + +// Bag::AttributeCollection +impl AttributeCollection for Bag { + fn new() -> Self { + Bag { vertices: Vec::new() } + } +} + + + +#[derive(Clone, Default)] +pub struct TreeDecomposition { + graph: graph!(A ---X--> A with VertexAttributeCollectionType = Bag), + pub max_bag_size: usize, + pub root_id: usize, +} + +// TreeDecomposition::TreeDecomposition +impl TreeDecomposition { + pub fn dfs_postordering(&self) -> Vec { + let mut answer: Vec = Vec::with_capacity(self.count_v()); + let mut bag_stack: VecDeque = VecDeque::from([self.root_id]); + while !bag_stack.is_empty() { + let last_bag = *bag_stack.back().unwrap(); + if self.v_degree_out(&last_bag).unwrap() > 0 && (answer.len() == 0 || !self.iter_adjacent_out(&last_bag).unwrap().contains(answer.last().unwrap())) { + bag_stack.extend(self.iter_adjacent_out(&last_bag).unwrap()); + continue; + } + answer.push(bag_stack.pop_back().unwrap()); + } + answer + } + + pub fn for_switch_selection_graph(graph: &SwitchSelectionGraph) -> Result { + let mut arboretum_graph = HashMapGraph::new(); + for id in graph.iter_v() { + arboretum_graph.add_vertex(id); + } + for edge in graph.iter_e() { + arboretum_graph.add_edge(edge.id1, edge.id2); + } + match TamakiPid::with_graph(&arboretum_graph).compute() { + ComputationResult::ComputedTreeDecomposition(td) => { + let mut answer: graph!(A ---X--> A with VertexAttributeCollectionType = Bag) = Graph::new(); + for bag in td.bags() { + let vertex_set = Vec::from_iter(bag.vertex_set.iter().sorted().cloned()); + answer.add_v(Some(bag.id)); + answer.v_attrs_mut(&bag.id).unwrap().vertices = vertex_set; + } + let tree_decomposition_root_id = answer.iter_v().next().unwrap(); + let mut unvisited_bags = VecDeque::from([tree_decomposition_root_id]); + while !unvisited_bags.is_empty() { + let curr_bag_id = unvisited_bags.pop_front().unwrap(); + for &adjacent_bag_id in td.bags[curr_bag_id].neighbors.iter() { + if answer.contains_e(&curr_bag_id, &adjacent_bag_id, &0).is_none() { + answer.add_e(&curr_bag_id, &adjacent_bag_id, true, None).unwrap(); + unvisited_bags.push_back(adjacent_bag_id); + } + } + } + Ok(TreeDecomposition { graph: answer, max_bag_size: td.max_bag_size, root_id: tree_decomposition_root_id }) + }, + ComputationResult::Bounds(bounds) => { + if bounds.lowerbound == graph.count_v() - 1 { + let mut answer: graph!(A ---X--> A with VertexAttributeCollectionType = Bag) = Graph::new(); + answer.add_v(Some(0)); + answer.v_attrs_mut(&0).unwrap().vertices = graph.iter_v().sorted().collect(); + Ok(TreeDecomposition { graph: answer, max_bag_size: graph.count_v(), root_id: 0 }) + } else { + Err(GraphError::from_str("Failed to compute a tree decomposition.")) + } + }, + } + } +} + +// TreeDecomposition::ImmutableGraphContainer +impl ImmutableGraphContainer for TreeDecomposition { + type EdgeAttributeCollectionType = (); + type EdgeIdType = u8; + type LocaleType = SimpleDirectedLocale<(), Bag, usize>; + type VertexAttributeCollectionType = Bag; + type VertexIdType = usize; + + fn unwrap(&self) -> &Graph { + &self.graph + } +}