diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 0000000..c9162b9 --- /dev/null +++ b/.editorconfig @@ -0,0 +1,12 @@ +root = true + +[*] +charset = utf-8 +indent_style = tab +indent_size = 4 +insert_final_newline = true +end_of_line = lf + +[*.{yml,yaml}] +indent_style = space +indent_size = 2 diff --git a/.github/.templateMarker b/.github/.templateMarker new file mode 100644 index 0000000..5e3a3e0 --- /dev/null +++ b/.github/.templateMarker @@ -0,0 +1 @@ +KOLANICH/python_project_boilerplate.py diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..89ff339 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,8 @@ +version: 2 +updates: + - package-ecosystem: "pip" + directory: "/" + schedule: + interval: "daily" + allow: + - dependency-type: "all" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..7fe33b3 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,15 @@ +name: CI +on: + push: + branches: [master] + pull_request: + branches: [master] + +jobs: + build: + runs-on: ubuntu-22.04 + steps: + - name: typical python workflow + uses: KOLANICH-GHActions/typical-python-workflow@master + with: + github_token: ${{ secrets.GITHUB_TOKEN }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ae045fa --- /dev/null +++ b/.gitignore @@ -0,0 +1,12 @@ +graph.dot + +__pycache__ +*.py[co] +/*.egg-info +*.srctrlbm +*.srctrldb +build +dist +.eggs +monkeytype.sqlite3 +/.ipynb_checkpoints diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml new file mode 100644 index 0000000..70a2f29 --- /dev/null +++ b/.gitlab-ci.yml @@ -0,0 +1,13 @@ +image: registry.gitlab.com/kolanich-subgroups/docker-images/fixed_python:latest + +variables: + DOCKER_DRIVER: overlay2 + SAST_ANALYZER_IMAGE_TAG: latest + SAST_DISABLE_DIND: "true" + SAST_CONFIDENCE_LEVEL: 5 + CODECLIMATE_VERSION: latest + +include: + - template: SAST.gitlab-ci.yml + - template: Code-Quality.gitlab-ci.yml + - template: License-Management.gitlab-ci.yml diff --git a/Code_Of_Conduct.md b/Code_Of_Conduct.md new file mode 100644 index 0000000..bcaa2bf --- /dev/null +++ b/Code_Of_Conduct.md @@ -0,0 +1 @@ +No codes of conduct! \ No newline at end of file diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..20f0fa8 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,4 @@ +include UNLICENSE +include *.md +include tests +include .editorconfig diff --git a/ReadMe.md b/ReadMe.md new file mode 100644 index 0000000..585f37d --- /dev/null +++ b/ReadMe.md @@ -0,0 +1,205 @@ +timberAllocation.py [![Unlicensed work](https://raw.githubusercontent.com/unlicense/unlicense.org/master/static/favicon.png)](https://unlicense.org/) +==================== +~~[wheel (GitLab)](https://gitlab.com/KOLANICH-research/timberAllocation.py/-/jobs/artifacts/master/raw/dist/timberAllocation-0.CI-py3-none-any.whl?job=build)~~ +[wheel (GHA via `nightly.link`)](https://nightly.link/KOLANICH-research/timberAllocation.py/workflows/CI/master/timberAllocation-0.CI-py3-none-any.whl) +~~![GitLab Build Status](https://gitlab.com/KOLANICH-research/timberAllocation.py/badges/master/pipeline.svg)~~ +~~![GitLab Coverage](https://gitlab.com/KOLANICH-research/timberAllocation.py/badges/master/coverage.svg)~~ +[![GitHub Actions](https://github.com/KOLANICH-research/timberAllocation.py/workflows/CI/badge.svg)](https://github.com/KOLANICH-research/timberAllocation.py/actions/) +[![Libraries.io Status](https://img.shields.io/librariesio/github/KOLANICH-research/timberAllocation.py.svg)](https://libraries.io/github/KOLANICH-research/timberAllocation.py) + +A not very working project. And completely unfinished paper for now full of bullshit. + +Sometimes one needs to make something not very loaded from timber, but has some wastes spared from previous construction. + +Assumptions: +* You have catalogued all the pieces. You have a TSV file, each row of which contains an id of a log and LENGTHS OF ITS PIECES (NOT full length, but lengths of all pieces sequentially, if the log is already made of multiple pieces and can be disassembled into them without cutting), sequentially. You can find a physical log given its id. For example ids are written on logs in large font with a permanent marker or paint (pen, nail scribing and pencil proven to be unreliable). +* You can cut any log into 2 pieces. +* Your construction tolerates pieces that are not pieces, but multiple logs stacked and fixed with nails, boards, metal plates and similar shit. +* All the pieces have the same section. + +This set of scripts aims to help in this situation. It tries to compute the optimal points to cut the existing timber. Though doesn't really help a lot: the results are suboptimal, inhuman (all the results were discarded by the human expert on the basis "I don't want to join these pieces") and practically not very usable. + +Optimality criteria (in an ideal case): +* The least of wastes is produced. A waste is an unused piece of a log that has been cut. +* The count of logs cut is minimal. +* The count of joins is minimal. + + +This is similar to both stock cutting problem (in which we cut pieces of identical size into pieces of different sizes in order to minimize waste and joins are not allowed) and ... problem (in which we join predefined different pieces to get ieces of the same size), both of which are known to be NP-hard. + + +Problem initial formalization as a search in a graph +-------------------------------------------------- + +We can formulate the problem as a search in the graph of states, where each state is: +* a set of available pieces of timber; +* slots. A slot is an amount of timber yet to be allocated. All the timber in a slot is joined into the same target piece. +* A `waste` variable. + +When there are no empty slots, we have found a solution. + +There are 2 possible kinds of actions in each state: +* `Cut` a piece into 2 pieces. Removes a piece from the set and adds 2 pieces which sum of lengths is equal to the length of the cut piece. Generates `waste` in amount of a cut piece. +* `Join` a piece and fill a slot with it. Removes a piece and decreases amount in the slot by it. Consumes `waste` in amount of a consumed piece. + +The initial state is when there is no waste and the pieces are the initial pieces. +The target state is when there is no empty slots. + +If all the cut pieces consumed there is no waste. + +This formulation of the task is useless because it is impossible to search in the graph for pieces > 3. But it gives us important insights we can use to reformulate the task. + +The insights: +* The problem always have a solution, if there is enough timber provided to occupy all the slots. We can always cut the logs into very (probably infinitily) small equal pieces and then reassemble new logs from them. +* All the pieces resulting from cuts persist untill consumed by joins. +* Thought we can cut any amount, we only need to cut the amounts which are differrences between the logs and the slots, because cuts that don't make us closer to filling a slot are useless. This allows us to search in a discrete space with a finite branching factor in each node! +* `Cut` and `Join` actions are MOSTLY (but not wholly) permuttible (when we say `permuttible` here, we mean that `a2(a1(s)) === a1(a2(s))`) + * `Join` can always follow a `Cut` producing a taken piece, but not opposite. + * `Cut` actions are always permuttible. + * `Join` actions can always be deferred to the end. So we ideally can put all the cuts to the beginning and joins to the end. But doing so is harmful in this formulation because a `Join` eliminates a piece and possibly a slot, reducing branching factor, but a `Cut` adds a piece increasing branching factor. +* Though our target solution requres the software to tell us the ids of pieces to cut and join, we cannot consider pieces as distinguishable in states because it creates redundant states. +* The most importantly, `Join` is an inverse of `Cut` in some sense. The algo can cut a piece and then join the cut pieces into a single one, essentially doing nothing except producing redundant cuts. + +All of these produce too many redundant paths, which are problematic to avoid, even if one avoids them, he misses the shortcuts they provide. But the fact that `Cut` and `Join` are inverses of each other and that we can defer `Join`s means we can do a bidirectional search. And in formulation for bidirectional search the problem is much easier! + +Problem reformalization 1 +------------------------- + +Each state is a set of pieces of timber. + +There are 2 initial states: +* All the pieces are the initial pieces. +* All the pieces are the target pieces and an additional piece for `unallocated`, if there may be any `unallocated` timber. + +There are only 1 possible kind of actions in each state: +* `Cut` a piece into 2 pieces. Removes a piece from the set and adds 2 pieces which sum of lengths is equal to the length of the cut piece. + +We run a bidirectional search from these 2 states, and somewhen they meet in the same state. They must meet - as we have seen the solution always exist. + +The insights: +* All the pieces resulting from cuts persist until consumed by joins. +* Let `M` is the middle state the searches meet in. It can contain either completely unknown pieces, or it can contain some pieces the same for either of initial states. +* If it contains some pieces from the initial states, since the pieces are indistinguishable and cuts are permuttible, we can claim they are the ones from the initial states. And we should do that, it will not make any cuts and any waste, and our goal is to minimize cuts and waste. +* If an unknown piece is longer than a piece in an initial state, it can be cut into the piece in the initial state and a remainder. Since the cuts are permuttible, this cut can be made from the initial state itself. Since we can only cut and the cut pieces persist unless further cut, this cut must be eventually made, sow we should to it now. The only variability here is the choice of the piece to cut. But the cut has to be done. +* So the optimal middle state minimizing cuts contains each piece from either one of initial states or both, if it is possible. + + +Problem reformalization 2 +------------------------- +Each state is 2 sets of pieces of timber + 2 sets of cuts. Pieces of timber sets are the difference between the middle state and the frontier states of the bidirectional search. Sets of cuts record the path. + +Initial state: sets of cuts are empty, sets of pieces are are initialized as in bidirectional search formalization. +Target state: sets of pieces are empty, sets of cuts contains the solution. + +There are 2 possible kind of actions in each state: +* `Annihilate` 2 pieces of equal length from the opposite sets, removing them from our consideration - they persist untill the middle state of the bidirectional search. +* `Cut` a piece of length `l0` into 2 pieces of lengths `l1` and `l2`, such as `l1 + l2 = l0` and `l1 >= l2`. Removes the piece `l0` from the set and adds 2 pieces `l1` and `l2`. + +Insigths: +* On each step sum of lengths of pieces in the sets is the same. +* `Annihilation`s are always permuttible with each other. +* `Cut`s are still permuttible and still can be postponed. Each `cut` can generate waste (if it is possible to have waste) and by itself is undesireable. +* `Cut`s can only contract lengths of individual pieces. +* Total length in opposite sets is equal and can only contract. So any `annihilation` brings us closer to the solution. +* Each `annihilation` also removes `cut` branches that can generate waste. +* Let `uncut` be amount of wood kept uncut, more precisely a piece is uncut if it is the same in the initial state and the final state. Then `waste = unallocated - uncut`. So to minimize waste we must maximize `uncut`. +* `Annihilate` nodes are idempotent and can always be applied. That allows to apply all annihilation operations possible in the curent state in a single node `AnnihilateAll`. +* If it is possible to make an `annihilation` of pieces or cutting a piece other than an annihilated one, this `cutting` and `annihilation` can be pemutted. It is though useful to do annihilation first. + + +Theorem: Only the steps having any of initial pieces determine the amount of waste. Whatever we do on the other steps, we cannot decrease amount of waste. +Proof: + While we have initial pieces, we can increase amount of waste by cutting them. While we don't, we cannot. +Corrolary: + As soon as we either annihilate or cut all the initial pieces we have fixed the amount of waste. + +Let's call a cut `necessary` if we have not solved the problem and we have no pieces to annihilate, so the only action we can do is to cut. +Let's call a cut `final` if after an annihilation has followed this cut, there will be no pieces left. + +Theorem: Each necessary nonfinal `cut` in one side can either increment count of `necessary` cuts either by 1 or by 2 +Proof: +We have 3 options: +1. If a necessary cut produces a piece of the same size that the one of the existing pieces, we would be able to annihilate them. The resulting piece will either be annihilated somewhen because there must be a solution (this means the cut will has to be done from the opposite side of the search to match it), or will be further cut. +2. If a necessary cut produces no pieces that can be annihilated, we have to cut 2 other pieces to match the sizes of the resulting pieces, as in 1. It is not possible to already have the pieces matching the sizes of any of the resulting pieces, otherwise this either would have been option 1 or it . +3. The final cut must produce the pieces that will be annihilated all, so it increments the count of further needed cuts by 0. + +So it makes sense to combine `Annihilate`s to `AnnihilateAll` and `AnnihilateAll` to the previous `Cut`, having only `CutAndAnnihilateAll` node. + + +Problem reformalization 3 +------------------------- +The same as in 2, but + +There are only 1 possible kind of actions in each state: +* `AnnihilateAndCut` a piece of length `l0` into 2 pieces of lengths `l1` and `l2`, such as `l1 + l2 = l0` and `l1 >= l2` and either `l1` or `l2` piece is present in the opposite set, so is guaranteedly annihilated. Further we call this operation a `cut`. Removes at least one piece from the opposite set and can also remove pieces from own set. + +Insights: +* upper bound on sum of counts of cuts from the both sides is ` - 2`. Proof: `(count of pieces from initial side - 1) + (count of pieces from target side - 1)` nonfinal cuts consuming 1 piece from either side and 1 final cut consuming both pieces at worst. +* We can also predict the upper bound on count of pieces in the middle state in bidi search: `(count of pieces from initial side - 1) + (count of pieces from target side - 1)` nonfinal cuts add 1 piece to the middle state, and the final cut adds 2, so `count of pieces from initial side + count of pieces from target side - 1`. + +We have a computation vs cuts vs waste tradeoff. Selecting the middle state to be: +* near the initial state minimizes count of cuts and possibly waste (in presence of redundant material) at expense of greater count of joins and computation. This can be achieved by generating cuts on initial side only if it is not possible to generate cuts on target side. +* near the target state minimizes count of joins at expense of greater count of cuts and waste (in presence of redundant material). This can be achieved by generating cuts on target side only if it is not possible to generate cuts on initial side. +* the one balancing cuts and joins from both sides. +* the one minimizing search cost at expense of both cuts and joins relatively to the minimal cuts/joins states. If we had a constant branching factor, this would have been the previous one, but count of children in each node is not easily predictable since cuts are subjet to constraints based on pieces sizes and they change. + +Depending on the problem, sometimes these variants of the middle state are the same, sometimes they are different. + +But we are still plagued with redundant paths because of permuttability, but we cannot prioritize cuts on size because some cuts possible in the state are mutually excluding. So we need to group cuts into the groups + + +A simple greedy nonoptimal non-search algo +------------------------------------------ + +A little step away from the solving of the problem with graph search is a simple algo derived from the reformalization 3. + +1. We have 2 arrays of parts as described. For example, `parts_I` is array of initial parts, `parts_T` is the array of target parts. Each part has an id, length and a direction from the set `{"I", "T"}` which means the direction of the search and matches the direction of the array. +2. We concatenate them: `parts_I || parts_T`. +3. We sort this array `O(n*log(n))`. +4. We crawl that array for `O(n)` in 2 passes, annihilating adjacent elements of opposite directions. + a. In the first pass we only annihilate elements that are equal. They must be adjacent since the arrays are sorted. + a. In the second one we partially annihilate elements that are not equal, leaving only the larger ones. + +5. We repeat this until the array is not empty. +6. Then we recover the solution from the set of cuts. + +Though it is still desireable to use a search algo since + a. it can find better solutions + b. can be reused for the modified problem, where pieces have not only lengths, but also fixed widths, and we want to make multiple rectangular frames of them having known dimensions. + + +Vector factorization +-------------------- + +Let's do some linear algebra. + +Let we have a row-vector of positive numbers `x`. If we can factor it into a product `b . U === x`, where `b` is a row-vector of minimal size, that will be our basis, `U` is the "upper-triangular-like" (mostly upper triangular, but with holes in the upper triangle and elements in the lower triangle when necessary) matrix, which elements are either 1 or 0, then we can probably have a common representation between the numbers we can operate in. + +The first insight, that since the middle state has at most `(count of pieces from initial side + count of pieces from target side) - 1` pieces, `b` is guaranteedly smaller than `x` at least by 1 if 'x' is made of numbers matching both sides of the bidi search. But in general case it is not the case and for vectors of `n` elements in the worst case we will have basis of `n` elements, if all the numbers are lineary independent. In the case of cuts we have `n-1` at worst since they are guaranteedly lineary dependent. + +Factoring is straightforward and is just a mixture of Gram-Schmidt orthogonalization, Gaussian elimination and Euclidus algorithm. + +1. The largest number should have all its components present. +2. The second largest should have all its components except the one differing it to the largest number present. + +So we repeat the following loop until the residues goes to 0. `i` is from `-1` to `-n`: +1. We find the largest unique number `L` and the second largest number `l`. If there is only one number, we set `L` to it and `l` to `0`; +2. We differ them `e[i] = L - l`. This will be our basis number. We assign `0` to the second largest and `1` to the largest for `i`th component. +3. We find all the elements in the vector matching the second largest element and just set their `i`th component to `0` by definition. +4. For other elements greater than the current basis number we subtract it and set `1` for the component, otherwise we set `0`. + +Some basis numbers can be duplicated. This is because we allow components to be only 1 and 0 - either taken by part or not. + +We can also apply the same procedure to the resulting basis again and multiply the matrices ... and again and again. This will give us a smaller basis, but the matrices will be non-binary. Though visualization of the matrix looks like it may be helpful for unsupervised learning. Needs some additional research. + + +Problem reformalization 4 +------------------------- + +OK, now instead of raw lengths of the pieces we have binary vectors. Taking one of the components from the piece doesn't exclude other components. So we can deal with permuttibility by enumerating all the combinations of components that can be taken. This increases counts of children in one node but completely eliminates permutations between nodes of the same part - once a part is cut, we can never cut it again anymore. Permutations between nodes of different parts can be eliminated too - now we always know which chunks we are allowed to cut from a part, independent from the history and side branches. + +** Fucking shit - it doesn't help. Probably a bug in the algo. BrFS (breadth-first search) that has worked for the small example has become intractable. GBeFS (greedy best-first search) doesn't achieve the solution. Even if we don't do combinations and mirror the old behavior with cutting part-by-part.** + +Solution strategy +----------------- + diff --git a/UNLICENSE b/UNLICENSE new file mode 100644 index 0000000..efb9808 --- /dev/null +++ b/UNLICENSE @@ -0,0 +1,24 @@ +This is free and unencumbered software released into the public domain. + +Anyone is free to copy, modify, publish, use, compile, sell, or +distribute this software, either in source code form or as a compiled +binary, for any purpose, commercial or non-commercial, and by any +means. + +In jurisdictions that recognize copyright laws, the author or authors +of this software dedicate any and all copyright interest in the +software to the public domain. We make this dedication for the benefit +of the public at large and to the detriment of our heirs and +successors. We intend this dedication to be an overt act of +relinquishment in perpetuity of all present and future rights to this +software under copyright law. + +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 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. + +For more information, please refer to diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..f16229b --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,40 @@ +[build-system] +requires = ["setuptools>=61.2.0", "setuptools_scm[toml]>=3.4.3"] +build-backend = "setuptools.build_meta" + +[project] +name = "timberAllocation" +readme = "ReadMe.md" +description = "Just a tool trying to solve the task of optimal cutting existing pieces of timber into pieces that when combined give pieces of needed length" +authors = [{name = "KOLANICH"}] +classifiers = [ + "Development Status :: 4 - Beta", + "Environment :: Other Environment", + "Intended Audience :: Developers", + "License :: Public Domain", + "Operating System :: OS Independent", + "Programming Language :: Python", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3 :: Only", + "Topic :: Software Development :: Libraries :: Python Modules", +] +keywords = ["algorythms", "greedy", "divide-and-conquer", "dynamic programming", "optimization"] +license = {text = "Unlicense"} +requires-python = ">=3.4" +dynamic = ["version"] +dependencies = [ + "RichConsole", # @ git+https://gitlab.com/KOLANICH/RichConsole.py.git +] +[project.optional-dependencies] +optimization = ["autograd"] + +[project.urls] +Homepage = "https://github.com/KOLANICH-research/timberAllocation.py" + +[tool.setuptools] +zip-safe = true + +[tool.setuptools.packages.find] +include = ["timberAllocation", "timberAllocation.*"] + +[tool.setuptools_scm] diff --git a/tests/47,31->51,22,5.dot b/tests/47,31->51,22,5.dot new file mode 100644 index 0000000..2e649c2 --- /dev/null +++ b/tests/47,31->51,22,5.dot @@ -0,0 +1,19 @@ +digraph "" { + 47 [color=green]; + 31 [color=green]; + 47 [color=green]; + 51 [color=red]; + 22 [color=red]; + 5 [color=red]; + 47 -> {25, 22} [color=green]; + 47 -> {20, 27} [color=green]; + 31 -> {5, 26} [color=green]; + 31 -> {9, 22} [color=green]; + 9 -> {5, 4} [color=cyan]; + 51 -> {47, 4} [color=red]; + 25 -> {20, 5} [color=cyan]; + 51 -> {20, 31} [color=red]; + 27 -> {22, 5} [color=cyan]; + 26 -> {22, 4} [color=cyan]; + 31 -> {4, 27} [color=green]; +} diff --git a/tests/47,31->51,22,5.tsv b/tests/47,31->51,22,5.tsv new file mode 100644 index 0000000..cd0b3be --- /dev/null +++ b/tests/47,31->51,22,5.tsv @@ -0,0 +1,3 @@ +#{51:1, 22:1, 5:1} +1 47 +2 31 diff --git a/tests/toy.tsv b/tests/toy.tsv new file mode 100644 index 0000000..7f5b6c6 --- /dev/null +++ b/tests/toy.tsv @@ -0,0 +1,5 @@ +#{6.7: 1, 4.3: 1} +1 3.87 +2 3.17 +3 4 +4 3.5 diff --git a/tests/wood.tsv b/tests/wood.tsv new file mode 100644 index 0000000..f06efb8 --- /dev/null +++ b/tests/wood.tsv @@ -0,0 +1,42 @@ +#{230: 9, 110: 10} +14 242 +32 234 +36 141 89 +27 230 +2 224 +43 219 +3 209 +33 209 +15 208 +24 178 29 +31 178 29 +29 154.5 50 +19 170 32 +41 186 +10 184 +42 180 +16 164.5 +6 163 +40 160 +39 103 56 +17 153 +25 152 +23 150 +28 150 +12 147 +8 142.5 +34 139 +26 136 +30 104 +13 99 +35 89 +20 86 +7 84.5 +9 80 +11 78 +1 31 +5 27 +37 23 +22 20 +18 14 +21 6.5 diff --git a/tests/wood1.tsv b/tests/wood1.tsv new file mode 100644 index 0000000..8ce7522 --- /dev/null +++ b/tests/wood1.tsv @@ -0,0 +1,6 @@ +#{46: 1, 44: 1, 26:1, 27:1, 110: 1, 80: 1, 11: 1, 21: 1, } +36 89 +39 103 +34 139 +22 20 +18 14 diff --git a/timberAllocation/__init__.py b/timberAllocation/__init__.py new file mode 100644 index 0000000..0c52590 --- /dev/null +++ b/timberAllocation/__init__.py @@ -0,0 +1,381 @@ +from collections import defaultdict +from defaultlist import defaultlist +from functools import partial +import itertools +import math +import typing +import csv +import ast +from pathlib import Path +from copy import deepcopy +from decimal import Decimal +import numpy as np + +import RichConsole + +from .utils import TupleHashable, StructuredComparable + +#FloatType = Decimal +FloatType = float +NumT = typing.Union[float, int, Decimal, np.float64] +TSVReadT = typing.List[typing.Tuple[NumT]] + + +def readFromTSV(path: Path) -> typing.Tuple[typing.Mapping, TSVReadT]: + with path.open("rt", encoding="utf-8") as f: + lsIter = iter(f) + firstLine = next(lsIter) + if firstLine[0] != "#": + raise ValueError("First line must be a comment specifying the problem") + + firstLine = firstLine[1:].lstrip() + task = ast.literal_eval(firstLine) + + decimalPrecision = FloatType("0.0001") + + if FloatType is Decimal: + task = {FloatType(k).quantize(decimalPrecision): v for k, v in task.items()} + + ds = [] + for r in csv.reader(lsIter, dialect=csv.excel_tab): + rI = r[0] + try: + rI = int(rI) + except ValueError: + pass + + rNew = [rI] + for i in range(1, len(r)): + n = r[i] + if n: + try: + n = int(n) + except ValueError: + n = FloatType(n) + rNew.append(n) + + r = tuple(rNew) + ds.append(r) + return task, getParts(ds) + + +class PartId(TupleHashable): + __slots__ = ("_origId", "chunkId") + + def _coIter(self, other) -> typing.Iterable[typing.Tuple[typing.Any, typing.Any]]: + yield "origId", self.origId, other.origId + yield "chunkId", self.chunkId, other.chunkId + + def __init__(self, origId: typing.Any, chunkId: typing.Any = ()) -> None: + if not isinstance(origId, tuple): + origId = (origId,) + self._origId = origId + self.chunkId = chunkId + + @property + def origId(self) -> "PartId": + return PartId(self._origId) + + def toTuple(self) -> typing.Union[typing.Tuple]: + return self._origId + (self.chunkId if self.chunkId else ()) + + def __repr__(self) -> str: + return repr(self.toTuple()) if self.chunkId else (repr(self._origId) if len(self._origId) > 1 else repr(self._origId[0])) + + def __add__(self, other: typing.Tuple) -> "PartId": + res = self.__class__(self._origId, self.chunkId + other) + return res + + def __lt__(self, other: "PartId") -> bool: + return self.toTuple() < other.toTuple() + + def __gt__(self, other: "PartId") -> bool: + return self.toTuple() > other.toTuple() + + def __le__(self, other: "PartId") -> bool: + return self.toTuple() <= other.toTuple() + + def __ge__(self, other: "PartId") -> bool: + return self.toTuple() >= other.toTuple() + + +class IPart: + __slots__ = ("id",) + + def __init__(self, id: PartId) -> None: + self.id = id + + def __lt__(self, other: "Part") -> bool: + return self.len < other.len + + def __gt__(self, other: "Part"): + return self.len > other.len + + def __ge__(self, other: "Part"): + return self.len >= other.len + + def __le__(self, other: "Part"): + return self.len <= other.len + + def toTuple(self): + return (self.id,) + + def __repr__(self): + return self.__class__.__name__[0] + "<" + repr(self.id) + ", " + str(self.len) + ">" + + +class Part(IPart): + __slots__ = ("id", "len") + + def __init__(self, id: PartId, len: NumT) -> None: + super().__init__(id) + self.len = len + + def toTuple(self): + return super().toTuple() + (self.len,) + + def split(self, splitPoint: NumT) -> typing.Tuple["Part", "Part"]: + secondChunkLen = self.len - splitPoint + firstChunk = deepcopy(self) + firstChunk.id += (0,) + firstChunk.len = splitPoint + + secondChunk = deepcopy(self) + secondChunk.id += (1,) + secondChunk.len = secondChunkLen + + return (firstChunk, secondChunk) + + +class DirPart(Part): + __slots__ = ("dir",) + + def __init__(self, id: PartId, len: NumT, dir: int) -> None: + super().__init__(id, len) + self.dir = dir + + def __repr__(self): + return str((RichConsole.groups.Fore.lightRed if self.dir else RichConsole.groups.Fore.lightGreen)(super().__repr__())) + + +class VectorPart(IPart): + __slots__ = ("basis", "pieces") + + def __init__(self, id: PartId, piecesVector: "np.ndarray", basis: typing.Optional[np.ndarray] = None) -> None: + if not isinstance(piecesVector, np.ndarray): + raise ValueError() + super().__init__(id) + self.pieces = piecesVector + self.basis = basis + + def toTuple(self): + return super().toTuple() + (self.pieces,) + + def split(self, *chunks: "np.ndarray") -> typing.Iterable["VectorPart"]: + restChunkPieces = deepcopy(self.pieces) + for i, piecesToRetain in enumerate(chunks): + if not isinstance(piecesToRetain, np.ndarray): + raise ValueError(piecesToRetain) + + restChunkPieces = restChunkPieces & ~piecesToRetain + chunk = deepcopy(self) + chunk.id += (i,) + chunk.pieces = piecesToRetain + + yield chunk + + if restChunkPieces.any(): + chunk = deepcopy(self) + chunk.id += (i,) + chunk.pieces = restChunkPieces + yield chunk + + @property + def cutted(self) -> bool: + return self.id != self.id.origId + + @classmethod + def _subParts(cls, pieces: np.ndarray) -> np.ndarray: + print(pieces) + return np.identity(pieces.shape[0])[pieces] + + def subParts(self): + return self.__class__._subParts(self.pieces) + + @classmethod + def _getCuttings(cls, pieces: np.ndarray) -> typing.Iterator[typing.Iterable[np.ndarray]]: + subParts = cls._subParts(pieces) + + for i in range(1, len(subParts) // 2 - 1): # the rest are the complements + for firstPieceParts in itertools.combinations(subParts, i): + firstPieceVec = np.sum(firstPieceParts, axis=0).astype(bool) + secondPieceVec = pieces & ~firstPieceVec + yield (firstPieceVec, secondPieceVec) + for secondPieceCuttings in cls._getCuttings(secondPieceVec): + yield (firstPieceVec,) + secondPieceCuttings + + def getCuttings(self) -> typing.Iterator[typing.Iterable[np.ndarray]]: + return self.__class__._getCuttings(self.pieces) + #for firstPieceVec in self.__class__._getCuttings(self.pieces): + # yield self.split(firstPieceVec) + + @property + def len(self) -> FloatType: + return self.basis @ self.pieces + + +class DirVectorPart(VectorPart): + __slots__ = ("dir",) + + def __init__(self, id: PartId, piecesVector: "np.ndarray", basis: np.ndarray, dir: int) -> None: + super().__init__(id, piecesVector, basis) + self.dir = dir + + def __repr__(self): + return str((RichConsole.groups.Fore.lightRed if self.dir else RichConsole.groups.Fore.lightGreen)(super().__repr__())) + + +def getParts(ds: TSVReadT) -> typing.Dict[PartId, NumT]: + parts = {} + for r in ds: + rId = r[0] + r = r[1:] + if len(r) == 1: + # parts[rId] = r[0] + parts[PartId(rId,)] = r[0] + else: + for i, rp in enumerate(r): + parts[PartId((rId, i))] = rp + return parts + + +class Slot(TupleHashable): + __slots__ = ("origL", "inSizeId", "l") + + def __init__(self, inSizeId: int, l: NumT) -> None: + self.origL = l + self.inSizeId = inSizeId + self.l = l + + def toTuple(self) -> typing.Tuple[NumT, int, NumT]: + return (self.origL, self.inSizeId, self.l) + + def __repr__(self) -> str: + return repr(self.toTuple()) + + +def getSlots(task: typing.Dict[int, int]) -> typing.Dict[int, Slot]: + slots = {} + j = 0 + for slotL, count in task.items(): + for i in range(count): + slots[j] = Slot(i, slotL) + j += 1 + return slots + + +class HashableDefaultDict(TupleHashable, defaultdict): + __slots__ = () + + def _toTuple(self) -> typing.Iterator[typing.Tuple]: + for k, v in self.items(): + yield (k, v.toTuple()) + + def toTuple(self) -> typing.Tuple: + return tuple(self._toTuple()) + + +class HashableDefaultList(TupleHashable, defaultlist): + __slots__ = () + + def _toTuple(self) -> typing.Iterator[typing.Tuple]: + for v in self: + yield v.toTuple() + + def toTuple(self) -> typing.Tuple: + return tuple(self._toTuple()) + + +class HashableList(TupleHashable, list): + __slots__ = () + + def toTuple(self) -> typing.Tuple: + return tuple(self) + + +def initSolution() -> typing.DefaultDict[NumT, typing.List[NumT]]: + return HashableDefaultDict(partial(HashableDefaultList, HashableList)) + + +def greedySolution(task: typing.Dict[NumT, int], parts: typing.Dict[PartId, NumT]) -> typing.Iterator[HashableDefaultDict]: + from .greedy import getCandidates, iteration + + opParts = deepcopy(parts) + slots = getSlots(task) + solution = initSolution() + + cumDelta = 0 + + while slots: + cumDelta = iteration(cumDelta, slots, slots, opParts, opParts, solution) + + yield solution + + +def greedyGenerationsSolution(task, parts): + from .greedy import getCandidates, iteration + + opParts = deepcopy(parts) + slots = getSlots(task) + solution = initSolution() + + cumDelta = 0 + + garbageSlots = {} + garbageParts = {} + + while slots: + cumDelta = iteration(cumDelta, slots, garbageSlots, opParts, garbageParts) + + while garbageSlots: + cumDelta = iteration(cumDelta, garbageSlots, slots, garbageParts, garbageParts) + + yield solution + + +def optimizeSolutionUsingGradientDescent(alpha, L, T, parts): + from .optimization import opt, alphaToSolution + from matplotlib import pyplot as plt + + plt.matshow(alpha) + plt.colorbar() + plt.show() + + res = opt(L, alpha, T) + + plt.matshow(res) + plt.colorbar() + plt.show() + + solution = alphaToSolution(L, res, T, parts) + yield solution + + +def gradientDescentSolution(task, parts): + from .optimization import genAlpha, genL, genTarg + + T = genTarg(getSlots(task)) + L = genL(parts) + alpha = genAlpha(L, T) + return optimizeSolutionUsingGradientDescent(alpha, L, T, parts) + + +def hybridSolution(task, parts): + from .optimization import solutionToAlphaMatrix, genL, genTarg + + sol = next(iter(greedySolution(task, parts))) + + T = genTarg(getSlots(task)) + L = genL(parts) + alpha = solutionToAlphaMatrix(task, parts, sol) + return optimizeSolutionUsingGradientDescent(alpha, L, T, parts) diff --git a/timberAllocation/__main__.py b/timberAllocation/__main__.py new file mode 100644 index 0000000..c7f9cbf --- /dev/null +++ b/timberAllocation/__main__.py @@ -0,0 +1,59 @@ +import itertools +import json +import math +from ast import literal_eval +from copy import deepcopy +from decimal import Decimal +from pathlib import Path +from pprint import pformat + +import numpy as np +from plumbum import cli + +from . import getParts, gradientDescentSolution, readFromTSV +from .sortSolution import sortSolution +from .stateGraph import Algo +from .vis import generateColors, visualizeAll + + +def algoSolution(task, parts): + a = Algo(task, parts) + solutionNodes = a() + for solutionNode in solutionNodes: + path = a.getPath(solutionNode) + print(path) + yield a.getSolutionFromPath(path) + + +methods = { + "sort": sortSolution, + "graph": algoSolution, +} + +from .z3Solution import z3Solution +methods["z3"] = z3Solution + + +thisDir = Path(__file__).absolute().parent +moduleParentDir = thisDir.parent +testsDir = moduleParentDir / "tests" + + +class AllocationCLI(cli.Application): + #method = cli.SwitchAttr(["-m", "--method"], argtype=cli.Set(*methods), argname="method", help="Method of solution.", default="sort") + method = cli.SwitchAttr(["-m", "--method"], argtype=cli.Set(*methods), argname="method", help="Method of solution.", default="z3") + + #def main(self, input=testsDir / "47,31->51,22,5.tsv"): + #def main(self, input=testsDir / "wood1.tsv"): + #def main(self, input=testsDir / "wood.tsv"): + def main(self, input=testsDir / "toy.tsv"): + task, parts = readFromTSV(Path(input)) + partsColors = generateColors(parts) + + for solution in methods[self.method](task, parts): + print("\n") + visualizeAll(solution, parts, partsColors) + + +if __name__ == "__main__": + AllocationCLI.run() diff --git a/timberAllocation/middleRepr.py b/timberAllocation/middleRepr.py new file mode 100644 index 0000000..2454679 --- /dev/null +++ b/timberAllocation/middleRepr.py @@ -0,0 +1,78 @@ +import typing + +import numpy as np + + +def getResultVector(initialLengths, finalLengths): + return np.array(tuple(initialLengths) + tuple(finalLengths)) + + +def getMiddleStateMaxSize(initialLengths, finalLengths): + return len(initialLengths) + len(finalLengths) - 1 + + +def vectorFactorization(vec: np.ndarray) -> typing.Tuple[np.ndarray, np.ndarray]: + """I have probably reinvented a wheel. I have searched the Internet and haven't found this kind of factorization. + Factors a vector vec into the product of an "upper-triangular-like" (mostly upper triangular, but with holes and elements in the lower triangle when necessary) made of 1es matrix `U` and a basis row-vector `b`, so b @ U === vec. + """ + + vec = np.array(vec) + + dtype = vec[0].__class__ # so works with Decimal too. In this case arrays dtype is `object`. + + basis = np.full(vec.shape[0], dtype(np.nan)) + matrix = np.zeros((basis.shape[0], vec.shape[0]), dtype=bool) + + i = -1 + + basisSize = 0 + while (vec > 0).any(): + remainingSizes = sorted(set(vec)) + if len(remainingSizes) >= 2: + secondLargest, largest = remainingSizes[-2:] + basisVec = largest - secondLargest + else: # 1 size only + basisVec = remainingSizes[0] + secondLargest = 0 + + basis[i] = basisVec + + for j, s in enumerate(vec): + if s == secondLargest: + matrix[matrix.shape[0] + i, j] = 0 + else: + if s >= basisVec: + matrix[matrix.shape[0] + i, j] = 1 + vec[j] -= basisVec + else: + matrix[matrix.shape[0] + i, j] = 0 + + i -= 1 + basisSize += 1 + + return basis[-basisSize:], matrix[-basisSize:, :] + + +def minimalReprToGraph(shared, mat, initialLengths, finalLengths): + import networkx + + g = networkx.DiGraph() + for n in initialLengths: + if isinstance(n, float) and n.is_integer(): + n = int(n) + g.add_node(n, color="green") + for n in finalLengths: + if isinstance(n, float) and n.is_integer(): + n = int(n) + g.add_node(n, color="red") + + for i, l in enumerate(initialLengths + finalLengths): + if isinstance(l, float) and l.is_integer(): + l = int(l) + for j, sL in enumerate(shared): + sL = float(sL) + if sL.is_integer(): + sL = int(sL) + if sL != l and mat[j, i]: + g.add_edge(l, sL) + return g diff --git a/timberAllocation/sortSolution.py b/timberAllocation/sortSolution.py new file mode 100644 index 0000000..4200e8f --- /dev/null +++ b/timberAllocation/sortSolution.py @@ -0,0 +1,115 @@ +import sys +import typing +from decimal import Decimal +from functools import partial +import heapq + +import more_itertools + +from . import DirPart, HashableDefaultDict, NumT, Part, PartId, Slot, getParts, getSlots, initSolution + +__all__ = ("sortSolution",) + + +def initialize(task: typing.Dict[Decimal, int], parts: typing.Dict[PartId, NumT]): + firstTotal = Decimal(0) + + first = [] + for partId, partLen in parts.items(): + first.append(DirPart(partId, partLen, 0)) + firstTotal += partLen + + secondTotal = Decimal(0) + + second = [] + for i, (partLen, count) in enumerate(task.items()): + for j in range(count): + second.append(DirPart(PartId((partLen, j)), partLen, 1)) + secondTotal += partLen + + unusedLen = firstTotal - secondTotal + + if unusedLen: + if unusedLen > 0: + second.append(DirPart(PartId("U"), unusedLen, 1)) + else: + raise ValueError("The problem is unsolvable: length in the left is greater than the length in the right", unusedLen) + + return first, second + + +def annihilateAdjacent(joint, solution, allowSplitting): + secondConsumed = False + newJoint = [] + + haveStalled = True + + f = heapq.heappop(joint) + while joint: + s = heapq.heappop(joint) + if secondConsumed: + secondConsumed = False + else: + #print("annihilateInequalPairUnparam", f, s, f.dir != s.dir, f.len - s.len) + if f.dir == s.dir: + heapq.heappush(newJoint, f) + secondConsumed = False + else: + i, t = (f, s) if f.dir == 0 else (s, f) + + if abs(f.len - s.len) <= sys.float_info.epsilon: + l = f.len + haveStalled = False + secondConsumed = True + else: + if t.len > i.len: + l = i.len + toAppend = t.split(l)[1] + haveStalled = False + heapq.heappush(newJoint, toAppend) + if t.id._origId[0] != "U": + solution[t.id._origId[0]][t.id._origId[1]].append((i.id.origId, l)) + secondConsumed = True + else: # i.len > t.len + if allowSplitting: + l = t.len + toAppend = i.split(l)[1] + haveStalled = False + heapq.heappush(newJoint, toAppend) + if t.id._origId[0] != "U": + solution[t.id._origId[0]][t.id._origId[1]].append((i.id.origId, l)) + secondConsumed = True + else: + heapq.heappush(newJoint, f) + secondConsumed = False + + if secondConsumed: + newJoint.extend(joint) + return newJoint, False + f = s + + if not secondConsumed: + heapq.heappush(newJoint, s) + return newJoint, haveStalled + + +def solve(joint): + solution = initSolution() + + print(joint) + + allowCutting = False + + while joint: + joint, allowCutting = annihilateAdjacent(joint, solution, allowCutting) + joint = sorted(joint) + print(joint) + + yield solution + + +def sortSolution(task: typing.Dict[Decimal, int], parts: typing.Dict[PartId, NumT]) -> typing.Iterator[HashableDefaultDict]: + first, second = initialize(task, parts) + joint = first + second + joint = sorted(joint) + return solve(joint) diff --git a/timberAllocation/stateGraph.py b/timberAllocation/stateGraph.py new file mode 100644 index 0000000..3508a82 --- /dev/null +++ b/timberAllocation/stateGraph.py @@ -0,0 +1,898 @@ +import typing +from abc import ABC, abstractmethod +from collections import defaultdict, deque +from copy import deepcopy +from decimal import Decimal +from heapq import heappop, heappush +from math import inf +from weakref import ref + +import numpy as np + +import networkx as nx +import networkx.drawing.nx_agraph +import RichConsole + +from . import DirPart, DirVectorPart, FloatType, NumT, Part, PartId, Slot, getParts, getSlots, initSolution +from .middleRepr import vectorFactorization +from .utils import RichReprable, RichStringableAndReprable, StructuredComparable, TupleHashable + + +class Metrics(RichStringableAndReprable, TupleHashable, StructuredComparable): + __slots__ = () + + def _coIter(self, other: "Metrics") -> typing.Iterable[typing.Tuple[typing.Any, typing.Any]]: + for k in self.__class__.__slots__: + sv = getattr(self, k) + ov = getattr(other, k) + yield (k, sv, ov) + + def intize(self) -> None: + for k in self.__class__.__slots__: + v = getattr(self, k) + if isinstance(v, float) and v.is_integer(): + setattr(self, k, int(v)) + + def toTuple(self) -> typing.Tuple[typing.Type["State"], typing.Any, NumT]: + return tuple(getattr(self, k) for k in self.__class__.__slots__) + + +class PathMetrics(Metrics): + __slots__ = ("_piecesInInitial", "_piecesInTarget", "lengthToAllocate", "chunksToAnnihilate") + + def __init__(self) -> None: + self._piecesInInitial = None + self._piecesInTarget = None + self.lengthToAllocate = None + self.chunksToAnnihilate = None + + @property + def piecesInInitial(self): + return self._piecesInInitial + + @piecesInInitial.setter + def piecesInInitial(self, v): + if v < 0: + raise ValueError("piecesInInitial", v) + self._piecesInInitial = v + + @property + def piecesInTarget(self): + return self._piecesInTarget + + @piecesInTarget.setter + def piecesInTarget(self, v): + if v < 0: + raise ValueError("piecesInTarget", v) + self._piecesInTarget = v + + def __rich_str__(self) -> RichConsole.RichStr: + return "P(" + RichConsole.groups.Fore.lightGreen(str(self.piecesInInitial)) + ", " + RichConsole.groups.Fore.lightRed(str(self.piecesInTarget)) + "), L" + RichConsole.groups.Fore.lightBlue(str(self.lengthToAllocate)) + + +class QualityMetrics(Metrics): + __slots__ = ("cutsCount", "joinsCount") + + def __init__(self) -> None: + self.cutsCount = None + self.joinsCount = None + + @property + def cutsAndJoinsCount(self) -> int: + return self.cutsCount + self.joinsCount + + def _coIter(self, other: "QualityMetrics") -> typing.Iterable[typing.Tuple[typing.Any, typing.Any]]: + yield ("cutsAndJoinsCount", self.cutsAndJoinsCount, other.cutsAndJoinsCount) + + def __rich_str__(self) -> RichConsole.RichStr: + return "C(" + RichConsole.groups.Fore.lightCyan(str(self.cutsCount)) + "," + RichConsole.groups.Back.lightYellow(RichConsole.groups.Fore.black(str(self.joinsCount))) + ")" + + +class MetricsSet(RichStringableAndReprable, TupleHashable): + __slots__ = ("path", "quality") + + def __init__(self) -> None: + self.path = PathMetrics() + self.quality = QualityMetrics() + + def toTuple(self) -> typing.Tuple[typing.Type["State"], typing.Any, NumT]: + return tuple(getattr(self, k) for k in self.__class__.__slots__) + + def __rich_str__(self) -> RichConsole.RichStr: + return self.path.__rich_str__() + ", " + self.quality.__rich_str__() + + def intize(self) -> None: + for k in self.__class__.__slots__: + getattr(self, k).intize() + + def __eq__(self, other: "MetricsSet") -> bool: + return self.path == other.path + + def __lt__(self, other: "MetricsSet") -> bool: + return self.path < other.path and self.quality <= other.quality + + def __gt__(self, other: "MetricsSet") -> bool: + return self.path > other.path or self.quality > other.quality + + def __le__(self, other: "MetricsSet") -> bool: + return self.path < other.path or (self.path == other.path and self.quality <= other.quality) + + def __ge__(self, other: "MetricsSet") -> bool: + return self.path > other.path or (self.path == other.path and self.quality >= other.quality) + + +class Metered: + __slots__ = ("metrics",) + + def __init__(self) -> None: + self.metrics = None + + def __lt__(self, other: "Metered") -> bool: + return self.metrics < other.metrics + + def __gt__(self, other: "Metered") -> bool: + return self.metrics > other.metrics + + def __ge__(self, other: "Metered") -> bool: + return self.metrics >= other.metrics + + def __le__(self, other: "Metered") -> bool: + return self.metrics <= other.metrics + + +class PartsStorage(TupleHashable): + """Stores only pats ids. Takes less memory but it is not always possible to use it""" + + __slots__ = ("sizes",) + + def __init__(self) -> None: + self.sizes = {} + + def add(self, part: Part) -> None: + s = self.sizes.get(part.len, None) + if s is None: + self.sizes[part.len] = s = set() + + s |= {part.id} + + def remove(self, part: Part) -> None: + self.sizes[part.len] -= {part.id} + if not self.sizes[part.len]: + del self.sizes[part.len] + + @classmethod + def _getFirstIdUncut(cls, idsOfTheLen): + for iD in idsOfTheLen: + if idsOfTheLen == idsOfTheLen.origId: + return iD + + def get(self, size: NumT, uncut=False): + idsOfTheLen = self.sizes[size] + if uncut: + iD = self.__class__._getFirstIdUncut(idsOfTheLen) + else: + iD = next(iter(idsOfTheLen)) + return Part(iD, size) + + def pop(self, size: NumT, uncut: bool = False) -> Part: + p = self.get(size, uncut) + self.remove(p) + return p + + def resize(self, part, newSize): + self.remove(part) + self.add(Part(part.id, newSize)) + + def parts(self) -> typing.Iterable[Part]: + for size, sizeGroup in self.sizes.items(): + for iD in sizeGroup: + yield Part(iD, size) + + def toTuple(self) -> typing.Tuple: + res = [] + for size, sizeGroup in self.sizes.items(): + res.append((size, tuple(sorted(sizeGroup)))) + return tuple(res) + + def __repr__(self) -> str: + els = [] + for el in sorted(self.sizes): + count = len(self.sizes[el]) + elStr = str(el.normalize() if isinstance(el, Decimal) else el) + if count > 1: + elStr += "×" + str(count) + els.append(elStr) + + return "{" + ", ".join(els) + "}" + + +class ObjectsPartsStorage(PartsStorage): + """Stores parts objects explicitly""" + + __slots__ = () + + def add(self, part: Part) -> None: + s = self.sizes.get(part.len, None) + if s is None: + self.sizes[part.len] = s = set() + + s.add(part) + + def remove(self, part: Part) -> None: + self.sizes[part.len] -= {part} + if not self.sizes[part.len]: + del self.sizes[part.len] + + @classmethod + def _getFirstUncut(cls, partsOfTheLen: typing.Set[DirVectorPart]) -> typing.Optional[DirVectorPart]: + for part in partsOfTheLen: + if not part.cutted: + return part + + def get(self, size: NumT, uncut: bool = False) -> DirVectorPart: + partsOfTheLen = self.sizes[size] + if uncut: + return self.__class__._getFirstUncut(partsOfTheLen) + else: + return next(iter(partsOfTheLen)) + + def resize(self, part, newSize): + raise NotImplementedError + + def parts(self): + for size, sizeGroup in self.sizes.items(): + for part in sizeGroup: + yield part + + +USE_VECTORIZED = True + + +class SolutionState(Metered, TupleHashable): + __slots__ = ("piecesInitial", "piecesTarget") + DONT_REPR = frozenset(("initialPieces", "piecesTarget")) + + def __init__(self, task: typing.Dict[float, int], parts: typing.Dict[PartId, NumT]) -> None: + super().__init__() + self.metrics = MetricsSet() + self.metrics.quality.cutsCount = 0 + self.metrics.quality.joinsCount = 0 + self.metrics.path.lengthToAllocate = 0 + self.metrics.path.piecesInInitial = 0 + self.metrics.path.piecesInTarget = 0 + + if USE_VECTORIZED: + self.piecesInitial = ObjectsPartsStorage() + self.piecesTarget = ObjectsPartsStorage() + else: + self.piecesInitial = PartsStorage() + self.piecesTarget = PartsStorage() + + + self.metrics.path.lengthToAllocate = FloatType(0) + + uniqueParts = defaultdict(list) + for partId, partLen in parts.items(): + uniqueParts[partLen].append(partId) + self.metrics.path.lengthToAllocate += partLen + self.metrics.path.piecesInInitial += 1 + + targetTotal = FloatType(0) + uniqueTargets = defaultdict(list) + for i, (partLen, count) in enumerate(task.items()): + for j in range(count): + partId = PartId((partLen, j)) + uniqueTargets[partLen].append(partId) + self.metrics.path.piecesInTarget += count + targetTotal += partLen + + unusedLen = self.metrics.path.lengthToAllocate - targetTotal + + if unusedLen: + if unusedLen > 0: + uniqueTargets[unusedLen].append(PartId((inf, 0))) + self.metrics.path.piecesInTarget += 1 + else: + raise ValueError("The problem is unsolvable: length in the left is greater than the length in the right", unusedLen, self.metrics.lengthToAllocate, targetTotal, self.piecesInitial, self.piecesTarget) + + if USE_VECTORIZED: + b, U = vectorFactorization(np.hstack([list(uniqueParts), list(uniqueTargets)])) + + UParts = U[:, : len(uniqueParts)] + UTargets = U[:, -len(uniqueTargets) :] + + for i, (partLen, partIds) in enumerate(uniqueParts.items()): + partVec = UParts[:, i] + for partId in partIds: + self.piecesInitial.add(DirVectorPart(partId, partVec, b, 0)) + + for i, (partLen, partsOfThisLenIds) in enumerate(uniqueTargets.items()): + partVec = UTargets[:, i] + for partId in partsOfThisLenIds: + self.piecesTarget.add(DirVectorPart(partId, partVec, b, 1)) + + self.metrics.path.chunksToAnnihilate = np.sum(UParts) + else: + for i, (partLen, partIds) in enumerate(uniqueParts.items()): + for partId in partIds: + self.piecesInitial.add(DirPart(partId, partLen, 0)) + + for i, (partLen, partsOfThisLenIds) in enumerate(uniqueTargets.items()): + for partId in partsOfThisLenIds: + self.piecesTarget.add(DirPart(partId, partLen, 1)) + + self.metrics.path.chunksToAnnihilate = 0 + + def toTuple(self) -> typing.Tuple: + return (self.piecesInitial, self.piecesTarget) + + def __repr__(self) -> str: + return repr(self.piecesInitial) + ", " + repr(self.piecesTarget) + ";" + self.metrics.__rich_str__().plain() + + +class State(RichReprable, Metered): + __slots__ = ("children", "sizeToConsume", "parent", "metrics") + nodeAbbrev = None + + def __init__(self, sizeToConsume: typing.Any, parent: typing.Optional["State"]) -> None: + super().__init__() + self.metrics = MetricsSet() + + self.sizeToConsume = sizeToConsume + self.children = None + self.parent = parent + + @abstractmethod + def apply(self, solutionState: SolutionState) -> None: + raise NotImplemented + + def toTuple(self) -> typing.Tuple[typing.Type["State"], typing.Any, NumT]: + return (self.__class__, self.sizeToConsume,) + tuple(getattr(self, k) for k in self.__class__.__slots__) + + def __eq__(self, other: "State") -> bool: + return self.toTuple() == other.toTuple() + + def __hash__(self) -> int: + return hash(self.toTuple()) + + def generateUnfilteredChildren(self, solutionState: SolutionState, statesObserved: typing.Set[SolutionState]) -> None: + self.children = [] + + annihilateNodesMade = 0 + + #for n in Annihilate.make(self, solutionState): + for n in AnnihilateAll.make(self, solutionState): + annihilateNodesMade += 1 + yield n + + if not annihilateNodesMade: + cutNodesMade = 0 + #for n in Cut.make(self, solutionState): + for n in CutAndAnnihilateAll.make(self, solutionState): + yield n + + def generateChildren(self, solutionState: SolutionState, statesObserved: typing.Set[SolutionState]) -> None: + filteredChildren = [] + for c in self.generateUnfilteredChildren(solutionState, statesObserved): + ss = deepcopy(solutionState) + c.apply(ss) + #if ss.metrics not in statesObserved: + if ss in statesObserved: + continue + filteredChildren.append(c) + + self.children = filteredChildren + + def __repr__(self) -> str: + return str(self.__rich_repr__()) + + +class Annihilate(State): + __slots__ = () + nodeAbbrev = "A" + + def __init__(self, sizeToConsume: typing.Any, parent: typing.Optional["State"]) -> None: + super().__init__(sizeToConsume, parent) + + @classmethod + def updateMetrics(cls, metrics, partLen, piecesConsumedFromEach): + metrics.quality.cutsCount += 0 + metrics.quality.joinsCount += 0 + metrics.path.piecesInInitial -= piecesConsumedFromEach + metrics.path.piecesInTarget -= piecesConsumedFromEach + metrics.path.lengthToAllocate -= partLen + metrics.intize() + + @classmethod + def make(cls, parent: State, solutionState: SolutionState) -> typing.Iterable[State]: + for partLen in set(solutionState.piecesInitial.sizes) & set(solutionState.piecesTarget.sizes): + newState = cls(partLen, parent) + + newState.metrics = deepcopy(solutionState.metrics) + cls.updateMetrics(newState.metrics, partLen, 1) + + #solutionStateNew = deepcopy(solutionState) + #newState.apply(solutionStateNew) + yield newState + + def apply(self, solutionState: SolutionState) -> None: + solutionState.piecesTarget.pop(self.sizeToConsume) + solutionState.piecesInitial.pop(self.sizeToConsume) + + self.updateMetrics(solutionState.metrics, self.sizeToConsume, 1) + + def __rich_repr__(self) -> RichConsole.RichStr: + return RichConsole.groups.Fore.lightGreen(self.__class__.nodeAbbrev) + RichConsole.groups.Fore.lightGreen("<") + RichConsole.groups.Fore.cyan(repr(self.sizeToConsume)) + ";" + self.metrics.__rich_str__() + RichConsole.groups.Fore.lightGreen(">") + + +class AnnihilateAll(Annihilate): + __slots__ = () + nodeAbbrev = "Aa" + + def __init__(self, sizeToConsume: typing.Any, parent: typing.Optional["State"]) -> None: + super().__init__(sizeToConsume, parent) + + @classmethod + def make(cls, parent: State, solutionState: SolutionState) -> typing.Iterable[State]: + duplicated = set(solutionState.piecesInitial.sizes) & set(solutionState.piecesTarget.sizes) + if duplicated: + sumLen = sum(partLen for partLen in duplicated) + + newState = cls(None, parent) + newState.metrics = deepcopy(solutionState.metrics) + try: + cls.updateMetrics(newState.metrics, sumLen, len(duplicated)) + except BaseException: + print(solutionState, solutionState.piecesInitial.sizes, newState.metrics) + raise + + yield newState + + def getNodesToAnnihilate(self, solutionState): + print(solutionState.piecesInitial.sizes, solutionState.piecesTarget.sizes) + return set(solutionState.piecesInitial.sizes) & set(solutionState.piecesTarget.sizes) + + def apply(self, solutionState: SolutionState): + duplicated = self.getNodesToAnnihilate(solutionState) + sumLen = sum(partLen for partLen in duplicated) + + res = [] + for l in duplicated: + a = solutionState.piecesInitial.pop(l) + b = solutionState.piecesTarget.pop(l) + res.append((a, b)) + + self.updateMetrics(solutionState.metrics, sumLen, len(duplicated)) + + solutionState.metrics.intize() + return res + + def __rich_repr__(self) -> RichConsole.RichStr: + return RichConsole.groups.Fore.lightGreen(self.__class__.nodeAbbrev) + RichConsole.groups.Fore.lightGreen("<") + self.metrics.__rich_str__() + RichConsole.groups.Fore.lightGreen(">") + + +class CutAndAnnihilateAll(State): + """For visualization purposes only for now""" + + __slots__ = ("chunks", "direction") + nodeAbbrev = "CA" + + @classmethod + def updateMetricsCut(cls, metrics: MetricsSet, partLen: FloatType, direction: int) -> None: + if direction: + metrics.quality.joinsCount += 1 + metrics.path.piecesInTarget += 1 + else: + metrics.quality.cutsCount += 1 + metrics.path.piecesInInitial += 1 + metrics.intize() + + @classmethod + def updateMetricsAnnihilate(cls, metrics: MetricsSet, partLen: int, piecesConsumedFromEach: int, annihilatedChunks: int) -> None: + metrics.quality.cutsCount += 0 + metrics.quality.joinsCount += 0 + metrics.path.piecesInInitial -= piecesConsumedFromEach + metrics.path.piecesInTarget -= piecesConsumedFromEach + metrics.path.lengthToAllocate -= partLen + metrics.path.chunksToAnnihilate -= annihilatedChunks + metrics.intize() + + @classmethod + def enumerateVectorizedCuts(cls, parent: State, solutionState: SolutionState, storage: ObjectsPartsStorage, direction: int) -> typing.Iterator[State]: + for partLenT, partsOfTheLen in storage.sizes.items(): + part = ObjectsPartsStorage._getFirstUncut(partsOfTheLen) + if part is None: + continue + + for chunks in part.getCuttings(): + #print("firstPieceVec", firstPieceVec, "second", second) + newState = cls(part.len, parent, direction, chunks) + ss1 = deepcopy(solutionState) + newState.apply(ss1) + newState.metrics = ss1.metrics + yield newState + + @classmethod + def makeVectorized(cls, parent: State, solutionState: SolutionState) -> typing.Iterable[State]: + for partLenI, iPartsOfSize in solutionState.piecesInitial.sizes.items(): + iPartFirst = next(iter(iPartsOfSize)) + iPartUncut = ObjectsPartsStorage._getFirstUncut(iPartsOfSize) + + for partLenT, tPartsOfSize in solutionState.piecesTarget.sizes.items(): + if partLenI < partLenT: + direction = 1 + partLen = partLenT + tPart = ObjectsPartsStorage._getFirstUncut(tPartsOfSize) + if tPart is None: + continue + iPart = iPartFirst + cutPart = tPart + elif partLenI > partLenT: + direction = 0 + partLen = partLenI + tPart = next(iter(tPartsOfSize)) + iPart = iPartUncut + cutPart = iPart + else: + continue + + if iPart is None: + continue + + commonPieces = iPart.pieces & tPart.pieces + if commonPieces.any(): + #for firstPieceVec in DirVectorPart._getCuttings(commonPieces): + # #print(firstPieceL, partLen, secondPieceL) + + # #metric = addedWaste + + # print("firstPieceVec", firstPieceVec) + # newState = cls(partLen, parent, direction, firstPieceVec) + # ss1 = deepcopy(solutionState) + # newState.apply(ss1) + # newState.metrics = ss1.metrics + # yield newState + + + #print("firstPieceVec", firstPieceVec) + newState = cls(partLen, parent, direction, (commonPieces, cutPart.pieces & ~commonPieces)) + ss1 = deepcopy(solutionState) + newState.apply(ss1) + newState.metrics = ss1.metrics + yield newState + + @classmethod + def makeVectorized(cls, parent: State, solutionState: SolutionState) -> typing.Iterable[State]: + yield from cls.enumerateVectorizedCuts(parent, solutionState, solutionState.piecesInitial, 0) + yield from cls.enumerateVectorizedCuts(parent, solutionState, solutionState.piecesTarget, 1) + + @classmethod + def makeNonVectorized(cls, parent: State, solutionState: SolutionState) -> typing.Iterable[State]: + for partLenI in solutionState.piecesInitial.sizes: + for partLenT in solutionState.piecesTarget.sizes: + if partLenI < partLenT: + direction = 1 + elif partLenI > partLenT: + direction = 0 + else: + continue + + firstPieceL = min(partLenI, partLenT) + partLen = max(partLenI, partLenT) + secondPieceL = partLen - firstPieceL + + if firstPieceL < secondPieceL: + # continue + firstPieceL, secondPieceL = secondPieceL, firstPieceL + + #print(firstPieceL, partLen, secondPieceL) + + #metric = addedWaste + + newState = cls(partLen, parent, direction, (firstPieceL,)) + ss1 = deepcopy(solutionState) + newState.apply(ss1) + newState.metrics = ss1.metrics + + #solutionStateNew = deepcopy(solutionState) + #newState.apply(solutionStateNew) + yield newState + + @classmethod + def make(cls, parent: State, solutionState: SolutionState) -> typing.Iterable[State]: + if USE_VECTORIZED: + return cls.makeVectorized(parent, solutionState) + else: + return cls.makeNonVectorized(parent, solutionState) + + def __init__(self, sizeToConsume: typing.Any, parent: State, direction: bool, chunks: NumT) -> None: + super().__init__(sizeToConsume, parent) + self.chunks = chunks + self.direction = direction + + def apply(self, solutionState: SolutionState) -> None: + piecesDb = solutionState.piecesTarget if self.direction else solutionState.piecesInitial + part = piecesDb.pop(self.sizeToConsume) + self.updateMetricsCut(solutionState.metrics, part.len, self.direction) + + chunks = part.split(*self.chunks) + + for c in chunks: + piecesDb.add(c) + + # annihilate step + + sumLen = 0 + res = [] + annihilatedChunks = 0 + + duplicated0 = set(solutionState.piecesInitial.sizes) & set(solutionState.piecesTarget.sizes) + + duplicated = duplicated0 + + while duplicated: + sumLen += sum(partLen for partLen in duplicated) + + for l in duplicated: + a = solutionState.piecesInitial.pop(l) + b = solutionState.piecesTarget.pop(l) + if USE_VECTORIZED: + annihilatedChunks += np.sum(a.pieces) + res.append((a, b)) + + duplicated = set(solutionState.piecesInitial.sizes) & set(solutionState.piecesTarget.sizes) + + self.updateMetricsAnnihilate(solutionState.metrics, sumLen, len(duplicated0), annihilatedChunks) + + solutionState.metrics.intize() + return res + + def __rich_repr__(self) -> RichConsole.RichStr: + return RichConsole.groups.Fore.lightGreen(self.__class__.nodeAbbrev) + RichConsole.groups.Fore.lightGreen("<") + RichConsole.groups.Fore.cyan(repr(self.sizeToConsume)) + RichConsole.groups.Fore.lightGreen(" <- " if self.direction else " -> ") + RichConsole.groups.Fore.cyan(repr(self.chunks)) + ";" + self.metrics.__rich_str__() + RichConsole.groups.Fore.lightGreen(">") + + +class Root(State): + __slots__ = () + + def __init__(self) -> None: + super().__init__(None, None) + + def __rich_repr__(self) -> RichConsole.RichStr: + return RichConsole.groups.Fore.lightGreen(self.__class__.__name__) + + def apply(self, solutionState: SolutionState) -> None: + pass + + +# from matplotlib import pyplot as plt +_debugNx = nx.DiGraph() + + +class SearchStrategy: + __slots__ = ("frontier",) + + def __init__(self, initialState): + raise NotImplementedError + + def push(self, v): + raise NotImplementedError + + def pop(self): + raise NotImplementedError + + def prune(self, bestSolutionMetrics): + raise NotImplementedError + + def __bool__(self) -> bool: + return bool(self.frontier) + + def __repr__(self): + return repr(self.frontier) + + +class BeFS(SearchStrategy): + __slots__ = () + + def __init__(self, initialState: Root) -> None: + self.frontier = [initialState] + + def push(self, v: State) -> None: + heappush(self.frontier, v) + + def pop(self) -> State: + return heappop(self.frontier) + + @classmethod + def findMiddleElement(cls, pivotA, pivotB, predicate): + while pivotA + 1 < pivotB: + pivotMiddle = (pivotA + pivotB) // 2 + if predicate(pivotMiddle): + pivotB = pivotA + else: + pivotA = pivotMiddle + return pivotMiddle + + def prune(self, bestSolutionMetrics): + self.frontier = sorted(self.frontier, key=lambda el: el.metrics.quality) + pivotMiddle = self.findMiddleElement(0, len(self.frontier) - 1, (lambda pivotMiddle: self.frontier[pivotMiddle].metrics.quality > bestSolutionMetrics.quality)) + self.frontier = sorted(self.frontier[:pivotMiddle]) + + +class GBeFS(BeFS): + __slots__ = () + + def pop(self) -> State: + res = self.frontier[0] + self.frontier = [] + return res + + def prune(self, bestSolutionMetrics): + pass + + +class BrFS(SearchStrategy): + __slots__ = () + + def __init__(self, initialState: Root) -> None: + self.frontier = deque((initialState,)) + + def push(self, v: State) -> None: + self.frontier.append(v) + + def pop(self) -> State: + return self.frontier.popleft() + + def prune(self, bestSolutionMetrics): + pass + + +class DFS(SearchStrategy): + __slots__ = () + + def __init__(self, initialState: Root) -> None: + self.frontier = deque((initialState,)) + + def push(self, v: State) -> None: + self.frontier.append(v) + + def pop(self) -> State: + return self.frontier.pop() + + def prune(self, bestSolutionMetrics): + pass + + +class Algo: + __slots__ = ("initialState", "initialSolutionState", "statesObserved", "solutions", "bestSolutionMetrics") + + def __init__(self, task: typing.Dict[int, int], parts: typing.Dict[typing.Any, NumT]) -> None: + self.initialState = Root() + self.initialSolutionState = SolutionState(task, parts) + self.statesObserved = set() + self.solutions = [] + self.bestSolutionMetrics = MetricsSet() + for smName in self.bestSolutionMetrics.__class__.__slots__: + subMetrics = getattr(self.bestSolutionMetrics, smName) + for k in subMetrics.__class__.__slots__: + setattr(subMetrics, k, inf) + + # from analysis + self.bestSolutionMetrics.quality.cutsCount = self.initialSolutionState.metrics.path.piecesInInitial - 1 + self.bestSolutionMetrics.quality.joinsCount = self.initialSolutionState.metrics.path.piecesInTarget - 1 + + def getPath(self, node: State) -> typing.Iterable[State]: + path = [] + while node: + path.append(node) + node = node.parent + + path = list(reversed(path[:-1])) # without root node + return path + + def getSolutionStateFromPath(self, path: typing.Iterable[State]) -> SolutionState: + fs = deepcopy(self.initialSolutionState) + + for node in path: + node.apply(fs) + return fs + + def getSolutionFromPath(self, path: typing.Iterable[State]): + solution = initSolution() + fs = deepcopy(self.initialSolutionState) + + for node in path: + if isinstance(node, (CutAndAnnihilateAll, Annihilate)): + for iP, tP in node.apply(fs): + print(tP.id, tP.id.origId, tP.id._origId) + origTargetSize = tP.id._origId[0] # we encode sizes of targets in their ids as the first component + if origTargetSize == inf: + continue + origTargetInSlotId = tP.id._origId[1] + solution[origTargetSize][origTargetInSlotId].append((iP.id.origId, iP.len)) + else: + node.apply(fs) + return solution + + def getSolutionState(self, node: State) -> SolutionState: + return self.getSolutionStateFromPath(self.getPath(node)) + + def goalTest(self, solutionState: SolutionState) -> bool: + #return not bool(set(solutionState.piecesInitial.sizes) - set(solutionState.piecesTarget.sizes)) + return not bool(solutionState.piecesInitial.sizes) + + def appendSolution(self, mostPromisingNode: State, fs: SolutionState, frontier: list) -> None: + heappush(self.solutions, mostPromisingNode) + needPruneFrontier = False + for smName in self.bestSolutionMetrics.__class__.__slots__: + subMetrics = getattr(fs.metrics, smName) + subMetricsBest = getattr(self.bestSolutionMetrics, smName) + for k in subMetrics.__class__.__slots__: + bestKnown = getattr(subMetricsBest, k) + newOne = getattr(subMetrics, k) + if bestKnown > newOne: + bestKnown = newOne + setattr(subMetricsBest, k, bestKnown) + needPruneFrontier = True + + if needPruneFrontier: + frontier.prune(self.bestSolutionMetrics) + + def search(self, strategy: SearchStrategy, limit: int, maxSolutions: int, debug: bool = False) -> None: + frontier = strategy(self.initialState) + + self.statesObserved = set() + + fs = deepcopy(self.initialSolutionState) + i = 0 + while True: + if not frontier: + break + mostPromisingNode = frontier.pop() + print(mostPromisingNode) + fs = self.getSolutionState(mostPromisingNode) + + #self.statesObserved |= {fs.metrics} + self.statesObserved |= {fs} + #_debugNx.add_node(mostPromisingNode.__rich_repr__().plain(), i=i) + + parentState = self.getSolutionState(mostPromisingNode.parent) + _debugNx.add_edge(repr(parentState), repr(fs), i=i) + #_debugNx.add_edge(parentState.metrics.__rich_str__().plain(), fs.metrics.__rich_str__().plain(), i=i) + + if self.goalTest(fs): + _debugNx.add_edge(repr(fs), "solved") + #_debugNx.add_edge(repr(fs.metrics.__rich_str__().plain()), "solved") + self.appendSolution(mostPromisingNode, fs, frontier) + + if len(self.solutions) >= maxSolutions: + break + + networkx.drawing.nx_agraph.write_dot(_debugNx, "./graph.dot") + + i += 1 + + mostPromisingNode.generateChildren(fs, self.statesObserved) + + for child in mostPromisingNode.children: + #print("path", child, child.metrics.path, ">", self.bestSolutionMetrics.path, child.metrics.path > self.bestSolutionMetrics.path) + #print("quality", child, child.metrics.quality, ">", self.bestSolutionMetrics.quality, child.metrics.quality > self.bestSolutionMetrics.quality) + if child.metrics.quality > self.bestSolutionMetrics.quality: + continue + frontier.push(child) + + if i > limit: + raise Exception("Fuck") + + print("Nodes expanded", i) + + networkx.drawing.nx_agraph.write_dot(_debugNx, "./graph.dot") + + def __call__(self) -> typing.List[State]: + self.search(DFS, 100, 1) + #print("DFS", self.solutions) + self.search(BeFS, 1000, 100, True) + print("BeFS", self.solutions) + return self.solutions + + def __repr__(self): + return self.__class__.__name__ + "(" + ", ".join((k + "=" + repr(getattr(self, k))) for k in self.__class__.__slots__) + ")" diff --git a/timberAllocation/utils/__init__.py b/timberAllocation/utils/__init__.py new file mode 100644 index 0000000..a5276b0 --- /dev/null +++ b/timberAllocation/utils/__init__.py @@ -0,0 +1,77 @@ +import typing +from abc import ABC, abstractmethod + +import RichConsole + + +class RichReprable(ABC): + __slots__ = () + + @abstractmethod + def __rich_repr__(self) -> RichConsole.RichStr: + raise NotImplementedError + + def __repr__(self) -> str: + return str(self.__rich_repr__()) + + +class RichStringable(ABC): + __slots__ = () + + @abstractmethod + def __rich_str__(self) -> RichConsole.RichStr: + raise NotImplementedError + + def __str__(self) -> str: + return str(self.__rich_str__()) + + +class RichStringableAndReprable(RichReprable, RichStringable): + __slots__ = () + + def __rich_repr__(self) -> RichConsole.RichStr: + return RichConsole.groups.Fore.lightGreen(self.__class__.__name__) + RichConsole.groups.Fore.lightGreen("<") + self.__rich_str__() + RichConsole.groups.Fore.lightGreen(">") + + +class TupleHashable(ABC): + __slots__ = () + + @abstractmethod + def toTuple(self) -> typing.Tuple: + raise NotImplementedError() + + def __eq__(self, other: "TupleHashable") -> bool: + #if not isinstance(other, tuple): + ot = other.toTuple() + #else: + # ot = other + return self.toTuple() == ot + + def __hash__(self) -> int: + return hash(self.toTuple()) + + +class StructuredComparable(ABC): + __slots__ = () + + @abstractmethod + def _coIter(self, other) -> typing.Iterable[typing.Tuple[str, typing.Any, typing.Any]]: + raise NotImplementedError + + def _doComparison(self, other: "StructuredComparable", predicate: typing.Callable) -> bool: + res = True + for k, sv, ov in self._coIter(other): + res = res and predicate(sv, ov) + return res + + def __lt__(self, other: "StructuredComparable") -> bool: + return self._doComparison(other, lambda sv, ov: sv < ov) + + def __gt__(self, other: "StructuredComparable") -> bool: + return self._doComparison(other, lambda sv, ov: sv > ov) + + def __le__(self, other: "StructuredComparable") -> bool: + return self._doComparison(other, lambda sv, ov: sv <= ov) + + def __ge__(self, other: "StructuredComparable") -> bool: + return self._doComparison(other, lambda sv, ov: sv >= ov) diff --git a/timberAllocation/utils/z3Numpy.py b/timberAllocation/utils/z3Numpy.py new file mode 100644 index 0000000..61d9147 --- /dev/null +++ b/timberAllocation/utils/z3Numpy.py @@ -0,0 +1,363 @@ +__all__ = ("ndarray", "array") + +import sys + +import numpy as np +import z3 + +z3Internal = sys.modules["z3.z3"] + +_get_argsBackup = z3Internal._get_args + + +def _get_argsMonkeyPatched(args): + try: + argsNew = [] + for arg in args: + if isinstance(arg, np.ndarray): + argsNew.extend(arg.flatten()) + else: + argsNew.append(arg) + return _get_argsBackup(argsNew) + except BaseException: + return _get_argsBackup(args) + + +z3Internal._get_args = _get_argsMonkeyPatched + + +def z3OptimizeDecorate(f): + def optimizeF(self, *args, **kwargs): + if isinstance(args[0], np.ndarray): + for el in args[0].flat: + f(self, el, *args[1:], **kwargs) + + return optimizeF + + +z3.Optimize.minimize = z3OptimizeDecorate(z3.Optimize.minimize) +z3.Optimize.maximize = z3OptimizeDecorate(z3.Optimize.maximize) + + +def generateFiniteLenInts(name, count, tp, ctx): + tp = np.dtype(tp) + bitsCount = tp.itemsize * 8 + return [z3.BitVec(name + "__" + str(i), bitsCount) for i in range(count)] + + +def generateFiniteLenFloats(name, count, tp, ctx): + tp = np.dtype(tp) + fpSort = floatTypes[tp.itemsize] + if isinstance(fpSort, tuple): + fpSort = z3.FPSort(*fpSort) + return [z3.FP(name + "__" + str(i), fpSort) for i in range(count)] + + +typesRemapping = { + np.bool_: lambda name, count, tp, ctx: z3.BoolVector(name, count, ctx), + bool: lambda name, count, tp, ctx: z3.BoolVector(name, count, ctx), + int: lambda name, count, tp, ctx: z3.IntVector(name, count, ctx), + float: lambda name, count, tp, ctx: z3.RealVector(name, count, ctx), +} + +floatTypes = { + 1: (4, 4), + #2: (5, 11), + 2: z3.FloatHalf(), + #4: (8, 24), + 4: z3.FloatSingle(), + #8: (11, 53), + 8: z3.FloatDouble(), + 10: (15, 63), + #16: (15, 111), + 16: z3.FloatQuadruple(), + 32: (19, 237), +} + +supertypesRemapping = { + np.integer: generateFiniteLenInts, + np.floating: generateFiniteLenFloats, +} + +backRemapSortsToTypes = { + z3.Z3_BOOL_SORT: bool, + z3.Z3_INT_SORT: int, + z3.Z3_REAL_SORT: float, +} + + +def _getBackRemapSortsToDTypes(x): + sort = x.sort() + if isinstance(x, z3.FPRef): + size = sort.sbits() + sort.ebits() + return np.dtype("float" + str(size)) + elif isinstance(x, z3.BitVecRef): + return np.dtype("int" + str(sort.size())) + else: + return np.dtype(backRemapSortsToTypes[sort.kind()]) + + +def _makeOurNdArray(underlyingArray, shape, order=None): + res = ourNdArray(shape=shape, order=order, dtype=object) + res.flat = np.array(underlyingArray, order=order, subok=True) + return res + + +class ndarray(np.ndarray): + __slots__ = () + + def _applyElementwise(self, another, pred): + return _makeOurNdArray([pred(s, a) for (s, a) in np.broadcast(self, another)], self.shape) + + def __eq__(self, another): + return self._applyElementwise(another, lambda s, a: s == a) + + def __ne__(self, another): + return self._applyElementwise(another, lambda s, a: s != a) + + def __ge__(self, another): + return self._applyElementwise(another, lambda s, a: s >= a) + + def __le__(self, another): + return self._applyElementwise(another, lambda s, a: s <= a) + + def __gt__(self, another): + return self._applyElementwise(another, lambda s, a: s > a) + + def __lt__(self, another): + return self._applyElementwise(another, lambda s, a: s < a) + + def __and__(self, another): + if isinstance(self.flat[0], z3.BoolRef): + return self._applyElementwise(another, z3.And) + else: + return self._applyElementwise(another, lambda s, a: s & a) + + def __or__(self, another): + if isinstance(self.flat[0], z3.BoolRef): + return self._applyElementwise(another, z3.Or) + else: + return self._applyElementwise(another, lambda s, a: s | a) + + def __xor__(self, another): + if isinstance(self.flat[0], z3.BoolRef): + return self._applyElementwise(another, z3.Xor) + else: + return self._applyElementwise(another, lambda s, a: s ^ a) + + def __invert__(self, another): + if isinstance(self.flat[0], z3.BoolRef): + return self._applyElementwise(None, lambda s, a: z3.Not(s)) + else: + return self._applyElementwise(another, lambda s, a: ~s) + + def toN(self): + dt = _getBackRemapSortsToDTypes(self.flat[0]) + res = np.ndarray(shape=self.shape, dtype=dt) + uintView = res.view(dtype=np.dtype("u" + str(dt.itemsize))) + for i, elRes in enumerate(self.flat): + if isinstance(elRes, z3.RatNumRef): + res.flat[i] = elRes.as_fraction() + elif isinstance(elRes, z3.z3.FPNumRef): + uintView.flat[i] = fpNumRefToFloatInt(elRes) + else: + res.flat[i] = elRes + return res + + +ourNdArray = ndarray + + +def arrayCreateFromNameAndShape(name, shape, dtype=bool, order=None, ctx=None): + ctor = typesRemapping.get(dtype, None) + if ctor is None: + for superType, ctorCandidate in supertypesRemapping.items(): + if issubclass(dtype, superType): + ctor = ctorCandidate + break + + flat = ctor(name, np.product(shape), dtype, ctx=ctx) + return _makeOurNdArray(flat, shape, order) + + +def arrayCreateFromExistingArray(*args, **kwargs): + res = np.array(*args, **kwargs) + res1 = ourNdArray(shape=res.shape, order=None, dtype=res.dtype) + res1.flat = res.flat + del res + return res1 + + +def array(*args, **kwargs): + if isinstance(args[0], str) or isinstance(kwargs.get("name"), str): + return arrayCreateFromNameAndShape(*args, **kwargs) + else: + return arrayCreateFromExistingArray(*args, **kwargs) + + +array.__wraps__ = arrayCreateFromNameAndShape + +simplifyBackup = z3.simplify + + +def newSimplify(a, *arguments, **keywords): + if isinstance(a, np.ndarray): + return _makeOurNdArray([simplifyBackup(e, *arguments, **keywords) for e in a.flat], a.shape) + else: + simplifyBackup(a, *arguments, **keywords) + + +newSimplify.__wraps__ = simplifyBackup +z3.simplify = newSimplify + + +def fpNumRefToFloatInt(a): + significandBitsWithoutSign = a.sbits() - 1 + signRemoveMask = (1 << significandBitsWithoutSign) - 1 + # sign exponent significand + return ( + ( + (int(a.isNegative()) << a.ebits()) + | + a.exponent_as_long(biased=True) + ) << significandBitsWithoutSign + ) | (a.significand_as_long() & signRemoveMask) + +if not hasattr(z3.IntNumRef, "__int__"): + z3.IntNumRef.__int__ = z3.IntNumRef.as_long + +if not hasattr(z3.FPNumRef, "__float__"): + + def FPNumRefToNumpyFloat(self): + dtF = np.dtype("f" + str((self.ebits() + self.sbits()) // 8)) + dtI = np.dtype("u" + str(dtF.itemsize)) + f = np.ndarray((1,), dtype=dtF) + i = f.view(dtype=dtI) + i[0] = fpNumRefToFloatInt(self) + return f[0] + + def numpyFPNumRefToFloat(self): + return float(self.toNumpyFloat()) + + z3.FPNumRef.toNumpyFloat = FPNumRefToNumpyFloat + z3.FPNumRef.__float__ = numpyFPNumRefToFloat + + +z3IntValBackup = z3.IntVal + + +def z3IntValPatched(v, *args, **kwargs): + if isinstance(v, np.integer): + v = int(v) + return z3IntValBackup(v, *args, **kwargs) + + +z3IntValPatched.__wraps__ = z3IntValBackup +z3Internal.IntVal = z3.IntVal = z3IntValPatched + +z3FPValBackup = z3.FPVal + + +def z3FPValPatched(v, *args, **kwargs): + if isinstance(v, np.integer): + v = int(v) + elif isinstance(v, np.floating): + v = float(v) + return z3FPValBackup(v, *args, **kwargs) + + +z3FPValPatched.__wraps__ = z3FPValBackup +z3Internal.FPVal = z3.FPVal = z3FPValPatched + +z3RealValBackup = z3.RealVal + + +def z3RealValPatched(v, *args, **kwargs): + if isinstance(v, np.integer): + v = int(v) + elif isinstance(v, np.floating): + v = float(v) + return z3RealValBackup(v, *args, **kwargs) + + +z3RealValPatched.__wraps__ = z3RealValBackup +z3Internal.RealVal = z3.RealVal = z3RealValPatched + + +def Abs(self): + return z3.If(self >= 0, self, -self) + + +if not hasattr(z3.ArithRef, "__abs__"): + z3.ArithRef.__abs__ = Abs + +if not hasattr(z3.FPRef, "__abs__"): + z3.FPRef.__abs__ = Abs + +if not hasattr(z3.BoolRef, "__int__"): + + def z3BoolRefToInt(self): + return z3.If(self, 1, 0) + + z3.BoolRef.__int__ = z3BoolRefToInt + +if not hasattr(z3.BoolRef, "__add__"): + + def z3BoolRefAdd(self, another): + return self.__int__() + another.__int__() + + z3.BoolRef.__add__ = z3BoolRefAdd + +ModelRefGetItemBackup = z3.ModelRef.__getitem__ + + +def __getitem__Patched(self, idx): + if isinstance(idx, np.ndarray): + dt = _getBackRemapSortsToDTypes(idx.flat[0]) + res = ourNdArray(shape=idx.shape, dtype=object, subok=True) # our ndarray + for i, el in enumerate(idx.flat): + res.flat[i] = ModelRefGetItemBackup(self, el) + return res + else: + return ModelRefGetItemBackup(self, idx) + + +z3.ModelRef.__getitem__ = __getitem__Patched + + +_to_float_str_backup = z3Internal._to_float_str + + +def _to_float_str(*args, **kwargs): + if isinstance(args[0], np.number): + args = list(args) + args[0] = float(args[0]) + + return _to_float_str_backup(*args, **kwargs) + + +z3Internal._to_float_str = _to_float_str + + +""" +import numpy as np +import z3 +from timberAllocation.utils import z3Numpy + +x = z3Numpy.array("x", (2, 2), dtype=np.float16) +y = z3Numpy.array("y", (2, 2), dtype=np.float16) + +s = z3.Solver() +s.add(x @ y == np.identity(x.shape[0])) +s.add(abs(x) < 10) +s.add(abs(y) < 10) +s.add(abs(x) > 0.1) +s.add(abs(y) > 0.1) + +if s.check() == z3.sat: + mod = s.model() + xx = mod[x] + yy = mod[y] + print(xx, yy, z3.simplify(xx @ yy)) + print(xx.toN(), yy.toN(), xx.toN() @ yy.toN()) +""" diff --git a/timberAllocation/vis.py b/timberAllocation/vis.py new file mode 100644 index 0000000..830f192 --- /dev/null +++ b/timberAllocation/vis.py @@ -0,0 +1,84 @@ +import colorsys +import typing +from collections import defaultdict +from decimal import Decimal + +from defaultlist import defaultlist +from numpy import float32, float64 + +import RichConsole +from RichConsole import RichStr + +from . import FloatType, NumT, PartId + + +def makeLogPart(partsColors, zoom: float, chunkL: NumT, label: str, partId: PartId, char: str = "-") -> typing.Union[RichStr, str]: + size = int(round(chunkL * zoom) - len(label)) + halfSize = int(round(size / 2)) + text = char * halfSize + label + char * (size - halfSize) + if partId is not None: + return RichConsole.RGBColor(None, *partsColors[partId.origId])(text) + else: + return text + + +def makeLog(parts: typing.Iterable[RichStr]) -> RichStr: + return "|" + RichConsole.rsjoin("|", parts) + "|" + + +def visualizeSolution(solution, partsColors, zoom=FloatType(0.75)) -> None: + for slotL, combineds in sorted(solution.items(), key=lambda x: x[0], reverse=True): + totalRealL = FloatType(0) + for combined in combineds: + l = [] + for chunk in combined: + chunkId, chunkL = chunk + totalRealL += chunkL + label = repr(chunkId) + "," + str(chunkL) + l.append(makeLogPart(partsColors, zoom, chunkL, label, chunkId)) + if totalRealL < slotL: + lack = slotL - totalRealL + l.append(rgbBlack(makeLogPart(partsColors, zoom, lack, str(lack), None, char="."))) + print(slotL, makeLog(l)) + + +rgbBlack = RichConsole.RGBColor(0, 0, 0, bg=True) # simple bg colors don't work in Jupyter + + +def visualizeCutting(parts, solution, partsColors, zoom=FloatType(0.75)) -> None: + cutting = defaultdict(list) + uncut = set(parts) + for slotL, combineds in solution.items(): + for combined in combineds: + for chunk in combined: + chunkL = chunk[1] + cutting[chunk[0]].append(chunk[1]) + + for partId, chunkLs in cutting.items(): + l = [] + for chunkL in chunkLs: + label = str(chunkL) + l.append(makeLogPart(partsColors, zoom, chunkL, label, partId)) + waste = parts[partId.origId] - sum(chunkLs) + if waste: + l.append(rgbBlack(makeLogPart(partsColors, zoom, waste, str(waste), partId, char="&"))) + print(partId, parts[partId.origId], makeLog(l)) + uncut -= {partId} + + for partId, partL in sorted(((partId, parts[partId.origId]) for partId in uncut), key=lambda x: x[1], reverse=True): + print(partId, parts[partId.origId], RichConsole.groups.Underline.underline((makeLog((makeLogPart(partsColors, zoom, partL, str(partL), partId, char="="),))))) + + +def generateColors(parts: typing.Dict[PartId, FloatType]) -> typing.Dict[PartId, typing.Tuple[int, int, int]]: + step = 1.0 / len(parts) + partsColors = {} + for i, p in enumerate(parts): + partsColors[p.origId] = tuple((int(round(255 * c)) for c in colorsys.hsv_to_rgb(i / 2 * step, 1 - 0.3 * (i % 3 == 0), 0.5 + 0.5 * (i % 2)))) + return partsColors + + +def visualizeAll(solution, parts, partsColors, zoom=FloatType(0.75)): + print("\nSolution") + visualizeSolution(solution, partsColors, zoom) + print("\nCutting") + visualizeCutting(parts, solution, partsColors, zoom) diff --git a/timberAllocation/z3Solution.py b/timberAllocation/z3Solution.py new file mode 100644 index 0000000..eedab09 --- /dev/null +++ b/timberAllocation/z3Solution.py @@ -0,0 +1,60 @@ +import typing + +import z3 +from timberAllocation.utils import z3Numpy +import numpy as np + +from . import DirPart, DirVectorPart, FloatType, NumT, Part, PartId, Slot, getParts, getSlots, initSolution +from .middleRepr import vectorFactorization + +def z3Solution(task: typing.Dict[NumT, int], parts: typing.Dict[PartId, NumT]): + intial = [] + initialTotal = 0 + for partId, partLen in parts.items(): + intial.append(partLen) + initialTotal += partLen + + target = [] + + targetTotal = 0 + for i, (partLen, count) in enumerate(task.items()): + for j in range(count): + partId = PartId((partLen, j)) + target.append(partLen) + targetTotal += partLen + + unusedLen = initialTotal - targetTotal + + if unusedLen: + if unusedLen > 0: + target.append(unusedLen) + else: + raise ValueError("The problem is unsolvable: length in the left is greater than the length in the right", unusedLen) + + intial = np.array(intial) + target = np.array(target) + + it = np.hstack([intial, target]) + print(it, intial.shape[0], target.shape[0]) + + middle = z3Numpy.array("middle", (it.shape[0] - 1,), dtype=float) + alphaI = z3Numpy.array("ai", (middle.shape[0], intial.shape[0]), dtype=bool) + alphaT = z3Numpy.array("at", (middle.shape[0], target.shape[0]), dtype=bool) + alpha = z3Numpy.array(np.hstack([alphaI,alphaT])) # np.hstack always returns np.array + print(alpha) + + + s = z3.Optimize() + + s.add(middle @ alpha == it) + print((alpha.sum(axis=0)).__class__, (alpha.sum(axis=0) == 1).__class__) + #s.add(alpha.sum(axis=0)) + + joins = alphaT.sum() + cuts = alphaI.sum() + s.minimize(cuts) + s.minimize(joins) + + res = s.check() + print(res) +