From adfd34e861f8427a3544438b1bd885c1f20c1c76 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Tue, 21 Jun 2022 17:59:52 +0200 Subject: [PATCH 01/20] C api documentation --- c-documentation.md | 182 ++++++++++++++++++++++++++++++++++++++++ python-documentation.md | 36 ++++++++ 2 files changed, 218 insertions(+) create mode 100644 c-documentation.md create mode 100644 python-documentation.md diff --git a/c-documentation.md b/c-documentation.md new file mode 100644 index 0000000..8c6cf11 --- /dev/null +++ b/c-documentation.md @@ -0,0 +1,182 @@ +# C-API for Pickhardt Payments + +## Overview + +This document describes *Pickhardt payments* C library API. + +## Primitive data types + +Nodes ID are encoded as 33 bytes numbers. +``` +typedef char[33] nodeID_t; +``` + +Channels ID are encoded as 8 bytes numbers. +``` +typedef uint64_t channelID_t; +``` + +## Payment Session + +A `paymentSession` is an object that contains all the information known about the Lightning Network +from the initial knowledge of the public channels and the accumulated information +resulting from successful and failed attempts to forward payments while the session is in use. +``` +typedef struct paymentSession paymentSession; +``` + +To create and initialize a `paymentSession` the function `paymentSession_new` must be called. +This function returns a pointer to a newly allocated `paymentSession`. +This function returns `NULL` if the necessary memory couldn't be allocated. +``` +paymentSession* paymentSession_new(); +``` + +To cleanup and release the memory one needs to call `paymentSession_delete`. +This function does not fail. +``` +void paymentSession_delete(paymentSession* s); +``` + +Internally a Session represents the Lightning-Network graph of nodes and channels. +Initially a newly created Session has no knowled of any channel nor nodes, so one has to provide +this information by means of the function `paymentSession_addChannel`. +This function returns an error code that could have the values: +`PAYMENTSESSION_SUCCESS` if no errors occurred or +`PAYMENTSESSION_BADALLOC` if this function failed to allocate the memory for the new channel. +``` +int paymentSession_addChannel(paymentSession* s, + channelID_t short_channel_id, + nodeID_t Source, + nodeID_t Dest, + int64_t capacity, + int64_t fee_rate, + int64_t base_fee); +``` + +A channel can be removed by calling `paymentSession_rmChannel`. +This function returns an error code. +``` +int paymentSession_rmChannel(paymentSession* s, + channelID_t short_channel_id, + nodeID_t Source, + nodeID_t Dest); +``` + +A Session can be used to compute a multipart-payment with maximum probability of success based on +the knowledge of the Lightning-Network topology, channels capacity and previous knowledge of the +state of the liquidity in the channels. +By calling `paymentSession_optimizedPayment` one obtains a pointer to a `multiPartPayment` object +that contains a multipart-payment canditate based on the `amount` requested to be forwarded between +nodes `Source` and `Dest`. +This function returns `NULL` if the `paymentSession_optimizedPayment` fails. +The variable `status` contains an error code to log the cause of the failure. +``` +multiPartPayment* paymentSession_optimizedPayment(paymentSession* s, + nodeID_t Source, + nodeID_t Dest, + int64_t amount, + char* status); +``` + +Once an onion is sent to the Lightning-Network from the success or failure of this attempt we can +update our knowledge of the liquidity of the channels involved in that onion. That information can +be used by the Session to compute future payments attempts. +To update the knowledge of the Network one must call `paymentSession_updateKnowledge`. +This function returns an error code. +``` +int paymentSession_updateKnowledge(paymentSession* s, + int64_t short_channel_id, + nodeID_t Source, + nodeID_t Dest, + int64_t amount, + char status); +``` + +The accumulated knowledge in the Session can be erased by calling +`paymentSession_forgetInformation`. +This function returns an error code. +``` +int paymentSession_forgetInformation(paymentSession* s); +``` + +## Multipart Payments + +The function `paymentSession_optimizedPayment` produces a multipart payment (MPP) based on the +knowledge accumulated in the Session in use. The object that implements the concept of that MPP is +the `multiPartPayment`. We do not expose a constructor for this object in the public interface, +because it will always represent the result of an `paymentSession_optimizedPayment` and nothing +else. The `multiPartPayment` exposes an interface to query information about the MPP that the +Session has computed, therefore all of the `multiPartPayment_*` functions, except for the +destructor, do not modify the state of the MPP, the `multiPartPayment` is read-only. +``` +typedef struct multiPartPayment multiPartPayment; +``` + +Once a `multiPartPayment` object is no longer used, it can be freed by calling +`multiPartPayment_delete`. This function never fails. +``` +void multiPartPayment_delete(multiPartPayment* mpp); +``` + +`multiPartPayment` can be queried to get the information about the MPP information they encode. +To query the amount of satoshis this payment sends one calls `multiPartPayment_amount`. +This function never fails. +``` +int64_t multiPartPayment_amount(const multiPartPayment* mpp); +``` + +To query the source and the destination of the payment the functions +`multiPartPayment_source` and `multiPartPayment_destination` are called respectively. +These functions never fail. +``` +nodeID_t multiPartPayment_source(const multiPartPayment* mpp); +``` +``` +nodeID_t multiPartPayment_destination(const multiPartPayment* mpp); +``` + +Each MPP consist of a list of paths in the channels/nodes graph. +By calling `multiPartPayment_numParts` one can query for the number of paths that constitute the +MPP. +This functions never fails. +``` +uint64_t multiPartPayment_numParts(const multiPartPayment* mpp); +``` + +By calling `multiPartPayment_pathLengths` one can query the length of the individual paths that make +the MPP. The pointer `lengths` must be previously allocated to hold at least as many paths as the +MPP. +This function never fails. +``` +void multiPartPayment_pathLengths(const multiPartPayment* mpp, + uint64_t* lengths); +``` + +Similar to the `multiPartPayment_pathLengths`, the function `multiPartPayment_pathValues` queries +the MPP to obtain the amount of satoshis commited to each individual path. +The pointer `values` must be previously allocated to hold at least as many paths as the +MPP. +This function never fails. +``` +void multiPartPayment_pathValues(const multiPartPayment* mpp, + int64_t* values); +``` + +The function `multiPartPayment_pathChannels` writes the list of channels for every path in the MPP +to the input array `short_channel_id`, and for each of these channels a value of the `orientation` +is set either to `CHANNEL_FORWARD` or `CHANNEL_BACKWARD` to indicate the direction the liquidity is +transfered. If the satoshis are transfered from a node with lowest lexicographical order to the +other then the `orientation` is set to `CHANNEL_FORWARD`, otherwise it is `CHANNEL_BACKWARD`. +This function never fails. +``` +void multiPartPayment_pathChannels(const multiPartPayment* mpp, + int64_t* short_channel_id, + char* orientation); +``` + +## Example + +``` +// TODO +``` diff --git a/python-documentation.md b/python-documentation.md new file mode 100644 index 0000000..f5d384c --- /dev/null +++ b/python-documentation.md @@ -0,0 +1,36 @@ +# Python-API for Pickhardt Payments + +## Overview + +This document describes *Pickhardt payments* Python library API. + +## Example + +Below there is an example use case of the Python API + +``` +# perform a payment from node A to node B, of amt satoshis + +import pickhardtpayments + +gossip = get_gossip() # whatever resource that provides a list of the channels in LN +oracle = new_oracle() # a simulation of LN or a direct connection to it + +pay_mpp = pickhardtpayments.paymentSession() + +for channel in gossip: + pay_mpp.add_channel(channel) + +res_amt = amt +while res_amt>0: + mpp = pay_mpp.optimizedPayment(A,B,res_amt) + + for path,value in mpp: + success,bad_channel = oracle.send_onion(path,value) + + if success: + res_amt -= value + pay_mpp.updateSuccess(path,value) + else: + pay_mpp.updateFailure(path,value,bad_channel) +``` From ce304f1d5c909ef235f03e1c78e0b6cc3409b253 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Wed, 22 Jun 2022 07:05:24 +0200 Subject: [PATCH 02/20] change variable name in python example --- python-documentation.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python-documentation.md b/python-documentation.md index f5d384c..808c71f 100644 --- a/python-documentation.md +++ b/python-documentation.md @@ -16,21 +16,21 @@ import pickhardtpayments gossip = get_gossip() # whatever resource that provides a list of the channels in LN oracle = new_oracle() # a simulation of LN or a direct connection to it -pay_mpp = pickhardtpayments.paymentSession() +session = pickhardtpayments.paymentSession() for channel in gossip: - pay_mpp.add_channel(channel) + session.add_channel(channel) res_amt = amt while res_amt>0: - mpp = pay_mpp.optimizedPayment(A,B,res_amt) + mpp = session.optimizedPayment(A,B,res_amt) for path,value in mpp: success,bad_channel = oracle.send_onion(path,value) if success: res_amt -= value - pay_mpp.updateSuccess(path,value) + session.updateSuccess(path,value) else: - pay_mpp.updateFailure(path,value,bad_channel) + session.updateFailure(path,value,bad_channel) ``` From a3c388c54ece4cc833a1cd74ceda0ed4b8574dbc Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Tue, 9 Aug 2022 10:04:59 +0200 Subject: [PATCH 03/20] python MinCostFlow into SyncSimulatedPaymentSession --- .../SyncSimulatedPaymentSession.py | 22 +++++++++---------- 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index 1be5c23..41c6b90 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -1,7 +1,7 @@ from .UncertaintyNetwork import UncertaintyNetwork from .OracleLightningNetwork import OracleLightningNetwork -from ortools.graph import pywrapgraph +from .MinCostFlow import MCFNetwork from typing import List @@ -56,7 +56,7 @@ def _prepare_mcf_solver(self, src, dest, amt: int = 1, mu: int = 100_000_000, ba returns the instantiated min_cost_flow object from the google OR-lib that contains the piecewise linearized problem """ - self._min_cost_flow = pywrapgraph.SimpleMinCostFlow() + self._min_cost_flow = MCFNetwork() self._arc_to_channel = {} for s, d, channel in self._uncertainty_network.network.edges(data="channel"): @@ -73,10 +73,10 @@ def _prepare_mcf_solver(self, src, dest, amt: int = 1, mu: int = 100_000_000, ba cnt = 0 # QUANTIZATION): for capacity, cost in channel.get_piecewise_linearized_costs(mu=mu): - index = self._min_cost_flow.AddArcWithCapacityAndUnitCost(self._mcf_id[s], - self._mcf_id[d], - capacity, - cost) + index = self._min_cost_flow.AddArc(self._mcf_id[s], + self._mcf_id[d], + capacity, + cost) self._arc_to_channel[index] = (s, d, channel, 0) if self._prune_network and cnt > 1: break @@ -197,12 +197,10 @@ def _generate_candidate_paths(self, src, dest, amt, mu: int = 100_000_000, base: start = time.time() #print("solving mcf...") - status = self._min_cost_flow.Solve() - - if status != self._min_cost_flow.OPTIMAL: - print('There was an issue with the min cost flow input.') - print(f'Status: {status}') - exit(1) + try: + self._min_cost_flow.Solve() + except: + raise BaseException('There was an issue with the min cost flow input.') paths = self._disect_flow_to_paths(src, dest) end = time.time() From 306492aae536c1b66e15f49b45f9372c18dc4eb2 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Tue, 9 Aug 2022 10:22:21 +0200 Subject: [PATCH 04/20] added MinCostFlow as submodule --- .gitmodules | 4 + subprojects/MinCostFlow/LICENSE | 190 +++++++ .../MinCostFlow/benchmark/benchmark-mcf.cpp | 299 ++++++++++ subprojects/MinCostFlow/benchmark/gen.py | 105 ++++ subprojects/MinCostFlow/benchmark/meson.build | 4 + subprojects/MinCostFlow/benchmark/run.sh | 20 + subprojects/MinCostFlow/c-api.md | 109 ++++ subprojects/MinCostFlow/doc/biblio.bib | 21 + subprojects/MinCostFlow/doc/info.tex | 56 ++ subprojects/MinCostFlow/doc/intro.tex | 36 ++ subprojects/MinCostFlow/doc/main.tex | 47 ++ subprojects/MinCostFlow/doc/makefile | 17 + subprojects/MinCostFlow/doc/proposal.tex | 47 ++ subprojects/MinCostFlow/doc/schedule.tex | 48 ++ subprojects/MinCostFlow/documentation.md | 78 +++ .../MinCostFlow/examples/kattis-maxflow.cpp | 72 +++ .../examples/kattis-mincostmaxflow.cpp | 76 +++ .../examples/kattis-shortestpath1.cpp | 64 +++ subprojects/MinCostFlow/examples/meson.build | 8 + subprojects/MinCostFlow/include/meson.build | 2 + .../MinCostFlow/include/mincostflow/graph.hpp | 412 ++++++++++++++ .../include/mincostflow/maxflow.hpp | 289 ++++++++++ .../include/mincostflow/meson.build | 1 + .../include/mincostflow/mincostflow.hpp | 479 ++++++++++++++++ .../include/mincostflow/scope_guard.hpp | 45 ++ .../include/mincostflow/shortestpath.hpp | 518 ++++++++++++++++++ .../include/mincostflow/vectorized_map.hpp | 259 +++++++++ subprojects/MinCostFlow/meson.build | 21 + subprojects/MinCostFlow/readme.md | 18 + subprojects/MinCostFlow/test/max-flow.cpp | 90 +++ subprojects/MinCostFlow/test/meson.build | 15 + .../MinCostFlow/test/min-cost-flow.cpp | 274 +++++++++ subprojects/MinCostFlow/test/short-path.cpp | 86 +++ subprojects/MinCostFlow/test/vec-map.cpp | 41 ++ 34 files changed, 3851 insertions(+) create mode 100644 .gitmodules create mode 100644 subprojects/MinCostFlow/LICENSE create mode 100644 subprojects/MinCostFlow/benchmark/benchmark-mcf.cpp create mode 100644 subprojects/MinCostFlow/benchmark/gen.py create mode 100644 subprojects/MinCostFlow/benchmark/meson.build create mode 100644 subprojects/MinCostFlow/benchmark/run.sh create mode 100644 subprojects/MinCostFlow/c-api.md create mode 100644 subprojects/MinCostFlow/doc/biblio.bib create mode 100644 subprojects/MinCostFlow/doc/info.tex create mode 100644 subprojects/MinCostFlow/doc/intro.tex create mode 100644 subprojects/MinCostFlow/doc/main.tex create mode 100644 subprojects/MinCostFlow/doc/makefile create mode 100644 subprojects/MinCostFlow/doc/proposal.tex create mode 100644 subprojects/MinCostFlow/doc/schedule.tex create mode 100644 subprojects/MinCostFlow/documentation.md create mode 100644 subprojects/MinCostFlow/examples/kattis-maxflow.cpp create mode 100644 subprojects/MinCostFlow/examples/kattis-mincostmaxflow.cpp create mode 100644 subprojects/MinCostFlow/examples/kattis-shortestpath1.cpp create mode 100644 subprojects/MinCostFlow/examples/meson.build create mode 100644 subprojects/MinCostFlow/include/meson.build create mode 100644 subprojects/MinCostFlow/include/mincostflow/graph.hpp create mode 100644 subprojects/MinCostFlow/include/mincostflow/maxflow.hpp create mode 100644 subprojects/MinCostFlow/include/mincostflow/meson.build create mode 100644 subprojects/MinCostFlow/include/mincostflow/mincostflow.hpp create mode 100644 subprojects/MinCostFlow/include/mincostflow/scope_guard.hpp create mode 100644 subprojects/MinCostFlow/include/mincostflow/shortestpath.hpp create mode 100644 subprojects/MinCostFlow/include/mincostflow/vectorized_map.hpp create mode 100644 subprojects/MinCostFlow/meson.build create mode 100644 subprojects/MinCostFlow/readme.md create mode 100644 subprojects/MinCostFlow/test/max-flow.cpp create mode 100644 subprojects/MinCostFlow/test/meson.build create mode 100644 subprojects/MinCostFlow/test/min-cost-flow.cpp create mode 100644 subprojects/MinCostFlow/test/short-path.cpp create mode 100644 subprojects/MinCostFlow/test/vec-map.cpp diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..e799283 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,4 @@ +[submodule "subprojects/MinCostFlow"] + path = subprojects/MinCostFlow + url = https://github.com/Lagrang3/mincostflow.git + branch = master diff --git a/subprojects/MinCostFlow/LICENSE b/subprojects/MinCostFlow/LICENSE new file mode 100644 index 0000000..5c67b5f --- /dev/null +++ b/subprojects/MinCostFlow/LICENSE @@ -0,0 +1,190 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + Copyright 2022 Eduardo Quintana Miranda + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/subprojects/MinCostFlow/benchmark/benchmark-mcf.cpp b/subprojects/MinCostFlow/benchmark/benchmark-mcf.cpp new file mode 100644 index 0000000..def00a7 --- /dev/null +++ b/subprojects/MinCostFlow/benchmark/benchmark-mcf.cpp @@ -0,0 +1,299 @@ +#include +#include +#include +#include +#include + +typedef long long value_type; +typedef int nodeID_type; +typedef int arcID_type; +typedef ln::digraph graph_type; + +// typedef ln::maxflow_augmenting_path maxflow_t; +// typedef ln::maxflow_augmenting_path maxflow_t; +// typedef ln::maxflow_scaling maxflow_t; +// typedef ln::maxflow_scaling maxflow_t; +// typedef ln::maxflow_preflow maxflow_t; + +//typedef ln::mincostflow_EdmondsKarp< +// value_type,ln::shortestPath_FIFO> mincostflow_t; // 2.18s +// typedef ln::mincostflow_EdmondsKarp< +// value_type,ln::shortestPath_BellmanFord> mincostflow_t; // TLE +// typedef ln::mincostflow_EdmondsKarp mincostflow_t; // expected failure + +// typedef ln::mincostflow_PrimalDual< +// ln::shortestPath_FIFO,maxflow_t> mincostflow_t; // 2.54s, 3.02s, 2.75s, 3.18s, TLE +// typedef ln::mincostflow_PrimalDual< +// ln::shortestPath_BellmanFord,maxflow_t> mincostflow_t; // TLE +//typedef ln::mincostflow_PrimalDual< +// ln::shortestPath_Dijkstra,maxflow_t> mincostflow_t; +// 1.34s, 1.76s, 1.52s, 2.42s, TLE + +// typedef ln::mincostflow_capacityScaling< +// ln::shortestPath_Dijkstra,maxflow_t> mincostflow_t; +// TLE + +struct integer_graph +{ + graph_type g; + std::vector capacity; + std::vector weight; + + integer_graph(int N_nodes) + { + for(nodeID_type i=0;i& cap) + { + capacity.resize(g.max_num_arcs()); + for(arcID_type i=0;i& cost) + { + weight.resize(g.max_num_arcs()); + for(arcID_type i=0;i solve_ortools( + const int N_nodes, + nodeID_type S, nodeID_type T, + const std::vector>& edges, + const std::vector& cap, + const std::vector& cost, + const std::string tname) +{ + auto start = std::chrono::high_resolution_clock::now(); + operations_research::SimpleMaxFlow max_flow; + operations_research::SimpleMinCostFlow mincost_flow; + + for(int i=0;i(stop-start) .count() + << std::endl; + + return {Flow,Cost}; +} + + +template +void solve(integer_graph& G, + const std::vector& capacity, + const std::vector& cost, + nodeID_type S, + nodeID_type T, + std::string tname) +{ + G.set_capacity(capacity); + G.set_cost(cost); + + auto& graph = G.g; + + auto start = std::chrono::high_resolution_clock::now(); + mincostflow_t f; + f.solve(graph,graph.get_node(S),graph.get_node(T),G.weight,G.capacity); + auto stop = std::chrono::high_resolution_clock::now(); + + std::cout << tname << " " << + std::chrono::duration_cast(stop-start) .count() + << std::endl; +} + +std::pair +check_constraints( + const integer_graph& G, + const std::vector>& ed_list, + const std::vector& capacity, + const std::vector& cost, + const nodeID_type S, + const nodeID_type T, + const int N_nodes, + const int N_edges) +{ + // check capacity constraint + for(int e=0;e=f && f>=0); + } + + // check balances + std::vector balance(N_nodes,0); + + for(int e=0;e=0); + assert(balance[S] == -balance[T]); + + value_type Cost = 0; + for(int e=0;e, + ln::maxflow_augmenting_path + > + >(G,capacity,weight,S,T,"Primal-dual"); + auto [flow,cost] = check_constraints(G,ed_list,capacity,weight,S,T,N,M); + + assert(flow_0==flow && cost_0==cost); + } + + // { + // solve< + // ln::mincostflow_capacityScaling< + // ln::shortestPath_Dijkstra, + // ln::maxflow_augmenting_path + // > + // >(G,capacity,weight,S,T,"Capacity-scaling"); + // auto [flow,cost] = check_constraints(G,ed_list,capacity,weight,S,T,N,M); + // assert(flow_0==flow && cost_0==cost); + // } + { + solve< + ln::mincostflow_costScaling< + ln::maxflow_augmenting_path + > + >(G,capacity,weight,S,T,"Cost-scaling"); + auto [flow,cost] = check_constraints(G,ed_list,capacity,weight,S,T,N,M); + assert(flow_0==flow && cost_0==cost); + } + { + auto [flow,cost] = solve_ortools(N,S,T,ed_list,capacity,weight,"Ortools"); + assert(flow_0==flow && cost_0==cost); + } + return 0; +} + diff --git a/subprojects/MinCostFlow/benchmark/gen.py b/subprojects/MinCostFlow/benchmark/gen.py new file mode 100644 index 0000000..facb244 --- /dev/null +++ b/subprojects/MinCostFlow/benchmark/gen.py @@ -0,0 +1,105 @@ +import networkx +import random +import sys,os +import tempfile +import numpy +import subprocess +import matplotlib.pyplot as plt +import json + +def generate(n,m,mcap,mwei,f): + S=0 + T=1 + g=networkx.gnm_random_graph(n,m,directed=True) + print(n,m,S,T,file=f) + for a,b in g.edges(): + print(a,b,random.choice(range(mcap)),random.choice(range(mwei)),file=f) + +def plot_data(n_array,series,outname): + fig = plt.figure() + ax = fig.subplots() + for x in series: + ax.loglog(n_array,series[x],'-o',label=x) + ax.legend() + ax.set_xlabel('N nodes') + ax.set_ylabel('Time (micro seconds)') + ax.grid() + fig.savefig(outname) + plt.show() + +def save_data(n_array,series,outname): + alldata={'N' : n_array} + alldata['series']=dict() + for name in series: + alldata['series'][name]=[ int(x) for x in series[name] ] + json.dump(alldata,open(outname,'w')) + +def read_data(filename): + alldata = json.load(open(filename,'r')) + return alldata['N'],alldata['series'] + +def cat(fname): + f = open(fname,'r') + for l in f: + print(l,end='') + +class TmpFile(): + def __init__(self): + self.name = tempfile.mktemp() + def __enter__(self): + return self + def __exit__(self,exc_type,exc_value,exc_traceback): + os.remove(self.name) + +def run_benchmark(): + n_array = [ 2**i for i in range(7,15)] + m_array = [ int(n*7.5) for n in n_array ] + cost_arr = [ 200 ]*len(n_array) + cap_arr = [ 200 ]*len(n_array) + rep_arr = [20]*len(n_array) + + names = [ + 'Ortools', + 'Edmonds-Karp', + 'Primal-dual', + # 'Capacity-scaling', + 'Cost-scaling', + ] + + series = dict() + for x in names: + series[x] = numpy.zeros(len(n_array)) + + for i in range(len(n_array)): + n,m,cap,cost = n_array[i],m_array[i],cap_arr[i],cost_arr[i] + print("N nodes = ",n) + for j in range(rep_arr[i]): + print(" rep = ",j) + with TmpFile() as tmp: + with open(tmp.name,'w') as fin: + generate(n,m,cap,cost,fin) + # cat(tmp.name) + p = subprocess.run( + ['./benchmark/benchmark-mcf'], + stdin=open(tmp.name,'r'),capture_output=True) + + data = p.stdout.decode().split('\n') + if p.returncode!=0: + print(p.stderr.decode()) + print("on test case") + cat(tmp.name) + raise "benchmark failed" + for d in data: + try: + name,val = d.split() + val = int(val) + series[name][i] = max(series[name][i],val) + except: + pass + return n_array, series + +if __name__ == "__main__": + #n_array,series = read_data('latest.json') + n_array,series = run_benchmark() + save_data(n_array,series,'latest.json') + plot_data(n_array,series,'latest.png') diff --git a/subprojects/MinCostFlow/benchmark/meson.build b/subprojects/MinCostFlow/benchmark/meson.build new file mode 100644 index 0000000..59965f0 --- /dev/null +++ b/subprojects/MinCostFlow/benchmark/meson.build @@ -0,0 +1,4 @@ +ortools = dependency('ortools') + +executable('benchmark-mcf','benchmark-mcf.cpp', + dependencies: [dep_mincostflow,ortools]) diff --git a/subprojects/MinCostFlow/benchmark/run.sh b/subprojects/MinCostFlow/benchmark/run.sh new file mode 100644 index 0000000..b2f23bf --- /dev/null +++ b/subprojects/MinCostFlow/benchmark/run.sh @@ -0,0 +1,20 @@ +tmp=`mktemp` +exe=benchmark/benchmark-mcf +gen=../benchmark/gen.py + +# lightning +N=2000 +M=150000 + +N=800 +M=6000 + +Mcap=200 +Mcost=200 + +for i in `seq 10`; do + echo "" + python $gen $N $M $Mcap $Mcost > $tmp + $exe < $tmp + +done diff --git a/subprojects/MinCostFlow/c-api.md b/subprojects/MinCostFlow/c-api.md new file mode 100644 index 0000000..ee002ef --- /dev/null +++ b/subprojects/MinCostFlow/c-api.md @@ -0,0 +1,109 @@ +Graph representation +=== +The problem graph is represented by the type `lngraph`. +This data structure assumes that the graph is directed and there can be multiple edges connecting +the same pair of vertices. + +The graph API is edge centered. Edges can be added, modified or removed. +Edges are refered to by their unique id, which is 64bit number that the API assigns internally. +This assignment is unspecified. + +There is no need now to define an API to do the same thing for nodes, they are just added +automatically and possibly never removed from the graph. + +The constant values `lngraph_SUCCESS` and `lngraph_FAIL` are generally used as return values to +indicate the success or failure in a function call. + +Memory allocation +--- +Memory is allocated in the heap for a `lngraph` by calling +``` +lngraph* lngraph_new(); +``` +This function returns a pointer to a newly created `lngraph` in a valid state. +The returned value is `NULL` if the allocation or the initialization fails. + +Memory is freed by calling +``` +int lngraph_free(lngraph* g); +``` +No fail. You may assume this always return `lngraph_SUCCESS`. + +Getters +--- +Returns the number of edges. +``` +int lngraph_numEdges(lngraph* g); +``` +No fail. + + +Gets the list of edges sorted by their ids and their capacity and cost properties. +`edgeID`, `capacity` and `cost` need to be allocated to hold at least as many elements as edges in +the graph. +``` +int lngraph_edges(lngraph* g, int64_t * edgeID, int64_t *capacity, int64_t * cost); +``` +No fail. You may assume this always return `lngraph_SUCCESS`. + + +Setters +--- +Edges are added to the graph by calling the function +``` +int64_t lngraph_addEdge(lngraph* g, const char* nodeA_id, const char* nodeB_id, int64_t capacity, +int64_t cost); +``` +The function returns a `int64_t` that corresponds to a unique identifier for the edge in the +database. +Nodes are added automatically if needed everytime `lngraph_addEdge` is called. +Possible errors: fails to add a new edge, the return value in that case is `lngraph_NOID`. +There will never be an edge whose id corresponds to `lngraph_NOID`. + +Edges can be removed by calling +``` +int lngraph_rmEdge(lngraph* g, int64_t edgeUniqueID); +``` +No fail. You may assume this always return `lngraph_SUCCESS`. + + +Min. Cost Flow solver +=== + +Node excess. A possitive value of the excess makes a node into a *source node* +and a negative value makes it a *sink node*. +``` +int lngraph_setExcess(lngraph* g, const char* nodeID, int64_t excess); +``` +If the node does not exist then it will be added automatically. +Possible errors: a new node is failed to be added, in that case the returned value is `lngraph_FAIL`. +Otherwise `lngraph_SUCCESS` is returned. + + +Sets the excess of every known node to 0. +``` +int lngraph_resetExcess(lngraph* g); +``` +No fail. You may assume this always return `lngraph_SUCCESS`. + +Solve the MCF problem with the excess previously set and returns the status. +``` +int lngraph_mincostflow(lngraph * g,int64_t * flow); +``` +The `flow` pointer points to a an allocated array of at least M elements, where M is the current +number of edges. Edges are ordered internally according to their ids and the `flow` value +corresponds to edges in that order. +Possible errors: it was not possible to find a feasible flow that satisfy the problem constraints, +the return value is `lngraph_FAIL`. +Otherwise `lngraph_SUCCESS` is returned. + +Solve the MCF problem where there is only a source and a sink node. +This is equivalent to setting the excess of the `source` to the value `excess` +and the excess of the `sink` to `-excess` and then calling `lngraph_mincostflow`. +``` +int lngraph_mincostflow_singleSinkSource(lngraph * g, + const char* source, + const char* sink, + int64_t excess, + int64_t * flow); +``` diff --git a/subprojects/MinCostFlow/doc/biblio.bib b/subprojects/MinCostFlow/doc/biblio.bib new file mode 100644 index 0000000..a3ab84a --- /dev/null +++ b/subprojects/MinCostFlow/doc/biblio.bib @@ -0,0 +1,21 @@ +@book{ahuja1993network, + title={Network Flows: Theory, Algorithms, and Applications}, + author={Ahuja, R.K. and Magnanti, T.L. and Orlin, J.B.}, + isbn={9780136175490}, + lccn={lc92026702}, + url={https://books.google.it/books?id=WnZRAAAAMAAJ}, + year={1993}, + publisher={Prentice Hall} +} + +@misc{pickhardt21, + doi = {10.48550/ARXIV.2107.05322}, + url = {https://arxiv.org/abs/2107.05322}, + author = {Pickhardt, Rene and Richter, Stefan}, + keywords = {Networking and Internet Architecture (cs.NI), FOS: Computer and information sciences, FOS: Computer and information sciences}, + title = {Optimally Reliable \& Cheap Payment Flows on the Lightning Network}, + publisher = {arXiv}, + year = {2021}, + copyright = {Creative Commons Attribution Share Alike 4.0 International} +} + diff --git a/subprojects/MinCostFlow/doc/info.tex b/subprojects/MinCostFlow/doc/info.tex new file mode 100644 index 0000000..6df4664 --- /dev/null +++ b/subprojects/MinCostFlow/doc/info.tex @@ -0,0 +1,56 @@ +\section{Biographical Information} + +\subsection*{Personal details} + +\begin{tabular}{ll} + Name:& Eduardo Quintana Miranda\\ + Affiliation:& University of Trieste, Astronomical Observatory of Trieste\\ + Course:& Astrophysics \\ + Degree Program:& PhD \\ + Country:& Italy\\ + E-mail:& \url{eduardo.quintana@pm.me} \\ + Github:& \url{github.com/Lagrang3} \\ + Discord:& \texttt{eduardo-qm} +\end{tabular} + +\subsection*{Educational background.} + +\begin{itemize} + \item 2008--2013. BSc. Nuclear Physics, University of Havana, + \item 2015--2018. MSc. Theoretical Physics, University of Trieste, + \item + 2018--2019. Master in High Performance Computing, International School of + Advanced Studies (SISSA), Trieste, + \item 2019--present. PhD in astrophysics, + University of Trieste. +\end{itemize} + +\subsection*{Programming background.} +My main programming language is C++, but I am also profiecient in C and python. + +I started programming back in 2009 during my 2nd year of BSc. The main drivers +were physics simulations and programming competitions. From 2009 until 2013 +I've participated intensively in those kind of competitions, +the most important were the +ACM-ICPC\footnote{\url{icpc.global}} competitions held every year, my team +managed to classify many times to the regional phase, which in our geographical +context where named \emph{Caribbean Finals of the ICPC}. Until this date I +participate once in a while on \url{codeforces.com} rounds by the username +\texttt{Lagrang3}\footnote{\url{https://codeforces.com/profile/Lagrang3}}. +I have a fair knowledge of a variety of classical algorithms somewhat seasoned +by my participation on programming competitions. I am familiar with shortest +path, max-flow and linear min-cost flow graph problems and textbook algorithms +for solving them. + +I've acquired some professional programming skills during the Master in High +Performance Computing. Advanced C++, Python, bash, git, and parallel programming +in OpenMP, MPI and Cuda where among the topics taught in that curriculum. + +In 2021 I've succesfully completed a Google Summer of Code project for the +proposal of FFT utilities in Boost Math +library\footnote{\url{https://github.com/BoostGSoC21/math-fft-report/releases/download/v1.1/gsoc-report.pdf}}. + +Today, I am doing a PhD in astrophysics working on the developement of a +simulation code +to study the effects of general relativity in the formation of cosmological +structures.\footnote{\url{https://github.com/Lagrang3/gevolution-1.2/tree/gev-api}} diff --git a/subprojects/MinCostFlow/doc/intro.tex b/subprojects/MinCostFlow/doc/intro.tex new file mode 100644 index 0000000..6c7d67b --- /dev/null +++ b/subprojects/MinCostFlow/doc/intro.tex @@ -0,0 +1,36 @@ +\section{Introduction} + +Payments in the Lightning Network are performed along shortest path on the +network that try to minimize the transfer fees. +Because the participating nodes have incomplete information about the balance of +the individual channels, it is not always possible to perform the payments along +the shortest path, and the back-end of lightning has to work out alternative +routes until one is found with enough liquidity to process the payment. + +\cite{pickhardt21} proposes an alternative algorithm for payments, by +introducing the following novel ideas: +\begin{enumerate} + \item Each payment can be routed through a combination of different paths, + in what is called a \emph{Mult-Part Payment} (MPP). + The transfer of money on each channel is bounded by the channel's capacity + (in one specific direction) and the sum of the inbound and outbound payments + for the nodes satisfy balance constraints, the MPP problem can be solved + with the algorithms that find maximum flows on a directed network. + \item The arcs of the network can be characterized by a probability + distribution function for the realization of the transfer (the value of the + flow on the arc that can be processed). + The probability of success of the MPP becomes the multiplication of the + probability of success of every non-trivial flow on the network. + \item By defining a new function of the MPP flow, + which is the negative logarithm of that probability of success, + one has that the probability of an MPP payment is maximized when this new + cost function is minimized. And since the logarithm of the product is the + sum of logarithms, this cost function is separable in the sum of the cost + for the individual costs of the flow through the arcs of the network. + Thus the problem of maximizing the probability of success of a MPP payment + becomes a Min Cost Flow problem with non-linear costs on arcs. + \item By assuming some constraint in the probability function of the arcs, + one can reduce to a particular class of problems for which the cost of the + flow is a convex function. For this class of problem efficient polynomial + algorithms are known. +\end{enumerate} diff --git a/subprojects/MinCostFlow/doc/main.tex b/subprojects/MinCostFlow/doc/main.tex new file mode 100644 index 0000000..27a4e29 --- /dev/null +++ b/subprojects/MinCostFlow/doc/main.tex @@ -0,0 +1,47 @@ +\documentclass[11pt,a4paper]{article} + +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{amsthm} +\usepackage[a4paper,margin=1in]{geometry} +\usepackage{graphicx} +\usepackage{subcaption} +\usepackage{hyperref} +\usepackage{listings} +\usepackage{bold-extra} + +\hypersetup{ + colorlinks=false, + hidelinks=true +} +\bibliographystyle{plain} % FIXME: styles are not found + +\newcommand{\Ring}{\mathcal{R}} +\newcommand{\Order}{\mathcal{O}} +\newcommand{\Complex}{\mathbb{C}} +\newcommand{\Real}{\mathbb{R}} +\newcommand{\Integer}{\mathbb{Z}} + +\newcommand{\comment}[1]{\textrm{\textit{#1}}} +\newcommand{\function}[1]{\textbf{#1}} +\newcommand{\variable}[1]{\textit{#1}} + +\title{Summer of Bitcoin 2022:\\ Efficient Min Convex Cost Flow Solver for Lightning Network Payment Delivery} +\author{Eduardo Quintana Miranda} + + +\begin{document} + +\maketitle + +%\tableofcontents +\input intro +\input proposal +\input schedule + +\bibliography{biblio} + +\appendix +\input info + +\end{document} diff --git a/subprojects/MinCostFlow/doc/makefile b/subprojects/MinCostFlow/doc/makefile new file mode 100644 index 0000000..084ab10 --- /dev/null +++ b/subprojects/MinCostFlow/doc/makefile @@ -0,0 +1,17 @@ +VPATH:=images +TARGETS:=main.pdf +LATEXEC:=rubber --pdf main.tex +IMAGES:= +TEX:=main.tex biblio.bib info.tex intro.tex proposal.tex schedule.tex + +export + +default: $(TARGETS) + +main.pdf : $(TEX) $(IMAGES) + -$(LATEXEC) + +clean: + -$(LATEXEC) --clean + +.PHONY: clean default diff --git a/subprojects/MinCostFlow/doc/proposal.tex b/subprojects/MinCostFlow/doc/proposal.tex new file mode 100644 index 0000000..e58b297 --- /dev/null +++ b/subprojects/MinCostFlow/doc/proposal.tex @@ -0,0 +1,47 @@ +\section{Proposal} + +% Based on the description of the problem in the introduction +This Summer of Bitcoin project, proposes the design and implementation of a C++ +mini-library with no external dependencies for solving the min cost flow problem +(MCF) with convex costs that can be used to optimize the probability of success +of multipart payments in the Lightning Network. + +Some time of the project will be allocated to study of the state of the art of +the problem. At the present we are aware of two algorithms that can be used to +solve the MCF with convex costs. +One that transforms approximates the cost function to a piecewise linear +function and splits the arc into several linear cost arcs accordingly (section +14.4 of \cite{ahuja1993network}). +And the so called \emph{Capacity} +algorithm described in section 14.5 of \cite{ahuja1993network}, +which is a generalization of a \emph{Excess Scaling} solver for the linear +MCF problem. +During the course of this project duration, we will evaluate which algorithm +will be best suited for the problem domain. + +Parallelization of the algorithms will be considered and discussed during the +development of the back-end. +We will propose to use standard library tools for that purpose +which are available in C++17. + +For the preparation of this Summer of Bitcoin project proposal, we have been +working on a proof-of-concept MCF C++ library\footnote{% +\url{https://github.com/Lagrang3/mincostflow}} publicly available on Github. +The API of this MCF library permits the representation of a directed graph +\texttt{digraph} in a space efficient data structure. +It also allows for the computation +shortest paths on a weighted network using Dijsktra's and Breadth-First-Search +algorithms, which serve as building blocks for flow solvers. +There are two different back-ends for the solution of the Max Flow problem: +using Edmond-Karp's \emph{Augmenting Paths} and Dinic's \emph{Push-Relabel}. +Also these constitute building blocks for the Min Cost Flow solvers. +The library contains a MCF solver based on the \emph{Successive Shortest Path} +algorithm presented in section 9.7 on \cite{ahuja1993network}. + +In the library's respository we provide examples of applications that we have also +used as test for correctedness and efficiency of our implementation. +These examples actually solve two problems from the online judge Kattis: +\texttt{maxflow}\footnote{\url{https://open.kattis.com/problems/maxflow}} +and +\texttt{mincostmaxflow}\footnote{\url{https://open.kattis.com/problems/mincostmaxflow}}, +and they do pass the test cases within the time and memory constraints. diff --git a/subprojects/MinCostFlow/doc/schedule.tex b/subprojects/MinCostFlow/doc/schedule.tex new file mode 100644 index 0000000..329c540 --- /dev/null +++ b/subprojects/MinCostFlow/doc/schedule.tex @@ -0,0 +1,48 @@ +\section{Timeline} + + \subsection*{Before May 9} + Meet the mentors, discuss with them the proposed timeline, + deliverables and share ideas. + + \subsection*{May 9 -- May 23} + Search in the literature for algorithms that solve the min-cost flow problem + with convex cost functions. + Follow the bitcoin and lightning seminars. + + \subsection*{May 23 -- June 6} + Design of a standalone C++ library for solving the min-cost flow problem. + Implementation of a textbook efficient solution to the linear min-cost flow + problem, the \emph{Excess Scaling} algorithm described in + \cite{ahuja1993network} section 7.9. + This linear MCF solver can be used to solve convex cost problems with a + tuneable linear approximation by decomposing the non-linear cost arcs into + several linear cost arcs. + Set up a simple set of unit tests for correctedness and CI for the library + repository. + + \subsection*{June 6 -- June 20} + Use the code to reproduce the results + of Pickhardt's \emph{Payment Simulation} Jupyter Notebook. + Implementation of a textbook solution: the capacity scaling algorithm described in + \cite{ahuja1993network} section 14.5. This algorithm is a generalization of + the \emph{Excess Scaling}, hence it shouldn't be too difficult to adapt + the linear MCF already implemented to the general convex cost case. + + \subsection*{June 20 -- July 4} + Benchmark the canditate solutions. + Optimize the code to find a good tradeoff of runtimes and probability of + success. + Parallelization of the time critical routines. + + \subsection*{July 4 -- July 18} + Implementation of a small C-lightning plugin to demonstrate the use of the + library. + + \subsection*{July 18 -- August 1} + Document the library and buffer time for unfinished issues. + + \subsection*{August 1 -- August 15} + Preparation of the final report and draft the research article. + + \subsection*{August 15 -- August 22} + Submission of the final project report. diff --git a/subprojects/MinCostFlow/documentation.md b/subprojects/MinCostFlow/documentation.md new file mode 100644 index 0000000..e143c6b --- /dev/null +++ b/subprojects/MinCostFlow/documentation.md @@ -0,0 +1,78 @@ +Graph representation +=== +``` + class digraph; +``` + +Path solvers +=== + +Finds a path in a directed graph with un-weighted edges. +``` + class pathSearch_BFS; +``` + +Finds a path in a directed graph with un-weighted edges. +This label-relabel search algorithm has meaninful state. +Figure 7.6 of Ahuja 93. +``` + class pathSearch_labeling; +``` + +Pseudo-polynomial generic path optimization. +``` + class shortestPath_FIFO; +``` + +Bellman-Ford path optimization +``` + class shortestPath_BellmanFord; +``` + +Dijkstra path optimization. +``` + class shortestPath_Dijkstra; +``` + +Max-flow +=== + +Generic augmenting path algorithm, template on the path finder algorithm. +Ahuja figure 6.12. +``` + template + class maxflow_augmenting_path; +``` + +Capacity scaling algorithm template on the path finder. +Ahuja figure 7.3. +``` + template + class maxflow_scaling; +``` + +Preflow-push algorithm. +Ahuja figure 7.12. +``` + class maxflow_preflow; +``` + +Min-Cost-Flow +=== +Min-cost-max-flow based on Edmonds-Karp augmenting path algorithm, it greedily +searches for the smallest cost routes. The path optimizer could be any shortest +path algorithm that deals with negative weights. +``` + template + class mincostflow_EdmondsKarp; +``` + +Min-cost-max-flow Primal Dual algorithm. +It uses a potential function to set the edges costs to zero and it pushes flow +along cost-zero edges. +It is a template on a path optimizer and a maxflow engines. +The path optimizer could be any shortest path algorithm including Dijkstra. +``` + template + class mincostflow_PrimalDual; +``` diff --git a/subprojects/MinCostFlow/examples/kattis-maxflow.cpp b/subprojects/MinCostFlow/examples/kattis-maxflow.cpp new file mode 100644 index 0000000..ea5d066 --- /dev/null +++ b/subprojects/MinCostFlow/examples/kattis-maxflow.cpp @@ -0,0 +1,72 @@ +// https://open.kattis.com/problems/maxflow + +#include +#include +#include +#include + +using value_type = int; +// typedef ln::maxflow_augmenting_path maxflow_t; // 1.04s +// typedef ln::maxflow_augmenting_path maxflow_t; // 0.04s +// typedef ln::maxflow_scaling maxflow_t; // 0.42s +typedef ln::maxflow_scaling maxflow_t; // 0.07s +// typedef ln::maxflow_preflow maxflow_t; // 0.53s + +int main() +{ + int N,M,S,T; + std::cin >> N >> M >> S >> T; + + + + ln::digraph Graph; + std::vector capacity; + + for(int e=0;e>a>>b>>c; + auto [arc,arc2] = Graph.add_arc(a,b,e); + + if(capacity.size() 0 ? 1 : 0); + } + + + std::cout << N << ' ' << max_flow << ' ' << M_count << '\n'; + + for(int i=0;i +#include + +typedef long long value_type; + +// typedef ln::maxflow_augmenting_path maxflow_t; +typedef ln::maxflow_augmenting_path maxflow_t; +// typedef ln::maxflow_scaling maxflow_t; +// typedef ln::maxflow_scaling maxflow_t; +// typedef ln::maxflow_preflow maxflow_t; + +//typedef ln::mincostflow_EdmondsKarp< +// value_type,ln::shortestPath_FIFO> mincostflow_t; // 2.18s +// typedef ln::mincostflow_EdmondsKarp< +// value_type,ln::shortestPath_BellmanFord> mincostflow_t; // TLE +// typedef ln::mincostflow_EdmondsKarp mincostflow_t; // expected failure + +// typedef ln::mincostflow_PrimalDual< +// ln::shortestPath_FIFO,maxflow_t> mincostflow_t; // 2.54s, 3.02s, 2.75s, 3.18s, TLE +// typedef ln::mincostflow_PrimalDual< +// ln::shortestPath_BellmanFord,maxflow_t> mincostflow_t; // TLE +//typedef ln::mincostflow_PrimalDual< +// ln::shortestPath_Dijkstra,maxflow_t> mincostflow_t; +// 1.34s, 1.76s, 1.52s, 2.42s, TLE + +//typedef ln::mincostflow_capacityScaling< +// ln::shortestPath_Dijkstra,maxflow_t> mincostflow_t; +// TLE +typedef ln::mincostflow_costScaling< + maxflow_t> mincostflow_t; //0.37s, 0.06s, 0.10s, 0.07s, 0.15s + +int main() +{ + int N,M,S,T; + std::cin >> N >> M >> S >> T; + + ln::digraph G; + std::vector capacity; + std::vector weight; + + G.add_node(S); + G.add_node(T); + + for(int e=0;e>a>>b>>c>>w; + auto [arc,arc2] = G.add_arc(a,b,e); + + capacity.resize(G.max_num_arcs()); + weight.resize(G.max_num_arcs()); + + capacity[arc] = c; + capacity[arc2]=0; + + weight[arc] = w; + weight[arc2]=-w; + } + + mincostflow_t f; + + auto Flow = f.solve(G,G.get_node(S),G.get_node(T),weight,capacity); + long long Cost = 0; + + for(int e=0;e(weight[arc])*f.flow_at(G,arc,capacity); + } + + std::cout << Flow << " " << Cost << '\n'; + + return 0; +} + diff --git a/subprojects/MinCostFlow/examples/kattis-shortestpath1.cpp b/subprojects/MinCostFlow/examples/kattis-shortestpath1.cpp new file mode 100644 index 0000000..7da53aa --- /dev/null +++ b/subprojects/MinCostFlow/examples/kattis-shortestpath1.cpp @@ -0,0 +1,64 @@ +// https://open.kattis.com/problems/shortestpath1 + +#include "mincostflow/graph.hpp" +#include "mincostflow/shortestpath.hpp" +#include + +using length_type = int; +typedef ln::shortestPath_Dijkstra pathsolver_t; // 0.28s +//typedef ln::shortestPath_FIFO pathsolver_t; // 1.14s +//typedef ln::shortestPath_BellmanFord pathsolver_t; // >3.0s + +int main() +{ + + while(1) + { + int N_vertex,N_edges,Q,S; + std::cin>>N_vertex>>N_edges>>Q>>S; + if(N_vertex==0)break; + + ln::digraph Graph; + std::vector weights(N_edges); + pathsolver_t solver; + + + for(int e=0;e> a>>b>>w; + auto [arc,arc2] = Graph.add_arc(a,b,e); + + if(weights.size()>v; + auto pos = Graph.get_node(v); + if(Graph.is_valid(pos) && solver.distance.at(pos) + +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace ln +{ + + class digraph_types + { + public: + + // a lot more efficient if we can internally identify arcs and nodes by their unique + // position in the buffer array + typedef std::size_t pos_type; + static constexpr pos_type NONE = std::numeric_limits::max(); + + struct node_pos_t + { + pos_type x{NONE}; + bool operator<(const node_pos_t& that)const + { + return x < that.x; + } + bool operator==(const node_pos_t& that)const + { + return x==that.x; + } + + operator pos_type() const + { + return x; + } + node_pos_t& operator++() + { + ++x; + return *this; + } + node_pos_t operator++(int) + { + ++x; + return node_pos_t{x-1}; + } + }; + struct arc_pos_t + { + pos_type x{NONE}; + bool operator<(const arc_pos_t& that)const + { + return x < that.x; + } + bool operator==(const arc_pos_t& that)const + { + return x==that.x; + } + operator pos_type() const + { + return x; + } + arc_pos_t& operator++() + { + ++x; + return *this; + } + arc_pos_t operator++(int) + { + ++x; + return arc_pos_t{x-1}; + } + }; + + struct arc_data_t + { + node_pos_t a{NONE},b{NONE}; + arc_pos_t dual{NONE}; + }; + + struct node_data_t + { + std::vector out_arcs,in_arcs; + + void rm_arc(arc_pos_t arc) + { + auto rm_arc_vec=[arc](std::vector& V) + { + if (auto ptr = std::find(V.begin(),V.end(),arc); + ptr!=V.end()) + { + *ptr = V.back(); + V.pop_back(); + } + }; + rm_arc_vec(in_arcs); + rm_arc_vec(out_arcs); + } + void add_in_arc(arc_pos_t arc) + { + in_arcs.push_back(arc); + } + void add_out_arc(arc_pos_t arc) + { + out_arcs.push_back(arc); + } + }; + }; + + + // TODO: template on custom allocator + template + class digraph : public digraph_types + /* + Represents: a directed graph with dual arcs to simulate the residual network. + This data structure represents only the topological information. + Nodes and arcs have fixed positions. + */ + /* + ideally I would like to: + + digraph g; + + g.add_arc(node_a,node_b,arc_ab); // adds also the dual + + g::node_id_t n_1 = g.source_node(arc_x); + g::node_id_t n_2 = g.dest_node(arc_x); + + g::node_id_t n = g.add_node(node_factory); // generates a new node + + + + // an optimized interface allows to pass additional structure as arrays + g.max_nodes(); // size of the node array + g.max_arcs() ; // size of arc array + + for() + */ + { + vectorized_map my_arcs; + vectorized_map my_nodes; + + std::unordered_map arcs_htable; + std::vector arcs_ids; + std::vector arcs_ids_flag; + + std::unordered_map nodes_htable; + std::vector nodes_ids; + std::vector nodes_ids_flag; + + public: + + const auto& arcs()const + { + return my_arcs; + } + const auto& nodes()const + { + return my_nodes; + } + + bool is_valid(arc_pos_t arc)const + { + return my_arcs.is_valid(arc); + } + bool is_valid(node_pos_t node)const + { + return my_nodes.is_valid(node); + } + bool has_id(node_pos_t node)const + { + return nodes_ids_flag.at(node); + } + bool has_id(arc_pos_t arc)const + { + return arcs_ids_flag.at(arc); + } + + auto arc_ends(arc_pos_t arc)const + { + return std::pair{ + my_arcs.at(arc).a, + my_arcs.at(arc).b + }; + } + arc_pos_t arc_dual(arc_pos_t arc)const + { + return my_arcs.at(arc).dual; + } + + auto arc_ends_nocheck(arc_pos_t arc)const noexcept + { + return std::pair{ + my_arcs[arc].a, + my_arcs[arc].b + }; + } + arc_pos_t arc_dual_nocheck(arc_pos_t arc)const noexcept + { + return my_arcs[arc].dual; + } + + void erase(arc_pos_t arc) + { + if(! is_valid(arc)) + return; + + auto [a,b] = arc_ends(arc); + + my_nodes.at(a).rm_arc(arc); + my_nodes.at(b).rm_arc(arc); + + if(has_id(arc)) + { + auto id = arcs_ids.at(arc); + arcs_htable.erase(id); + } + + my_arcs.erase(arc); + arcs_ids.resize(my_arcs.capacity()); + arcs_ids_flag.resize(my_arcs.capacity()); + } + void erase(node_pos_t node) + { + if(! is_valid(node)) + return; + std::vector ls_arcs; + std::copy(my_nodes.at(node).in_arcs.begin(), + my_nodes.at(node).in_arcs.end(), + std::back_inserter(ls_arcs)); + std::copy(my_nodes.at(node).out_arcs.begin(), + my_nodes.at(node).out_arcs.end(), + std::back_inserter(ls_arcs)); + + // first remove all incoming and outgoin arcs + for(auto arc: ls_arcs) + erase(arc); + + if(has_id(node)) + { + auto id = nodes_ids.at(node); + nodes_htable.erase(id); + } + my_nodes.erase(node); + nodes_ids.resize(my_nodes.capacity()); + nodes_ids_flag.resize(my_nodes.capacity()); + } + node_pos_t new_node() + { + node_pos_t node = my_nodes.insert(node_data_t{}); + + nodes_ids.resize(my_nodes.capacity()); + nodes_ids_flag.resize(my_nodes.capacity()); + + nodes_ids_flag.at(node)=false; + + return node; + } + + const auto& out_arcs(node_pos_t node)const + { + if(!is_valid(node)) + throw std::runtime_error( + "digraph::out_arcs invalid node"); + return my_nodes.at(node).out_arcs; + } + const auto& in_arcs(node_pos_t node)const + { + if(!is_valid(node)) + throw std::runtime_error( + "digraph::in_arcs invalid node"); + return my_nodes.at(node).in_arcs; + } + + arc_pos_t new_arc(node_pos_t a, node_pos_t b) + { + if(!is_valid(a) || !is_valid(b)) + throw std::runtime_error("digraph::new_arc add a new arc with invalid end nodes"); + + arc_pos_t arc = my_arcs.insert(arc_data_t{a,b,NONE}); + arcs_ids.resize(my_arcs.capacity()); + arcs_ids_flag.resize(my_arcs.capacity()); + + arcs_ids_flag.at(arc)=false; + + my_nodes.at(a).add_out_arc(arc); + my_nodes.at(b).add_in_arc(arc); + return arc; + } + void set_dual(arc_pos_t arc1, arc_pos_t arc2) + { + if(!is_valid(arc1) || !is_valid(arc2)) + throw std::runtime_error("digraph::set_dual invalid arcs"); + + my_arcs.at(arc1).dual = arc2; + my_arcs.at(arc2).dual = arc1; + } + + auto max_num_arcs()const + { + return my_arcs.capacity(); + } + auto num_arcs()const + { + return my_arcs.size(); + } + auto max_num_nodes()const + { + return my_nodes.capacity(); + } + auto num_nodes()const + { + return my_nodes.size(); + } + + // translation + node_id_t get_node_id(node_pos_t node)const + { + if(!is_valid(node)) + throw std::runtime_error("digraph::get_node_id invalid node"); + if(!has_id(node)) + throw std::runtime_error("digraph::get_node_id node without id"); + return nodes_ids.at(node); + } + arc_id_t get_arc_id(arc_pos_t arc)const + { + if(!is_valid(arc)) + throw std::runtime_error("digraph::get_arc_id invalid arc"); + if(!has_id(arc)) + throw std::runtime_error("digraph::get_arc_id arc without id"); + return arcs_ids.at(arc); + } + node_pos_t get_node(node_id_t id)const + { + node_pos_t node{NONE}; + if(auto ptr = nodes_htable.find(id); + ptr != nodes_htable.end()) + { + node = ptr->second; + } + return node; + } + arc_pos_t get_arc(arc_id_t id)const + { + arc_pos_t arc{NONE}; + if(auto ptr = arcs_htable.find(id); + ptr != arcs_htable.end()) + { + arc = ptr->second; + } + return arc; + } + node_pos_t add_node(node_id_t id) + { + auto node = get_node(id); + if(!is_valid(node)) + { + node = new_node(); + nodes_ids.at(node) = id; + nodes_ids_flag.at(node) = true; + nodes_htable[id] = node; + } + return node; + } + std::pair add_arc(node_id_t a, node_id_t b, arc_id_t id) + { + auto n_a = add_node(a); + auto n_b = add_node(b); + + if(auto arc = get_arc(id); is_valid(arc)) + throw std::runtime_error("digraph::add_arc arc id already exists"); + + auto arc1 = new_arc(n_a,n_b); + auto arc2 = new_arc(n_b,n_a); + set_dual(arc1,arc2); + + arcs_ids.at(arc1) = id; + arcs_ids_flag.at(arc1) = true; + arcs_htable[id] = arc1; + + return {arc1,arc2}; + } + void remove_node(node_id_t id) + { + auto node = get_node(id); + if(!is_valid(node)) + return; + nodes_htable.erase(id); + erase(node); + } + void remove_arc(arc_id_t id) + { + auto arc = get_arc(id); + if(!is_valid(arc)) + return; + arcs_htable.erase(id); + + auto arc2 = arc_dual(arc); + erase(arc); + erase(arc2); + } + + digraph() + {} + }; +} diff --git a/subprojects/MinCostFlow/include/mincostflow/maxflow.hpp b/subprojects/MinCostFlow/include/mincostflow/maxflow.hpp new file mode 100644 index 0000000..0f9da7c --- /dev/null +++ b/subprojects/MinCostFlow/include/mincostflow/maxflow.hpp @@ -0,0 +1,289 @@ +#pragma once + +#include + +namespace ln +{ + + template + class maxflow_base : public digraph_types + { + public: + using value_type = T; + static constexpr value_type INFINITY = std::numeric_limits::max(); + + template + value_type flow_at( + const graph_t& g, + const arc_pos_t e, + const std::vector& capacity) + { + auto e2 = g.arc_dual(e); + return capacity.at(e2.x); + } + }; + + template + class maxflow_augmenting_path : public maxflow_base + { + public: + using base_type = maxflow_base; + using value_type = typename base_type::value_type; + using node_pos_t = typename base_type::node_pos_t; + using arc_pos_t = typename base_type::arc_pos_t; + using base_type::flow_at; + using base_type::INFINITY; + + + + template + value_type solve( + const graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + std::vector& capacity, + condition_t valid_arc) + { + value_type sent=0; + path_solver_type path_solver; + + while(1) + { + bool found = path_solver.solve( + g, + Source,Dest, + [valid_arc,&capacity](arc_pos_t e) + { + return capacity.at(e)>0 && valid_arc(e); + }); + + + if(!found) + break; + + auto path = path_solver.get_path(g,Dest); + + value_type k = INFINITY; + for(auto e : path) + { + k = std::min(k,capacity.at(e)); + } + + for(auto e: path) + { + capacity.at(e) -= k; + capacity.at(g.arc_dual(e)) += k; + } + + sent += k; + } + return sent; + } + + maxflow_augmenting_path() + {} + }; + + template + class maxflow_scaling : public maxflow_base + { + public: + using base_type = maxflow_base; + using value_type = typename base_type::value_type; + using node_pos_t = typename base_type::node_pos_t; + using arc_pos_t = typename base_type::arc_pos_t; + using base_type::flow_at; + using base_type::INFINITY; + + template + value_type solve( + const graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + std::vector& residual_cap, + condition_t valid_arc) + // augmenting path + { + value_type sent=0; + path_solver_type search_algo; + + value_type cap_flow = 1; + for(auto e : g.out_arcs(Source)) + cap_flow = std::max(cap_flow,residual_cap.at(e)); + + cap_flow = lower_bound_power2(cap_flow); + + // int cycle=0; + for(;cap_flow>0;) + { + // cycle++; + // std::cerr << "augmenting path cycle: " << cycle << '\n'; + // std::cerr << "flow sent: " << sent << '\n'; + // std::cerr << "cap flow: " << cap_flow << '\n'; + + bool found = search_algo.solve( + g, + Source,Dest, + // edge is valid if + [this,valid_arc,cap_flow,&residual_cap](arc_pos_t e) + { + return residual_cap.at(e)>=cap_flow && valid_arc(e); + }); + + if(! found) + { + cap_flow/=2; + // std::cerr << "path not found!\n"; + search_algo.reset(); + continue; + } + + auto path = search_algo.get_path(g,Dest); + + // std::cerr << "path found!\n"; + + for(auto e: path) + { + residual_cap[e] -= cap_flow; + residual_cap[g.arc_dual(e)] += cap_flow; + } + + sent += cap_flow; + } + return sent; + } + + maxflow_scaling() + {} + }; + + + template + class maxflow_preflow : public maxflow_base, public distance_structure + { + public: + using base_type = maxflow_base; + using value_type = typename base_type::value_type; + using node_pos_t = typename base_type::node_pos_t; + using arc_pos_t = typename base_type::arc_pos_t; + using base_type::flow_at; + using base_type::INFINITY; + + static constexpr auto flow_INFINITY = base_type::INFINITY; + static constexpr auto dist_INFINITY = distance_structure::INFINITY; + + std::vector excess; + + template + void initialize_distance( + const graph_t& g, + const node_pos_t Dest, + condition_t valid_arc) + { + distance_structure::init(g); + + distance.at(Dest)=0; + + std::queue q; + q.push(Dest); + + while(!q.empty()) + { + auto n = q.front(); + q.pop(); + + for(auto e: g.in_arcs(n)) + if( valid_arc(e) ) + { + // assert b==n + auto [a,b] = g.arc_ends(e); + int dnew = distance[b] + 1; + + if(distance[a]==dist_INFINITY) + { + distance[a] = dnew; + q.push(a); + } + } + } + } + public: + + template + value_type solve( + const graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + std::vector& residual_cap, + condition_t valid_arc) + { + excess.resize(g.max_num_nodes()); + std::fill(excess.begin(),excess.end(),0); + + initialize_distance(g,Dest,valid_arc); + std::queue q; + + auto push = [&](arc_pos_t e) + { + auto [a,b] = g.arc_ends(e); + const auto delta = std::min(excess[a],residual_cap.at(e)); + residual_cap.at(e) -= delta; + residual_cap.at(g.arc_dual(e)) += delta; + + assert(delta>=0); + + excess.at(a) -= delta; + excess.at(b) += delta; + + if(delta>0 && excess.at(b)==delta) + q.push(b); + }; + + auto relabel = [&](node_pos_t v) + { + int hmin = dist_INFINITY; + for(auto e : g.out_arcs(v)) + if(valid_arc(e) && residual_cap.at(e)>0) + hmin = std::min(hmin,distance.at(g.arc_ends(e).second)); + if(hmin0) + { + auto b = g.arc_ends(e).second; + if(distance[a]== distance[b]+1) + push(e); + } + + if(excess.at(a)==0) + break; + + relabel(a); + } + }; + + excess.at(Source) = flow_INFINITY; + distance.at(Source) = g.num_nodes(); + + for(auto e : g.out_arcs(Source)) + if(valid_arc(e)) + push(e); + + while(!q.empty()) + { + auto node = q.front(); + q.pop(); + + if(node!=Dest && node!=Source) + discharge(node); + } + return excess.at(Dest); + } + maxflow_preflow() + {} + + }; +} diff --git a/subprojects/MinCostFlow/include/mincostflow/meson.build b/subprojects/MinCostFlow/include/mincostflow/meson.build new file mode 100644 index 0000000..42a02b3 --- /dev/null +++ b/subprojects/MinCostFlow/include/mincostflow/meson.build @@ -0,0 +1 @@ +my_headers = files(['graph.hpp','shortestpath.hpp','maxflow.hpp','mincostflow.hpp']) diff --git a/subprojects/MinCostFlow/include/mincostflow/mincostflow.hpp b/subprojects/MinCostFlow/include/mincostflow/mincostflow.hpp new file mode 100644 index 0000000..afc1530 --- /dev/null +++ b/subprojects/MinCostFlow/include/mincostflow/mincostflow.hpp @@ -0,0 +1,479 @@ +#pragma once + +#include +#include + +namespace ln +{ + template + class mincostflow_EdmondsKarp : public maxflow_base + { + public: + using base_type = maxflow_base; + using value_type = typename base_type::value_type; + using node_pos_t = typename base_type::node_pos_t; + using arc_pos_t = typename base_type::arc_pos_t; + using base_type::flow_at; + using base_type::INFINITY; + + template + value_type solve( + const graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + const std::vector& weight, + std::vector& residual_cap + ) + // augmenting path + { + value_type sent =0 ; + path_optimizer_type path_opt; + + while(true) + { + path_opt.solve( + g, + Source, + weight, + // edge is valid if + [&residual_cap](arc_pos_t e){ + return residual_cap.at(e)>0; + }); + + if(! path_opt.is_reacheable(Dest)) + break; + + auto path = path_opt.get_path(g,Dest); + + value_type k = INFINITY; + for(auto e : path) + { + k = std::min(k,residual_cap.at(e)); + } + + for(auto e: path) + { + residual_cap[e] -= k; + residual_cap[g.arc_dual(e)] += k; + } + + sent += k; + } + return sent; + } + + mincostflow_EdmondsKarp() + {} + }; + + + template + class mincostflow_PrimalDual : public maxflow_type + { + public: + using base_type = maxflow_type; + using value_type = typename base_type::value_type; + using node_pos_t = typename base_type::node_pos_t; + using arc_pos_t = typename base_type::arc_pos_t; + using base_type::flow_at; + using base_type::INFINITY; + + + template + value_type solve( + const graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + const std::vector& weight, + std::vector& residual_cap + ) + { + std::vector reduced_weight = weight; + + value_type sent =0 ; + path_optimizer_type path_opt; + + while(true) + { + path_opt.solve( + g, + Source, + reduced_weight, + // edge is valid if + [&residual_cap](arc_pos_t e) -> bool + { + return residual_cap.at(e)>0; + }); + + if(! path_opt.is_reacheable(Dest)) + break; + + const auto& distance{path_opt.distance}; + + for(auto e : g.arcs()) + { + + auto [a,b] = g.arc_ends(e); + if(distance[a]bool + { + return reduced_weight[e]==0; + }); + + sent += F; + } + return sent; + } + + mincostflow_PrimalDual() + {} + }; + + template + class mincostflow_capacityScaling : public maxflow_type + { + public: + using base_type = maxflow_type; + using value_type = typename base_type::value_type; + using node_pos_t = typename base_type::node_pos_t; + using arc_pos_t = typename base_type::arc_pos_t; + using base_type::flow_at; + using base_type::INFINITY; + + template + value_type solve( + graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + const std::vector& weight, + std::vector& residual_cap) + { + + std::vector reduced_weight = weight; + value_type maxflow{0}; + + // find the max-flow-anycost + maxflow = maxflow_type::solve( + g,Source,Dest, + residual_cap, + [](arc_pos_t)->bool{return true;}); + + value_type cap_flow = lower_bound_power2(maxflow); + + std::vector excess(g.max_num_nodes(),0); + + std::vector weight_ex = weight; + + auto update_reduced_costs = + [&](const std::vector& potential) + { + for(auto e : g.arcs()) + { + auto [src,dst] = g.arc_ends(e); + auto p_src = potential.at(src), p_dst = potential.at(dst); + + p_src = p_src == INFINITY ? 0 : p_src; + p_dst = p_dst == INFINITY ? 0 : p_dst; + + weight_ex.at(e) += p_src - p_dst; + } + }; + + auto push_flow = + [&](arc_pos_t e,value_type delta) + { + // std::cerr << " push flow at " << e << " delta = " << delta << "\n"; + auto [src,dst] = g.arc_ends(e); + + residual_cap[e]-=delta; + residual_cap[g.arc_dual(e)]+=delta; + + excess.at(src) -= delta; + excess.at(dst) += delta; + + // std::cerr << "push " << delta << " over " << e << '\n'; + }; + + // auto report = + // [&]() + // { + // std::cerr << "residual cap + mod. costs\n"; + // for(auto e : g.arcs()) + // { + // std::cerr << " " << e << " -> " << residual_cap[e] << " " << weight_ex[e] << "\n"; + // } + // std::cerr << "potential + excess\n"; + // for(auto v : g.nodes()) + // { + // std::cerr << " " << v << " -> " << excess[v] << "\n"; + // } + // }; + // + // std::cerr << " maxflow = " << maxflow << "\n"; + + // int cycle=0; + for(;cap_flow>0;cap_flow/=2) + { + // cycle++; + // std::cerr << "cycle " << cycle << " cap_flow = " << cap_flow << '\n'; + // report(); + + // saturate edges with negative cost + for(auto e : g.arcs()) + while(residual_cap.at(e)>=cap_flow && weight_ex.at(e)<0) + { + push_flow(e,cap_flow); + } + + path_optimizer_type path_opt; + + // build S and T + std::set Sset,Tset; + for(auto v : g.nodes()) + { + if(excess.at(v)>=cap_flow) + Sset.insert(v); + if(excess.at(v)<=-cap_flow) + Tset.insert(v); + } + + const auto multi_source_node = g.new_node(); + excess.resize(g.max_num_nodes()); + excess.at(multi_source_node) = 0; + const Scope_guard rm_node = [&](){ g.erase(multi_source_node);}; + + + + for(auto v : Sset) + { + auto arc1 = g.new_arc(multi_source_node,v); + auto arc2 = g.new_arc(v,multi_source_node); + + g.set_dual(arc1,arc2); + + weight_ex.resize(g.max_num_arcs()); + residual_cap.resize(g.max_num_arcs()); + + weight_ex.at(arc1) = 0; + residual_cap.at(arc1) = excess.at(v); + + weight_ex.at(arc2) = 0; + residual_cap.at(arc2) = 0; + + excess.at(multi_source_node) += excess.at(v); + excess.at(v) = 0; + } + + const Scope_guard restore_excess = [&]() + { + for(auto e : g.out_arcs(multi_source_node)) + { + auto [src,dst] = g.arc_ends(e); + excess.at(dst) = residual_cap.at(e); + } + }; + + while(!Sset.empty() && !Tset.empty()) + { + path_opt.solve( + g, multi_source_node, + weight_ex, + [cap_flow,&residual_cap](arc_pos_t e)->bool + { + return residual_cap.at(e)>=cap_flow; + } + ); + + const auto& distance{path_opt.distance}; + + auto it = std::find_if(Tset.begin(),Tset.end(), + [&](node_pos_t v)->bool { + return distance.at(v) " << distance[v]<<"\n"; + // } + + update_reduced_costs(distance); + + auto path = path_opt.get_path(g,dst); + for(auto e: path) + { + // auto [src,dst] = g.arc_ends(e); + push_flow(e,cap_flow); + } + + if(excess.at(dst)>-cap_flow) + Tset.erase(dst); + } + } + + return maxflow; + } + + public: + mincostflow_capacityScaling() + {} + }; + + template + class mincostflow_costScaling : public maxflow_type + { + public: + using base_type = maxflow_type; + using value_type = typename base_type::value_type; + using node_pos_t = typename base_type::node_pos_t; + using arc_pos_t = typename base_type::arc_pos_t; + using base_type::flow_at; + using base_type::INFINITY; + + template + value_type solve( + graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + const std::vector& weight, + std::vector& residual_cap) + { + value_type maxflow{0}; + + // find the max-flow-anycost + maxflow = maxflow_type::solve( + g,Source,Dest, + residual_cap, + [](arc_pos_t)->bool{return true;}); + + // return maxflow; + + std::vector reduced_weight = weight; + std::vector potential(g.max_num_nodes(),0); + std::vector excess(g.max_num_nodes(),0); + + auto relabel = + [&](const node_pos_t x,value_type eps) + { + // std::cerr << " relabel " << x << " delta = " << eps << "\n"; + + potential[x] -= eps; + for(auto e : g.out_arcs(x)) + { + reduced_weight[e] -= eps; + } + for(auto e : g.in_arcs(x)) + { + reduced_weight[e] += eps; + } + }; + + auto push_flow = + [&](arc_pos_t e,value_type delta) + { + auto [src,dst] = g.arc_ends(e); + // std::cerr << " push flow at " << e + // << " (" << src << "," << dst << ")" + // << " delta = " << delta << "\n"; + + residual_cap[e]-=delta; + residual_cap[g.arc_dual(e)]+=delta; + + excess[src] -= delta; + excess[dst] += delta; + + // std::cerr << "push " << delta << " over " << e << '\n'; + }; + + value_type maxC = 0; + const int N = g.num_nodes(); + for(auto e : g.arcs()) + { + reduced_weight[e] *= N; + maxC = std::max(maxC,reduced_weight[e]); + } + maxC = lower_bound_power2(maxC); + + + // int cycle=0; + for(;maxC>0;maxC/=2){ + // cycle++; + // std::cerr << "cycle " << cycle << " eps = " << maxC << '\n'; + + // improve + for(auto e : g.arcs()) + { + if(reduced_weight[e]<0 && residual_cap[e]>0) + { + push_flow(e,residual_cap[e]); + } + // this also does the trick of making the flow = 0 for arcs with + // reduced_weight>0, + } + std::set active; + for(auto n : g.nodes()) + if(excess[n]>0) + { + active.insert(n); + // std::cerr << n << " is added to active\n"; + } + // int ct=0; + while(!active.empty()) + { + // ct++; + // if(ct>20)exit(1); + + auto t = *active.begin(); + + // std::cerr << t << " is active\n"; + + bool pushed = false; + + for(auto e : g.out_arcs(t)) + { + auto rw = reduced_weight[e]; + auto rc = residual_cap[e]; + if(rw<0 && rw>=-maxC && rc>0) + { + pushed = true; + auto [a,b] = g.arc_ends(e); + auto d = std::min(excess[a],rc); + + push_flow(e,d); + + if(excess[a]<=0) + active.erase(a); + if(excess[b]>0) + active.insert(b); + + break; + } + } + + if(!pushed) + relabel(t,maxC); + } + } + // std::cerr << "done\n"; + return maxflow; + } + + public: + mincostflow_costScaling() + {} + }; + +} diff --git a/subprojects/MinCostFlow/include/mincostflow/scope_guard.hpp b/subprojects/MinCostFlow/include/mincostflow/scope_guard.hpp new file mode 100644 index 0000000..3624bd8 --- /dev/null +++ b/subprojects/MinCostFlow/include/mincostflow/scope_guard.hpp @@ -0,0 +1,45 @@ +#pragma once + +/* + credits to + https://stackoverflow.com/questions/31365013/what-is-scopeguard-in-c +*/ + +#include // std::function +#include // std::move + +namespace ln +{ + class Non_copyable + { + private: + auto operator=( Non_copyable const& ) -> Non_copyable& = delete; + Non_copyable( Non_copyable const& ) = delete; + public: + auto operator=( Non_copyable&& ) -> Non_copyable& = default; + Non_copyable() = default; + Non_copyable( Non_copyable&& ) = default; + }; + + class Scope_guard + : public Non_copyable + { + private: + std::function cleanup_; + + public: + friend + void dismiss( Scope_guard& g ) { g.cleanup_ = []{}; } + + ~Scope_guard() { cleanup_(); } + + template< class Func > + Scope_guard( Func const& cleanup ) + : cleanup_( cleanup ) + {} + + Scope_guard( Scope_guard&& other ) + : cleanup_( std::move( other.cleanup_ ) ) + { dismiss( other ); } + }; +} diff --git a/subprojects/MinCostFlow/include/mincostflow/shortestpath.hpp b/subprojects/MinCostFlow/include/mincostflow/shortestpath.hpp new file mode 100644 index 0000000..127726b --- /dev/null +++ b/subprojects/MinCostFlow/include/mincostflow/shortestpath.hpp @@ -0,0 +1,518 @@ +#pragma once + +#include + +#include +#include +#include +#include +#include +#include + + +namespace ln +{ + inline long long int lower_bound_power2(long long int n) + { + if(n<=2) return n; + while(n != (n & -n)) + n -= (n & -n); + return n; + } + + + class parent_structure : public digraph_types + { + public: + std::vector parent; + + bool has_parent(node_pos_t x)const + { + return parent.at(x)!=arc_pos_t{NONE}; + } + bool is_reacheable(node_pos_t x)const + { + return has_parent(x); + } + + template + void init(const graph_t& g) + { + parent.resize(g.max_num_nodes()); + std::fill(parent.begin(),parent.end(),arc_pos_t{NONE}); + } + + + template + auto get_path(const graph_t& g,node_pos_t last)const + { + std::vector path; + while(1) + { + auto e = parent.at(last); + if(!g.is_valid(e)) + break; + + path.push_back(e); + auto [a,b] = g.arc_ends(e); + last = a; + + } + return path; + } + }; + + template + class distance_structure + { + public: + using value_type = T; + static constexpr value_type INFINITY = std::numeric_limits::max(); + std::vector distance; + + template + void init(const graph_t& g) + { + distance.resize(g.max_num_nodes()); + std::fill(distance.begin(),distance.end(),INFINITY); + } + }; + + + class pathSearch_BFS : public parent_structure, public distance_structure + /* + Represents: weigthless path using BFS + Invariant: + + User interface: + Complexity: |E|+|V| + */ + { + public: + using value_type = int; + using parent_structure::parent; + using parent_structure::init; + using distance_structure::distance; + using distance_structure::init; + using distance_structure::INFINITY; + + pathSearch_BFS() + { + } + + void reset() + {} + + template + bool solve ( + const graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + condition_t valid_arc) + // each call resets the state of the tree and distances + // O(|E|+|V|) + { + bool found = false; + + if(!g.is_valid(Source)) + throw std::runtime_error( + "pathSearch_BFS::solve source node is not valid"); + + if(!g.is_valid(Dest)) + throw std::runtime_error( + "pathSearch_BFS::solve destination node is not valid"); + + parent_structure::init(g); + distance_structure::init(g); + + distance.at(Source) = 0; + + std::queue q; + q.push(Source); + + while(!q.empty()) + { + auto node = q.front(); + q.pop(); + + if(node==Dest) + { + found = true; + break; + } + for(auto e: g.out_arcs(node)) + if( valid_arc(e) ) + { + auto [a,b] = g.arc_ends(e); + + if(distance.at(b)==INFINITY) + { + distance.at(b) = distance.at(a)+1; + parent.at(b) = e; + q.push(b); + } + } + } + return found; + } + }; + + class pathSearch_labeling : public parent_structure, + public distance_structure + /* + Represents: shortest path with labeling + Invariant: + + User interface: + Complexity: + */ + { + public: + using value_type = distance_structure::value_type; + using parent_structure::parent; + using parent_structure::init; + using distance_structure::distance; + using distance_structure::init; + using distance_structure::INFINITY; + + node_pos_t last_source{NONE},last_dest{NONE}; + std::vector dist_freq; + + template + void initialize ( + const graph_t &g, + condition_t valid_arc) + { + parent_structure::init(g); + distance_structure::init(g); + + dist_freq.resize(g.num_nodes()+1); + std::fill(dist_freq.begin(),dist_freq.end(),0); + + std::queue q; + distance.at(last_dest)=0; + + // TODO: write a general purpose BFS label solver + q.push(last_dest); + + while(!q.empty()) + { + auto n = q.front(); + q.pop(); + + for(auto e: g.in_arcs(n)) + if( valid_arc(e) ) + { + auto [a,b] = g.arc_ends(e); + value_type dnew = distance[b] + 1; + + if(distance[a]==INFINITY) + { + distance[a] = dnew; + dist_freq.at(dnew)++; + q.push(a); + } + } + } + } + + + public: + + pathSearch_labeling() + {} + + + void reset() + { + last_source = last_dest = node_pos_t{NONE}; + } + template + bool solve ( + const graph_t& g, + const node_pos_t Source, const node_pos_t Dest, + condition_t valid_arc) + { + if(last_source!=Source || last_dest!=Dest) + { + last_source = Source; + last_dest = Dest; + initialize(g,valid_arc); + } + + parent_structure::init(g); + + + for(auto current = Source; + distance.at(Source) + class shortestPath_FIFO : public parent_structure, public distance_structure + /* + Represents: shortest path label-correcting FIFO + Invariant: + + User interface: + Complexity: pseudo-polynomial + */ + { + public: + using value_type = T; + using parent_structure::parent; + using parent_structure::init; + using distance_structure::distance; + using distance_structure::init; + using distance_structure::INFINITY; + + shortestPath_FIFO() + {} + + template + void solve( + const graph_t& g, + const node_pos_t Source, + const std::vector& weight, + condition_t valid_arc) + // shortest path FIFO + // each call resets the state of the tree and distances + // O( pseudo-polynomial ) + { + parent_structure::init(g); + distance_structure::init(g); + + if(!g.is_valid(Source)) + throw std::runtime_error( + "shortestPath_FIFO::solve source node is not valid"); + + if(weight.size() q; + q.push(Source); + distance.at(Source)=0; + + while(!q.empty()) + { + auto node = q.front(); + q.pop(); + + for(auto e: g.out_arcs(node)) + if( valid_arc(e) ) + { + auto [a,b] = g.arc_ends(e); + const value_type dnew = distance.at(a)+weight.at(e); + + if(distance.at(b)>dnew) + { + distance.at(b) = dnew; + parent.at(b) = e; + q.push(b); + } + } + } + } + }; + + template + class shortestPath_BellmanFord : public parent_structure, public distance_structure + /* + Represents: shortest path using Bellman-Ford + Invariant: + + User interface: + Complexity: |V| |E| + */ + { + public: + using value_type = T; + using parent_structure::parent; + using parent_structure::init; + using distance_structure::distance; + using distance_structure::init; + using distance_structure::INFINITY; + + shortestPath_BellmanFord() + {} + + template + void solve ( + const graph_t& g, + const node_pos_t Source, + const std::vector& weight, + condition_t valid_arc) + // shortest path Bellman-Ford + // each call resets the state of the tree and distances + // O(|V||E|) + { + parent_structure::init(g); + distance_structure::init(g); + + if(!g.is_valid(Source)) + throw std::runtime_error( + "shortestPath_BellmanFord::solve source node is not valid"); + + if(weight.size()dnew) + { + distance[b]=dnew; + parent.at(b) = e; + updates = true; + } + } + } + if(! updates) + break; + } + // TODO: check for negative cycles + } + }; + + template + class shortestPath_Dijkstra : public parent_structure, public distance_structure + /* + Represents: shortest path with weights using Dijkstra + Invariant: + + User interface: + Complexity: |E| + |V| log |V| + */ + { + public: + using value_type = T; + using parent_structure::parent; + using parent_structure::init; + using distance_structure::distance; + using distance_structure::init; + using distance_structure::INFINITY; + + shortestPath_Dijkstra() + {} + + template + void solve ( + const graph_t& g, + const node_pos_t Source, + const std::vector& weight, + condition_t valid_arc) + // Dijkstra algorithm + // precondition: doesnt work with negative weights! + // O( |E|+|V| log |V| ) + { + parent_structure::init(g); + distance_structure::init(g); + + if(!g.is_valid(Source)) + throw std::runtime_error( + "shortestPath_Dijkstra::solve source node is not valid"); + + if(weight.size() visited(g.max_num_nodes(),false); + + distance.at(Source) = 0; + std::priority_queue< + std::pair, + std::vector< std::pair >, + std::greater > + > q; + q.push( {0,Source} ); + + while(!q.empty()) + { + const auto [dist,node] = q.top(); + q.pop(); + + if(visited.at(node)) + continue; + + visited[node]=true; + + for(auto e: g.out_arcs(node)) + if( valid_arc(e) ) + { + auto [a,b] = g.arc_ends(e); + + if(weight.at(e)<0) + throw std::runtime_error( + "shortestPath_Dijkstra::solve found a negative edge"); + + value_type dnew = dist + weight.at(e); + if(distance.at(b)>dnew) + { + distance[b] = dnew; + parent[b] = e; + q.push({dnew,b}); + } + } + } + } + }; +} diff --git a/subprojects/MinCostFlow/include/mincostflow/vectorized_map.hpp b/subprojects/MinCostFlow/include/mincostflow/vectorized_map.hpp new file mode 100644 index 0000000..6ec9071 --- /dev/null +++ b/subprojects/MinCostFlow/include/mincostflow/vectorized_map.hpp @@ -0,0 +1,259 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace ln +{ + template + class vectorized_map + /* + a map-like data structure + that holds elements in a vector, hence access is O(1) + new elements are added to the smallest available slot and the key is returned, + we don't choose the key associated to an element, but once the key is set it will + remain valid until the element is removed. + iterators point to key values, not data, + though data can be accessed through the key + */ + { + using size_t = std::size_t; + + static constexpr index_t MAX_IDX = std::numeric_limits::max(); + std::vector< bool > valid_flag; + std::vector< data_t > data; + std::set< index_t > free_slots; + + void check_invariants()const + { + assert(valid_flag.size()==data.size()); + } + + void free_space() + { + while(!data.empty() && !valid_flag.back()) + // eliminate unused elements from the back of the buffer + { + auto ptr = index_t{data.size()-1}; + assert(free_slots.find(ptr)!=free_slots.end()); + + free_slots.erase(ptr); + valid_flag.pop_back(); + data.pop_back(); + } + check_invariants(); + } + + class base_iterator + { + protected: + const vectorized_map& const_ref; + index_t pos; + + void check_invariants()const + { + bool is_valid = const_ref.is_valid(pos); + bool is_infty = size_t(pos)==MAX_IDX; + + assert(is_valid ^ is_infty); + } + + + void next_valid() + { + while(size_t(pos)=const_ref.capacity()) + pos = index_t{MAX_IDX}; + } + + public: + base_iterator(const vectorized_map& c, index_t x): + const_ref{c}, + pos{x} + {} + + }; + + class iterator : public base_iterator + { + using base_iterator::pos; + using base_iterator::next_valid; + using base_iterator::check_invariants; + using base_iterator::const_ref; + + public: + iterator(vectorized_map& c, index_t x): + base_iterator{c} + { + next_valid(); + check_invariants(); + } + + bool operator==(const iterator& that)const + { + return pos == that.pos; + } + bool operator != (const iterator& that)const + { + return !(pos==that.pos); + } + index_t operator * ()const + { + return pos; + } + + iterator& operator++() + { + if(!const_ref.is_valid(pos)) + return *this; + ++pos; + next_valid(); + check_invariants(); + return *this; + } + }; + class const_iterator : public base_iterator + { + using base_iterator::pos; + using base_iterator::next_valid; + using base_iterator::check_invariants; + using base_iterator::const_ref; + + public: + + const_iterator(const vectorized_map& c, index_t x): + base_iterator{c,x} + { + next_valid(); + check_invariants(); + } + + bool operator==(const const_iterator& that)const + { + return pos == that.pos; + } + bool operator != (const const_iterator& that)const + { + return !(pos==that.pos); + } + index_t operator * ()const + { + return pos; + } + + const_iterator& operator++() + { + if(!const_ref.is_valid(pos)) + return *this; + ++pos; + next_valid(); + check_invariants(); + return *this; + } + }; + + public: + bool is_valid(index_t pos)const + { + return size_t(pos) +#include + +#define CHECK(cond,mesg) \ + if(!(cond)) throw std::runtime_error(mesg); + +template +void test_case(const std::vector>& arcs, + const std::vector & capacity, + const int source, + const int sink, + const std::vector& sol) +{ + ln::digraph graph; + pathsolver_t solver; + + graph.add_node(source); + graph.add_node(sink); + + std::vector res_cap; + for(auto i=0UL;i +void test() +{ + { // case 1 + int source=0; + int sink = 1; + std::vector flow{1,0,0,0,0,0}; + + std::vector< std::pair > arc{{0,1},{0,2},{1,3},{1,2},{1,0},{3,1}}; + std::vector capacity{1,9,5,1,7,4}; + + test_case(arc,capacity,source,sink,flow); + } + { // case 2 + int source=0; + int sink = 1; + std::vector flow{1,2,0,1,2}; + + std::vector< std::pair > arc{{0,2},{0,3}, {3,2},{2,1},{3,1}}; + std::vector capacity{1,2,2,2,2}; + + test_case(arc,capacity,source,sink,flow); + } +} + +int main() +{ + using value_type = int; + + try + { + test< ln::maxflow_augmenting_path >(); + test< ln::maxflow_augmenting_path >(); + test< ln::maxflow_preflow >(); + test< ln::maxflow_scaling >(); + test< ln::maxflow_scaling >(); + }catch(...) + { + return 1; + } + return 0; +} + + diff --git a/subprojects/MinCostFlow/test/meson.build b/subprojects/MinCostFlow/test/meson.build new file mode 100644 index 0000000..89d1e07 --- /dev/null +++ b/subprojects/MinCostFlow/test/meson.build @@ -0,0 +1,15 @@ +test('Vectorized Map', + executable('vec-map','vec-map.cpp', + dependencies: dep_mincostflow)) + +test('Shortest Path', + executable('short-path','short-path.cpp', + dependencies: dep_mincostflow)) + +test('Max Flow', + executable('max-flow','max-flow.cpp', + dependencies: dep_mincostflow)) + +test('Min Cost Max Flow', + executable('min-cost-flow','min-cost-flow.cpp', + dependencies: dep_mincostflow)) diff --git a/subprojects/MinCostFlow/test/min-cost-flow.cpp b/subprojects/MinCostFlow/test/min-cost-flow.cpp new file mode 100644 index 0000000..bd8f3fb --- /dev/null +++ b/subprojects/MinCostFlow/test/min-cost-flow.cpp @@ -0,0 +1,274 @@ +#include "mincostflow/graph.hpp" +#include "mincostflow/maxflow.hpp" +#include "mincostflow/mincostflow.hpp" +#include +#include + +#define CHECK(cond,mesg) \ + if(!(cond)) throw std::runtime_error(mesg); + +template +void test_case(const std::vector>& arcs, + const std::vector & capacity, + const std::vector & weight, + const int source, + const int sink, + const std::vector& sol) +{ + ln::digraph graph; + pathsolver_t solver; + + graph.add_node(source); + graph.add_node(sink); + + std::vector res_cap; + std::vector res_cost; + for(auto i=0UL;i +void test() +{ + { // case 1 + int source=0; + int sink = 1; + std::vector flow{1,0,0,0,0,0}; + + std::vector< std::pair > arc{{0,1},{0,2},{1,3},{1,2},{1,0},{3,1}}; + std::vector capacity{1,9,5,1,7,4}; + std::vector weight{1,1,1,1,1,1}; + + test_case(arc,capacity,weight,source,sink,flow); + } + { // case 2 + int source=0; + int sink = 1; + std::vector flow{1,2,0,1,2}; + + std::vector< std::pair > arc{{0,2},{0,3}, {3,2},{2,1},{3,1}}; + std::vector capacity{1,2,2,2,2}; + std::vector weight{1,1,1,1,1,1}; + + test_case(arc,capacity,weight,source,sink,flow); + } + { + // case 3 + int source=0; + int sink=1; + std::vector flow{2,5,2,0,0}; + std::vector< std::pair > arc{{0,2},{0,1},{2,1},{1,3},{0,3}}; + std::vector capacity{2,5,7,8,6}; + std::vector weight{1,3,2,2,6}; + + test_case(arc,capacity,weight,source,sink,flow); + } + { + // case 4 + int source = 0; + int sink = 1; + std::vector flow{0,4,1,0,0,1,1,0}; + std::vector< std::pair > arc{{0,2},{0,1},{0,3},{1,3},{2,3},{2,1},{3,2},{3,0}}; + std::vector capacity{2,4,3,3,3,1,1,4}; + std::vector weight{2,3,1,0,2,0,0,4}; + + test_case(arc,capacity,weight,source,sink,flow); + } + { + // case 5 + int source = 0; + int sink = 1; + std::vector flow{1,1,0,0,1,2}; + std::vector< std::pair > arc{{0,3},{0,2},{1,2},{1,0},{2,3},{3,1}}; + std::vector capacity{2,1,1,1,4,2}; + std::vector weight{4,1,0,1,2,0}; + + test_case(arc,capacity,weight,source,sink,flow); + + } +} + +int main() +{ + using value_type = int; + + try + { + test< ln::mincostflow_EdmondsKarp< + value_type, + ln::shortestPath_FIFO > >(); + test< ln::mincostflow_EdmondsKarp< + value_type, + ln::shortestPath_BellmanFord > >(); + + // expected failure + // test< ln::mincostflow_EdmondsKarp< + // value_type, + // ln::shortestPath_Dijkstra > >(); + + + test< ln::mincostflow_PrimalDual< + ln::shortestPath_FIFO, + ln::maxflow_augmenting_path >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_FIFO, + ln::maxflow_augmenting_path >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_FIFO, + ln::maxflow_scaling >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_FIFO, + ln::maxflow_scaling >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_FIFO, + ln::maxflow_preflow >>(); + + + test< ln::mincostflow_PrimalDual< + ln::shortestPath_BellmanFord, + ln::maxflow_augmenting_path >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_BellmanFord, + ln::maxflow_augmenting_path >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_BellmanFord, + ln::maxflow_scaling >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_BellmanFord, + ln::maxflow_scaling >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_BellmanFord, + ln::maxflow_preflow >>(); + + + test< ln::mincostflow_PrimalDual< + ln::shortestPath_Dijkstra, + ln::maxflow_augmenting_path >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_Dijkstra, + ln::maxflow_augmenting_path >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_Dijkstra, + ln::maxflow_scaling >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_Dijkstra, + ln::maxflow_scaling >>(); + test< ln::mincostflow_PrimalDual< + ln::shortestPath_Dijkstra, + ln::maxflow_preflow >>(); + + + + + + test< ln::mincostflow_capacityScaling< + ln::shortestPath_FIFO, + ln::maxflow_augmenting_path + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_FIFO, + ln::maxflow_augmenting_path + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_FIFO, + ln::maxflow_scaling + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_FIFO, + ln::maxflow_scaling + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_FIFO, + ln::maxflow_preflow + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_BellmanFord, + ln::maxflow_augmenting_path + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_BellmanFord, + ln::maxflow_augmenting_path + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_BellmanFord, + ln::maxflow_scaling + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_BellmanFord, + ln::maxflow_scaling + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_BellmanFord, + ln::maxflow_preflow + >> (); + + test< ln::mincostflow_capacityScaling< + ln::shortestPath_Dijkstra, + ln::maxflow_augmenting_path + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_Dijkstra, + ln::maxflow_augmenting_path + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_Dijkstra, + ln::maxflow_scaling + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_Dijkstra, + ln::maxflow_scaling + >> (); + test< ln::mincostflow_capacityScaling< + ln::shortestPath_Dijkstra, + ln::maxflow_preflow + >> (); + + test< ln::mincostflow_costScaling< + ln::maxflow_augmenting_path + >> (); + test< ln::mincostflow_costScaling< + ln::maxflow_augmenting_path + >> (); + test< ln::mincostflow_costScaling< + ln::maxflow_scaling + >> (); + test< ln::mincostflow_costScaling< + ln::maxflow_scaling + >> (); + test< ln::mincostflow_costScaling< + ln::maxflow_preflow + >> (); + }catch(std::exception& e) + { + std::cerr << e.what() << '\n'; + return 1; + } + return 0; +} + + diff --git a/subprojects/MinCostFlow/test/short-path.cpp b/subprojects/MinCostFlow/test/short-path.cpp new file mode 100644 index 0000000..4650435 --- /dev/null +++ b/subprojects/MinCostFlow/test/short-path.cpp @@ -0,0 +1,86 @@ +#include "mincostflow/graph.hpp" +#include "mincostflow/shortestpath.hpp" +#include +#include + +#define CHECK(cond,mesg) \ + if(!(cond)) throw std::runtime_error(mesg); + +template +void test_case(const std::vector>& arcs, + const std::vector & length, + const int source, + const std::vector& sol) +{ + ln::digraph graph; + pathsolver_t solver; + + for(auto i=0UL;i weights; + for(auto i=0UL;i +void test() +{ + { // case 1 + int source=0; + std::vector distance{0,1,2,6}; + + std::vector< std::pair > arc{{0,1},{0,2},{1,3},{1,2},{1,0},{3,1}}; + std::vector length{1,9,5,1,7,4}; + + test_case(arc,length,source,distance); + } + { // case 2 + int source=0; + std::vector distance{0,4,11,9}; + + std::vector< std::pair > arc{{0,1},{1,3},{1,0},{1,2},{2,1},{3,2}}; + std::vector length{4,5,4,7,7,3}; + + test_case(arc,length,source,distance); + } +} + +int main() +{ + using length_type = int; + + try + { + test< ln::shortestPath_Dijkstra >(); + test< ln::shortestPath_FIFO >(); + test< ln::shortestPath_BellmanFord >(); + }catch(...) + { + return 1; + } + return 0; +} + diff --git a/subprojects/MinCostFlow/test/vec-map.cpp b/subprojects/MinCostFlow/test/vec-map.cpp new file mode 100644 index 0000000..b50e920 --- /dev/null +++ b/subprojects/MinCostFlow/test/vec-map.cpp @@ -0,0 +1,41 @@ +#include + +#define CHECK(cond) \ + if(!(cond)) return 1; + + +int main() +{ + ln::vectorized_map vm; + + // initial state + CHECK(vm.capacity()==0 && vm.size()==0); + + // add three elements + vm.insert(1); + vm.insert(2); + vm.insert(3); + CHECK(vm.capacity()==3 && vm.size()==3); + + // remove the first + vm.erase(0); + CHECK(vm.capacity()==3 && vm.size()==2); + + // remove twice, no effect + vm.erase(0); + CHECK(vm.capacity()==3 && vm.size()==2); + + // add a new one, to the first empty slot + vm.insert(11); + CHECK(vm.capacity()==3 && vm.size()==3); + + // remove a non-existent index, no effect + vm.erase(4); + CHECK(vm.capacity()==3 && vm.size()==3); + + // remove from the tail, effect capacity changes + vm.erase(1); + vm.erase(2); + CHECK(vm.capacity()==1 && vm.size()==1); + return 0; +} From 5d392d0f1e530e3aeddb3d4d7b2e444ec73c63fa Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Mon, 15 Aug 2022 09:37:15 +0200 Subject: [PATCH 05/20] Min-Cost flow using eduardo's library --- .../SyncSimulatedPaymentSession.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index 41c6b90..6e836d7 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -1,7 +1,7 @@ from .UncertaintyNetwork import UncertaintyNetwork from .OracleLightningNetwork import OracleLightningNetwork -from .MinCostFlow import MCFNetwork +from MinCostFlow import MCFNetwork from typing import List @@ -72,9 +72,14 @@ def _prepare_mcf_solver(self, src, dest, amt: int = 1, mu: int = 100_000_000, ba continue cnt = 0 # QUANTIZATION): - for capacity, cost in channel.get_piecewise_linearized_costs(mu=mu): - index = self._min_cost_flow.AddArc(self._mcf_id[s], - self._mcf_id[d], + direction = 0 + if s>d: + direction = 1 + for part,(capacity, cost) in enumerate(channel.get_piecewise_linearized_costs(mu=mu)): + index = self._min_cost_flow.AddArc(s, + d, + channel.short_channel_id, + direction + part*2, capacity, cost) self._arc_to_channel[index] = (s, d, channel, 0) @@ -84,15 +89,15 @@ def _prepare_mcf_solver(self, src, dest, amt: int = 1, mu: int = 100_000_000, ba # Add node supply to 0 for all nodes for i in self._uncertainty_network.network.nodes(): - self._min_cost_flow.SetNodeSupply(self._mcf_id[i], 0) + self._min_cost_flow.SetNodeSupply(i, 0) # add amount to sending node self._min_cost_flow.SetNodeSupply( - self._mcf_id[src], int(amt)) # /QUANTIZATION)) + src, int(amt)) # /QUANTIZATION)) # add -amount to recipient nods self._min_cost_flow.SetNodeSupply( - self._mcf_id[dest], -int(amt)) # /QUANTIZATION)) + dest, -int(amt)) # /QUANTIZATION)) def _next_hop(self, path): """ From 1f4671e3d2871fc7fc981bc36a9832529b3db4b2 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Mon, 15 Aug 2022 10:55:24 +0200 Subject: [PATCH 06/20] removed module --- .gitmodules | 4 - subprojects/MinCostFlow/LICENSE | 190 ------- .../MinCostFlow/benchmark/benchmark-mcf.cpp | 299 ---------- subprojects/MinCostFlow/benchmark/gen.py | 105 ---- subprojects/MinCostFlow/benchmark/meson.build | 4 - subprojects/MinCostFlow/benchmark/run.sh | 20 - subprojects/MinCostFlow/c-api.md | 109 ---- subprojects/MinCostFlow/doc/biblio.bib | 21 - subprojects/MinCostFlow/doc/info.tex | 56 -- subprojects/MinCostFlow/doc/intro.tex | 36 -- subprojects/MinCostFlow/doc/main.tex | 47 -- subprojects/MinCostFlow/doc/makefile | 17 - subprojects/MinCostFlow/doc/proposal.tex | 47 -- subprojects/MinCostFlow/doc/schedule.tex | 48 -- subprojects/MinCostFlow/documentation.md | 78 --- .../MinCostFlow/examples/kattis-maxflow.cpp | 72 --- .../examples/kattis-mincostmaxflow.cpp | 76 --- .../examples/kattis-shortestpath1.cpp | 64 --- subprojects/MinCostFlow/examples/meson.build | 8 - subprojects/MinCostFlow/include/meson.build | 2 - .../MinCostFlow/include/mincostflow/graph.hpp | 412 -------------- .../include/mincostflow/maxflow.hpp | 289 ---------- .../include/mincostflow/meson.build | 1 - .../include/mincostflow/mincostflow.hpp | 479 ---------------- .../include/mincostflow/scope_guard.hpp | 45 -- .../include/mincostflow/shortestpath.hpp | 518 ------------------ .../include/mincostflow/vectorized_map.hpp | 259 --------- subprojects/MinCostFlow/meson.build | 21 - subprojects/MinCostFlow/readme.md | 18 - subprojects/MinCostFlow/test/max-flow.cpp | 90 --- subprojects/MinCostFlow/test/meson.build | 15 - .../MinCostFlow/test/min-cost-flow.cpp | 274 --------- subprojects/MinCostFlow/test/short-path.cpp | 86 --- subprojects/MinCostFlow/test/vec-map.cpp | 41 -- 34 files changed, 3851 deletions(-) delete mode 100644 .gitmodules delete mode 100644 subprojects/MinCostFlow/LICENSE delete mode 100644 subprojects/MinCostFlow/benchmark/benchmark-mcf.cpp delete mode 100644 subprojects/MinCostFlow/benchmark/gen.py delete mode 100644 subprojects/MinCostFlow/benchmark/meson.build delete mode 100644 subprojects/MinCostFlow/benchmark/run.sh delete mode 100644 subprojects/MinCostFlow/c-api.md delete mode 100644 subprojects/MinCostFlow/doc/biblio.bib delete mode 100644 subprojects/MinCostFlow/doc/info.tex delete mode 100644 subprojects/MinCostFlow/doc/intro.tex delete mode 100644 subprojects/MinCostFlow/doc/main.tex delete mode 100644 subprojects/MinCostFlow/doc/makefile delete mode 100644 subprojects/MinCostFlow/doc/proposal.tex delete mode 100644 subprojects/MinCostFlow/doc/schedule.tex delete mode 100644 subprojects/MinCostFlow/documentation.md delete mode 100644 subprojects/MinCostFlow/examples/kattis-maxflow.cpp delete mode 100644 subprojects/MinCostFlow/examples/kattis-mincostmaxflow.cpp delete mode 100644 subprojects/MinCostFlow/examples/kattis-shortestpath1.cpp delete mode 100644 subprojects/MinCostFlow/examples/meson.build delete mode 100644 subprojects/MinCostFlow/include/meson.build delete mode 100644 subprojects/MinCostFlow/include/mincostflow/graph.hpp delete mode 100644 subprojects/MinCostFlow/include/mincostflow/maxflow.hpp delete mode 100644 subprojects/MinCostFlow/include/mincostflow/meson.build delete mode 100644 subprojects/MinCostFlow/include/mincostflow/mincostflow.hpp delete mode 100644 subprojects/MinCostFlow/include/mincostflow/scope_guard.hpp delete mode 100644 subprojects/MinCostFlow/include/mincostflow/shortestpath.hpp delete mode 100644 subprojects/MinCostFlow/include/mincostflow/vectorized_map.hpp delete mode 100644 subprojects/MinCostFlow/meson.build delete mode 100644 subprojects/MinCostFlow/readme.md delete mode 100644 subprojects/MinCostFlow/test/max-flow.cpp delete mode 100644 subprojects/MinCostFlow/test/meson.build delete mode 100644 subprojects/MinCostFlow/test/min-cost-flow.cpp delete mode 100644 subprojects/MinCostFlow/test/short-path.cpp delete mode 100644 subprojects/MinCostFlow/test/vec-map.cpp diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index e799283..0000000 --- a/.gitmodules +++ /dev/null @@ -1,4 +0,0 @@ -[submodule "subprojects/MinCostFlow"] - path = subprojects/MinCostFlow - url = https://github.com/Lagrang3/mincostflow.git - branch = master diff --git a/subprojects/MinCostFlow/LICENSE b/subprojects/MinCostFlow/LICENSE deleted file mode 100644 index 5c67b5f..0000000 --- a/subprojects/MinCostFlow/LICENSE +++ /dev/null @@ -1,190 +0,0 @@ - Apache License - Version 2.0, January 2004 - http://www.apache.org/licenses/ - - TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION - - 1. Definitions. - - "License" shall mean the terms and conditions for use, reproduction, - and distribution as defined by Sections 1 through 9 of this document. - - "Licensor" shall mean the copyright owner or entity authorized by - the copyright owner that is granting the License. - - "Legal Entity" shall mean the union of the acting entity and all - other entities that control, are controlled by, or are under common - control with that entity. For the purposes of this definition, - "control" means (i) the power, direct or indirect, to cause the - direction or management of such entity, whether by contract or - otherwise, or (ii) ownership of fifty percent (50%) or more of the - outstanding shares, or (iii) beneficial ownership of such entity. - - "You" (or "Your") shall mean an individual or Legal Entity - exercising permissions granted by this License. - - "Source" form shall mean the preferred form for making modifications, - including but not limited to software source code, documentation - source, and configuration files. - - "Object" form shall mean any form resulting from mechanical - transformation or translation of a Source form, including but - not limited to compiled object code, generated documentation, - and conversions to other media types. - - "Work" shall mean the work of authorship, whether in Source or - Object form, made available under the License, as indicated by a - copyright notice that is included in or attached to the work - (an example is provided in the Appendix below). - - "Derivative Works" shall mean any work, whether in Source or Object - form, that is based on (or derived from) the Work and for which the - editorial revisions, annotations, elaborations, or other modifications - represent, as a whole, an original work of authorship. For the purposes - of this License, Derivative Works shall not include works that remain - separable from, or merely link (or bind by name) to the interfaces of, - the Work and Derivative Works thereof. - - "Contribution" shall mean any work of authorship, including - the original version of the Work and any modifications or additions - to that Work or Derivative Works thereof, that is intentionally - submitted to Licensor for inclusion in the Work by the copyright owner - or by an individual or Legal Entity authorized to submit on behalf of - the copyright owner. For the purposes of this definition, "submitted" - means any form of electronic, verbal, or written communication sent - to the Licensor or its representatives, including but not limited to - communication on electronic mailing lists, source code control systems, - and issue tracking systems that are managed by, or on behalf of, the - Licensor for the purpose of discussing and improving the Work, but - excluding communication that is conspicuously marked or otherwise - designated in writing by the copyright owner as "Not a Contribution." - - "Contributor" shall mean Licensor and any individual or Legal Entity - on behalf of whom a Contribution has been received by Licensor and - subsequently incorporated within the Work. - - 2. Grant of Copyright License. Subject to the terms and conditions of - this License, each Contributor hereby grants to You a perpetual, - worldwide, non-exclusive, no-charge, royalty-free, irrevocable - copyright license to reproduce, prepare Derivative Works of, - publicly display, publicly perform, sublicense, and distribute the - Work and such Derivative Works in Source or Object form. - - 3. Grant of Patent License. Subject to the terms and conditions of - this License, each Contributor hereby grants to You a perpetual, - worldwide, non-exclusive, no-charge, royalty-free, irrevocable - (except as stated in this section) patent license to make, have made, - use, offer to sell, sell, import, and otherwise transfer the Work, - where such license applies only to those patent claims licensable - by such Contributor that are necessarily infringed by their - Contribution(s) alone or by combination of their Contribution(s) - with the Work to which such Contribution(s) was submitted. If You - institute patent litigation against any entity (including a - cross-claim or counterclaim in a lawsuit) alleging that the Work - or a Contribution incorporated within the Work constitutes direct - or contributory patent infringement, then any patent licenses - granted to You under this License for that Work shall terminate - as of the date such litigation is filed. - - 4. Redistribution. You may reproduce and distribute copies of the - Work or Derivative Works thereof in any medium, with or without - modifications, and in Source or Object form, provided that You - meet the following conditions: - - (a) You must give any other recipients of the Work or - Derivative Works a copy of this License; and - - (b) You must cause any modified files to carry prominent notices - stating that You changed the files; and - - (c) You must retain, in the Source form of any Derivative Works - that You distribute, all copyright, patent, trademark, and - attribution notices from the Source form of the Work, - excluding those notices that do not pertain to any part of - the Derivative Works; and - - (d) If the Work includes a "NOTICE" text file as part of its - distribution, then any Derivative Works that You distribute must - include a readable copy of the attribution notices contained - within such NOTICE file, excluding those notices that do not - pertain to any part of the Derivative Works, in at least one - of the following places: within a NOTICE text file distributed - as part of the Derivative Works; within the Source form or - documentation, if provided along with the Derivative Works; or, - within a display generated by the Derivative Works, if and - wherever such third-party notices normally appear. The contents - of the NOTICE file are for informational purposes only and - do not modify the License. You may add Your own attribution - notices within Derivative Works that You distribute, alongside - or as an addendum to the NOTICE text from the Work, provided - that such additional attribution notices cannot be construed - as modifying the License. - - You may add Your own copyright statement to Your modifications and - may provide additional or different license terms and conditions - for use, reproduction, or distribution of Your modifications, or - for any such Derivative Works as a whole, provided Your use, - reproduction, and distribution of the Work otherwise complies with - the conditions stated in this License. - - 5. Submission of Contributions. Unless You explicitly state otherwise, - any Contribution intentionally submitted for inclusion in the Work - by You to the Licensor shall be under the terms and conditions of - this License, without any additional terms or conditions. - Notwithstanding the above, nothing herein shall supersede or modify - the terms of any separate license agreement you may have executed - with Licensor regarding such Contributions. - - 6. Trademarks. This License does not grant permission to use the trade - names, trademarks, service marks, or product names of the Licensor, - except as required for reasonable and customary use in describing the - origin of the Work and reproducing the content of the NOTICE file. - - 7. Disclaimer of Warranty. Unless required by applicable law or - agreed to in writing, Licensor provides the Work (and each - Contributor provides its Contributions) on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or - implied, including, without limitation, any warranties or conditions - of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A - PARTICULAR PURPOSE. You are solely responsible for determining the - appropriateness of using or redistributing the Work and assume any - risks associated with Your exercise of permissions under this License. - - 8. Limitation of Liability. In no event and under no legal theory, - whether in tort (including negligence), contract, or otherwise, - unless required by applicable law (such as deliberate and grossly - negligent acts) or agreed to in writing, shall any Contributor be - liable to You for damages, including any direct, indirect, special, - incidental, or consequential damages of any character arising as a - result of this License or out of the use or inability to use the - Work (including but not limited to damages for loss of goodwill, - work stoppage, computer failure or malfunction, or any and all - other commercial damages or losses), even if such Contributor - has been advised of the possibility of such damages. - - 9. Accepting Warranty or Additional Liability. While redistributing - the Work or Derivative Works thereof, You may choose to offer, - and charge a fee for, acceptance of support, warranty, indemnity, - or other liability obligations and/or rights consistent with this - License. However, in accepting such obligations, You may act only - on Your own behalf and on Your sole responsibility, not on behalf - of any other Contributor, and only if You agree to indemnify, - defend, and hold each Contributor harmless for any liability - incurred by, or claims asserted against, such Contributor by reason - of your accepting any such warranty or additional liability. - - END OF TERMS AND CONDITIONS - - Copyright 2022 Eduardo Quintana Miranda - - Licensed under the Apache License, Version 2.0 (the "License"); - you may not use this file except in compliance with the License. - You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. diff --git a/subprojects/MinCostFlow/benchmark/benchmark-mcf.cpp b/subprojects/MinCostFlow/benchmark/benchmark-mcf.cpp deleted file mode 100644 index def00a7..0000000 --- a/subprojects/MinCostFlow/benchmark/benchmark-mcf.cpp +++ /dev/null @@ -1,299 +0,0 @@ -#include -#include -#include -#include -#include - -typedef long long value_type; -typedef int nodeID_type; -typedef int arcID_type; -typedef ln::digraph graph_type; - -// typedef ln::maxflow_augmenting_path maxflow_t; -// typedef ln::maxflow_augmenting_path maxflow_t; -// typedef ln::maxflow_scaling maxflow_t; -// typedef ln::maxflow_scaling maxflow_t; -// typedef ln::maxflow_preflow maxflow_t; - -//typedef ln::mincostflow_EdmondsKarp< -// value_type,ln::shortestPath_FIFO> mincostflow_t; // 2.18s -// typedef ln::mincostflow_EdmondsKarp< -// value_type,ln::shortestPath_BellmanFord> mincostflow_t; // TLE -// typedef ln::mincostflow_EdmondsKarp mincostflow_t; // expected failure - -// typedef ln::mincostflow_PrimalDual< -// ln::shortestPath_FIFO,maxflow_t> mincostflow_t; // 2.54s, 3.02s, 2.75s, 3.18s, TLE -// typedef ln::mincostflow_PrimalDual< -// ln::shortestPath_BellmanFord,maxflow_t> mincostflow_t; // TLE -//typedef ln::mincostflow_PrimalDual< -// ln::shortestPath_Dijkstra,maxflow_t> mincostflow_t; -// 1.34s, 1.76s, 1.52s, 2.42s, TLE - -// typedef ln::mincostflow_capacityScaling< -// ln::shortestPath_Dijkstra,maxflow_t> mincostflow_t; -// TLE - -struct integer_graph -{ - graph_type g; - std::vector capacity; - std::vector weight; - - integer_graph(int N_nodes) - { - for(nodeID_type i=0;i& cap) - { - capacity.resize(g.max_num_arcs()); - for(arcID_type i=0;i& cost) - { - weight.resize(g.max_num_arcs()); - for(arcID_type i=0;i solve_ortools( - const int N_nodes, - nodeID_type S, nodeID_type T, - const std::vector>& edges, - const std::vector& cap, - const std::vector& cost, - const std::string tname) -{ - auto start = std::chrono::high_resolution_clock::now(); - operations_research::SimpleMaxFlow max_flow; - operations_research::SimpleMinCostFlow mincost_flow; - - for(int i=0;i(stop-start) .count() - << std::endl; - - return {Flow,Cost}; -} - - -template -void solve(integer_graph& G, - const std::vector& capacity, - const std::vector& cost, - nodeID_type S, - nodeID_type T, - std::string tname) -{ - G.set_capacity(capacity); - G.set_cost(cost); - - auto& graph = G.g; - - auto start = std::chrono::high_resolution_clock::now(); - mincostflow_t f; - f.solve(graph,graph.get_node(S),graph.get_node(T),G.weight,G.capacity); - auto stop = std::chrono::high_resolution_clock::now(); - - std::cout << tname << " " << - std::chrono::duration_cast(stop-start) .count() - << std::endl; -} - -std::pair -check_constraints( - const integer_graph& G, - const std::vector>& ed_list, - const std::vector& capacity, - const std::vector& cost, - const nodeID_type S, - const nodeID_type T, - const int N_nodes, - const int N_edges) -{ - // check capacity constraint - for(int e=0;e=f && f>=0); - } - - // check balances - std::vector balance(N_nodes,0); - - for(int e=0;e=0); - assert(balance[S] == -balance[T]); - - value_type Cost = 0; - for(int e=0;e, - ln::maxflow_augmenting_path - > - >(G,capacity,weight,S,T,"Primal-dual"); - auto [flow,cost] = check_constraints(G,ed_list,capacity,weight,S,T,N,M); - - assert(flow_0==flow && cost_0==cost); - } - - // { - // solve< - // ln::mincostflow_capacityScaling< - // ln::shortestPath_Dijkstra, - // ln::maxflow_augmenting_path - // > - // >(G,capacity,weight,S,T,"Capacity-scaling"); - // auto [flow,cost] = check_constraints(G,ed_list,capacity,weight,S,T,N,M); - // assert(flow_0==flow && cost_0==cost); - // } - { - solve< - ln::mincostflow_costScaling< - ln::maxflow_augmenting_path - > - >(G,capacity,weight,S,T,"Cost-scaling"); - auto [flow,cost] = check_constraints(G,ed_list,capacity,weight,S,T,N,M); - assert(flow_0==flow && cost_0==cost); - } - { - auto [flow,cost] = solve_ortools(N,S,T,ed_list,capacity,weight,"Ortools"); - assert(flow_0==flow && cost_0==cost); - } - return 0; -} - diff --git a/subprojects/MinCostFlow/benchmark/gen.py b/subprojects/MinCostFlow/benchmark/gen.py deleted file mode 100644 index facb244..0000000 --- a/subprojects/MinCostFlow/benchmark/gen.py +++ /dev/null @@ -1,105 +0,0 @@ -import networkx -import random -import sys,os -import tempfile -import numpy -import subprocess -import matplotlib.pyplot as plt -import json - -def generate(n,m,mcap,mwei,f): - S=0 - T=1 - g=networkx.gnm_random_graph(n,m,directed=True) - print(n,m,S,T,file=f) - for a,b in g.edges(): - print(a,b,random.choice(range(mcap)),random.choice(range(mwei)),file=f) - -def plot_data(n_array,series,outname): - fig = plt.figure() - ax = fig.subplots() - for x in series: - ax.loglog(n_array,series[x],'-o',label=x) - ax.legend() - ax.set_xlabel('N nodes') - ax.set_ylabel('Time (micro seconds)') - ax.grid() - fig.savefig(outname) - plt.show() - -def save_data(n_array,series,outname): - alldata={'N' : n_array} - alldata['series']=dict() - for name in series: - alldata['series'][name]=[ int(x) for x in series[name] ] - json.dump(alldata,open(outname,'w')) - -def read_data(filename): - alldata = json.load(open(filename,'r')) - return alldata['N'],alldata['series'] - -def cat(fname): - f = open(fname,'r') - for l in f: - print(l,end='') - -class TmpFile(): - def __init__(self): - self.name = tempfile.mktemp() - def __enter__(self): - return self - def __exit__(self,exc_type,exc_value,exc_traceback): - os.remove(self.name) - -def run_benchmark(): - n_array = [ 2**i for i in range(7,15)] - m_array = [ int(n*7.5) for n in n_array ] - cost_arr = [ 200 ]*len(n_array) - cap_arr = [ 200 ]*len(n_array) - rep_arr = [20]*len(n_array) - - names = [ - 'Ortools', - 'Edmonds-Karp', - 'Primal-dual', - # 'Capacity-scaling', - 'Cost-scaling', - ] - - series = dict() - for x in names: - series[x] = numpy.zeros(len(n_array)) - - for i in range(len(n_array)): - n,m,cap,cost = n_array[i],m_array[i],cap_arr[i],cost_arr[i] - print("N nodes = ",n) - for j in range(rep_arr[i]): - print(" rep = ",j) - with TmpFile() as tmp: - with open(tmp.name,'w') as fin: - generate(n,m,cap,cost,fin) - # cat(tmp.name) - p = subprocess.run( - ['./benchmark/benchmark-mcf'], - stdin=open(tmp.name,'r'),capture_output=True) - - data = p.stdout.decode().split('\n') - if p.returncode!=0: - print(p.stderr.decode()) - print("on test case") - cat(tmp.name) - raise "benchmark failed" - for d in data: - try: - name,val = d.split() - val = int(val) - series[name][i] = max(series[name][i],val) - except: - pass - return n_array, series - -if __name__ == "__main__": - #n_array,series = read_data('latest.json') - n_array,series = run_benchmark() - save_data(n_array,series,'latest.json') - plot_data(n_array,series,'latest.png') diff --git a/subprojects/MinCostFlow/benchmark/meson.build b/subprojects/MinCostFlow/benchmark/meson.build deleted file mode 100644 index 59965f0..0000000 --- a/subprojects/MinCostFlow/benchmark/meson.build +++ /dev/null @@ -1,4 +0,0 @@ -ortools = dependency('ortools') - -executable('benchmark-mcf','benchmark-mcf.cpp', - dependencies: [dep_mincostflow,ortools]) diff --git a/subprojects/MinCostFlow/benchmark/run.sh b/subprojects/MinCostFlow/benchmark/run.sh deleted file mode 100644 index b2f23bf..0000000 --- a/subprojects/MinCostFlow/benchmark/run.sh +++ /dev/null @@ -1,20 +0,0 @@ -tmp=`mktemp` -exe=benchmark/benchmark-mcf -gen=../benchmark/gen.py - -# lightning -N=2000 -M=150000 - -N=800 -M=6000 - -Mcap=200 -Mcost=200 - -for i in `seq 10`; do - echo "" - python $gen $N $M $Mcap $Mcost > $tmp - $exe < $tmp - -done diff --git a/subprojects/MinCostFlow/c-api.md b/subprojects/MinCostFlow/c-api.md deleted file mode 100644 index ee002ef..0000000 --- a/subprojects/MinCostFlow/c-api.md +++ /dev/null @@ -1,109 +0,0 @@ -Graph representation -=== -The problem graph is represented by the type `lngraph`. -This data structure assumes that the graph is directed and there can be multiple edges connecting -the same pair of vertices. - -The graph API is edge centered. Edges can be added, modified or removed. -Edges are refered to by their unique id, which is 64bit number that the API assigns internally. -This assignment is unspecified. - -There is no need now to define an API to do the same thing for nodes, they are just added -automatically and possibly never removed from the graph. - -The constant values `lngraph_SUCCESS` and `lngraph_FAIL` are generally used as return values to -indicate the success or failure in a function call. - -Memory allocation ---- -Memory is allocated in the heap for a `lngraph` by calling -``` -lngraph* lngraph_new(); -``` -This function returns a pointer to a newly created `lngraph` in a valid state. -The returned value is `NULL` if the allocation or the initialization fails. - -Memory is freed by calling -``` -int lngraph_free(lngraph* g); -``` -No fail. You may assume this always return `lngraph_SUCCESS`. - -Getters ---- -Returns the number of edges. -``` -int lngraph_numEdges(lngraph* g); -``` -No fail. - - -Gets the list of edges sorted by their ids and their capacity and cost properties. -`edgeID`, `capacity` and `cost` need to be allocated to hold at least as many elements as edges in -the graph. -``` -int lngraph_edges(lngraph* g, int64_t * edgeID, int64_t *capacity, int64_t * cost); -``` -No fail. You may assume this always return `lngraph_SUCCESS`. - - -Setters ---- -Edges are added to the graph by calling the function -``` -int64_t lngraph_addEdge(lngraph* g, const char* nodeA_id, const char* nodeB_id, int64_t capacity, -int64_t cost); -``` -The function returns a `int64_t` that corresponds to a unique identifier for the edge in the -database. -Nodes are added automatically if needed everytime `lngraph_addEdge` is called. -Possible errors: fails to add a new edge, the return value in that case is `lngraph_NOID`. -There will never be an edge whose id corresponds to `lngraph_NOID`. - -Edges can be removed by calling -``` -int lngraph_rmEdge(lngraph* g, int64_t edgeUniqueID); -``` -No fail. You may assume this always return `lngraph_SUCCESS`. - - -Min. Cost Flow solver -=== - -Node excess. A possitive value of the excess makes a node into a *source node* -and a negative value makes it a *sink node*. -``` -int lngraph_setExcess(lngraph* g, const char* nodeID, int64_t excess); -``` -If the node does not exist then it will be added automatically. -Possible errors: a new node is failed to be added, in that case the returned value is `lngraph_FAIL`. -Otherwise `lngraph_SUCCESS` is returned. - - -Sets the excess of every known node to 0. -``` -int lngraph_resetExcess(lngraph* g); -``` -No fail. You may assume this always return `lngraph_SUCCESS`. - -Solve the MCF problem with the excess previously set and returns the status. -``` -int lngraph_mincostflow(lngraph * g,int64_t * flow); -``` -The `flow` pointer points to a an allocated array of at least M elements, where M is the current -number of edges. Edges are ordered internally according to their ids and the `flow` value -corresponds to edges in that order. -Possible errors: it was not possible to find a feasible flow that satisfy the problem constraints, -the return value is `lngraph_FAIL`. -Otherwise `lngraph_SUCCESS` is returned. - -Solve the MCF problem where there is only a source and a sink node. -This is equivalent to setting the excess of the `source` to the value `excess` -and the excess of the `sink` to `-excess` and then calling `lngraph_mincostflow`. -``` -int lngraph_mincostflow_singleSinkSource(lngraph * g, - const char* source, - const char* sink, - int64_t excess, - int64_t * flow); -``` diff --git a/subprojects/MinCostFlow/doc/biblio.bib b/subprojects/MinCostFlow/doc/biblio.bib deleted file mode 100644 index a3ab84a..0000000 --- a/subprojects/MinCostFlow/doc/biblio.bib +++ /dev/null @@ -1,21 +0,0 @@ -@book{ahuja1993network, - title={Network Flows: Theory, Algorithms, and Applications}, - author={Ahuja, R.K. and Magnanti, T.L. and Orlin, J.B.}, - isbn={9780136175490}, - lccn={lc92026702}, - url={https://books.google.it/books?id=WnZRAAAAMAAJ}, - year={1993}, - publisher={Prentice Hall} -} - -@misc{pickhardt21, - doi = {10.48550/ARXIV.2107.05322}, - url = {https://arxiv.org/abs/2107.05322}, - author = {Pickhardt, Rene and Richter, Stefan}, - keywords = {Networking and Internet Architecture (cs.NI), FOS: Computer and information sciences, FOS: Computer and information sciences}, - title = {Optimally Reliable \& Cheap Payment Flows on the Lightning Network}, - publisher = {arXiv}, - year = {2021}, - copyright = {Creative Commons Attribution Share Alike 4.0 International} -} - diff --git a/subprojects/MinCostFlow/doc/info.tex b/subprojects/MinCostFlow/doc/info.tex deleted file mode 100644 index 6df4664..0000000 --- a/subprojects/MinCostFlow/doc/info.tex +++ /dev/null @@ -1,56 +0,0 @@ -\section{Biographical Information} - -\subsection*{Personal details} - -\begin{tabular}{ll} - Name:& Eduardo Quintana Miranda\\ - Affiliation:& University of Trieste, Astronomical Observatory of Trieste\\ - Course:& Astrophysics \\ - Degree Program:& PhD \\ - Country:& Italy\\ - E-mail:& \url{eduardo.quintana@pm.me} \\ - Github:& \url{github.com/Lagrang3} \\ - Discord:& \texttt{eduardo-qm} -\end{tabular} - -\subsection*{Educational background.} - -\begin{itemize} - \item 2008--2013. BSc. Nuclear Physics, University of Havana, - \item 2015--2018. MSc. Theoretical Physics, University of Trieste, - \item - 2018--2019. Master in High Performance Computing, International School of - Advanced Studies (SISSA), Trieste, - \item 2019--present. PhD in astrophysics, - University of Trieste. -\end{itemize} - -\subsection*{Programming background.} -My main programming language is C++, but I am also profiecient in C and python. - -I started programming back in 2009 during my 2nd year of BSc. The main drivers -were physics simulations and programming competitions. From 2009 until 2013 -I've participated intensively in those kind of competitions, -the most important were the -ACM-ICPC\footnote{\url{icpc.global}} competitions held every year, my team -managed to classify many times to the regional phase, which in our geographical -context where named \emph{Caribbean Finals of the ICPC}. Until this date I -participate once in a while on \url{codeforces.com} rounds by the username -\texttt{Lagrang3}\footnote{\url{https://codeforces.com/profile/Lagrang3}}. -I have a fair knowledge of a variety of classical algorithms somewhat seasoned -by my participation on programming competitions. I am familiar with shortest -path, max-flow and linear min-cost flow graph problems and textbook algorithms -for solving them. - -I've acquired some professional programming skills during the Master in High -Performance Computing. Advanced C++, Python, bash, git, and parallel programming -in OpenMP, MPI and Cuda where among the topics taught in that curriculum. - -In 2021 I've succesfully completed a Google Summer of Code project for the -proposal of FFT utilities in Boost Math -library\footnote{\url{https://github.com/BoostGSoC21/math-fft-report/releases/download/v1.1/gsoc-report.pdf}}. - -Today, I am doing a PhD in astrophysics working on the developement of a -simulation code -to study the effects of general relativity in the formation of cosmological -structures.\footnote{\url{https://github.com/Lagrang3/gevolution-1.2/tree/gev-api}} diff --git a/subprojects/MinCostFlow/doc/intro.tex b/subprojects/MinCostFlow/doc/intro.tex deleted file mode 100644 index 6c7d67b..0000000 --- a/subprojects/MinCostFlow/doc/intro.tex +++ /dev/null @@ -1,36 +0,0 @@ -\section{Introduction} - -Payments in the Lightning Network are performed along shortest path on the -network that try to minimize the transfer fees. -Because the participating nodes have incomplete information about the balance of -the individual channels, it is not always possible to perform the payments along -the shortest path, and the back-end of lightning has to work out alternative -routes until one is found with enough liquidity to process the payment. - -\cite{pickhardt21} proposes an alternative algorithm for payments, by -introducing the following novel ideas: -\begin{enumerate} - \item Each payment can be routed through a combination of different paths, - in what is called a \emph{Mult-Part Payment} (MPP). - The transfer of money on each channel is bounded by the channel's capacity - (in one specific direction) and the sum of the inbound and outbound payments - for the nodes satisfy balance constraints, the MPP problem can be solved - with the algorithms that find maximum flows on a directed network. - \item The arcs of the network can be characterized by a probability - distribution function for the realization of the transfer (the value of the - flow on the arc that can be processed). - The probability of success of the MPP becomes the multiplication of the - probability of success of every non-trivial flow on the network. - \item By defining a new function of the MPP flow, - which is the negative logarithm of that probability of success, - one has that the probability of an MPP payment is maximized when this new - cost function is minimized. And since the logarithm of the product is the - sum of logarithms, this cost function is separable in the sum of the cost - for the individual costs of the flow through the arcs of the network. - Thus the problem of maximizing the probability of success of a MPP payment - becomes a Min Cost Flow problem with non-linear costs on arcs. - \item By assuming some constraint in the probability function of the arcs, - one can reduce to a particular class of problems for which the cost of the - flow is a convex function. For this class of problem efficient polynomial - algorithms are known. -\end{enumerate} diff --git a/subprojects/MinCostFlow/doc/main.tex b/subprojects/MinCostFlow/doc/main.tex deleted file mode 100644 index 27a4e29..0000000 --- a/subprojects/MinCostFlow/doc/main.tex +++ /dev/null @@ -1,47 +0,0 @@ -\documentclass[11pt,a4paper]{article} - -\usepackage{amsmath} -\usepackage{amssymb} -\usepackage{amsthm} -\usepackage[a4paper,margin=1in]{geometry} -\usepackage{graphicx} -\usepackage{subcaption} -\usepackage{hyperref} -\usepackage{listings} -\usepackage{bold-extra} - -\hypersetup{ - colorlinks=false, - hidelinks=true -} -\bibliographystyle{plain} % FIXME: styles are not found - -\newcommand{\Ring}{\mathcal{R}} -\newcommand{\Order}{\mathcal{O}} -\newcommand{\Complex}{\mathbb{C}} -\newcommand{\Real}{\mathbb{R}} -\newcommand{\Integer}{\mathbb{Z}} - -\newcommand{\comment}[1]{\textrm{\textit{#1}}} -\newcommand{\function}[1]{\textbf{#1}} -\newcommand{\variable}[1]{\textit{#1}} - -\title{Summer of Bitcoin 2022:\\ Efficient Min Convex Cost Flow Solver for Lightning Network Payment Delivery} -\author{Eduardo Quintana Miranda} - - -\begin{document} - -\maketitle - -%\tableofcontents -\input intro -\input proposal -\input schedule - -\bibliography{biblio} - -\appendix -\input info - -\end{document} diff --git a/subprojects/MinCostFlow/doc/makefile b/subprojects/MinCostFlow/doc/makefile deleted file mode 100644 index 084ab10..0000000 --- a/subprojects/MinCostFlow/doc/makefile +++ /dev/null @@ -1,17 +0,0 @@ -VPATH:=images -TARGETS:=main.pdf -LATEXEC:=rubber --pdf main.tex -IMAGES:= -TEX:=main.tex biblio.bib info.tex intro.tex proposal.tex schedule.tex - -export - -default: $(TARGETS) - -main.pdf : $(TEX) $(IMAGES) - -$(LATEXEC) - -clean: - -$(LATEXEC) --clean - -.PHONY: clean default diff --git a/subprojects/MinCostFlow/doc/proposal.tex b/subprojects/MinCostFlow/doc/proposal.tex deleted file mode 100644 index e58b297..0000000 --- a/subprojects/MinCostFlow/doc/proposal.tex +++ /dev/null @@ -1,47 +0,0 @@ -\section{Proposal} - -% Based on the description of the problem in the introduction -This Summer of Bitcoin project, proposes the design and implementation of a C++ -mini-library with no external dependencies for solving the min cost flow problem -(MCF) with convex costs that can be used to optimize the probability of success -of multipart payments in the Lightning Network. - -Some time of the project will be allocated to study of the state of the art of -the problem. At the present we are aware of two algorithms that can be used to -solve the MCF with convex costs. -One that transforms approximates the cost function to a piecewise linear -function and splits the arc into several linear cost arcs accordingly (section -14.4 of \cite{ahuja1993network}). -And the so called \emph{Capacity} -algorithm described in section 14.5 of \cite{ahuja1993network}, -which is a generalization of a \emph{Excess Scaling} solver for the linear -MCF problem. -During the course of this project duration, we will evaluate which algorithm -will be best suited for the problem domain. - -Parallelization of the algorithms will be considered and discussed during the -development of the back-end. -We will propose to use standard library tools for that purpose -which are available in C++17. - -For the preparation of this Summer of Bitcoin project proposal, we have been -working on a proof-of-concept MCF C++ library\footnote{% -\url{https://github.com/Lagrang3/mincostflow}} publicly available on Github. -The API of this MCF library permits the representation of a directed graph -\texttt{digraph} in a space efficient data structure. -It also allows for the computation -shortest paths on a weighted network using Dijsktra's and Breadth-First-Search -algorithms, which serve as building blocks for flow solvers. -There are two different back-ends for the solution of the Max Flow problem: -using Edmond-Karp's \emph{Augmenting Paths} and Dinic's \emph{Push-Relabel}. -Also these constitute building blocks for the Min Cost Flow solvers. -The library contains a MCF solver based on the \emph{Successive Shortest Path} -algorithm presented in section 9.7 on \cite{ahuja1993network}. - -In the library's respository we provide examples of applications that we have also -used as test for correctedness and efficiency of our implementation. -These examples actually solve two problems from the online judge Kattis: -\texttt{maxflow}\footnote{\url{https://open.kattis.com/problems/maxflow}} -and -\texttt{mincostmaxflow}\footnote{\url{https://open.kattis.com/problems/mincostmaxflow}}, -and they do pass the test cases within the time and memory constraints. diff --git a/subprojects/MinCostFlow/doc/schedule.tex b/subprojects/MinCostFlow/doc/schedule.tex deleted file mode 100644 index 329c540..0000000 --- a/subprojects/MinCostFlow/doc/schedule.tex +++ /dev/null @@ -1,48 +0,0 @@ -\section{Timeline} - - \subsection*{Before May 9} - Meet the mentors, discuss with them the proposed timeline, - deliverables and share ideas. - - \subsection*{May 9 -- May 23} - Search in the literature for algorithms that solve the min-cost flow problem - with convex cost functions. - Follow the bitcoin and lightning seminars. - - \subsection*{May 23 -- June 6} - Design of a standalone C++ library for solving the min-cost flow problem. - Implementation of a textbook efficient solution to the linear min-cost flow - problem, the \emph{Excess Scaling} algorithm described in - \cite{ahuja1993network} section 7.9. - This linear MCF solver can be used to solve convex cost problems with a - tuneable linear approximation by decomposing the non-linear cost arcs into - several linear cost arcs. - Set up a simple set of unit tests for correctedness and CI for the library - repository. - - \subsection*{June 6 -- June 20} - Use the code to reproduce the results - of Pickhardt's \emph{Payment Simulation} Jupyter Notebook. - Implementation of a textbook solution: the capacity scaling algorithm described in - \cite{ahuja1993network} section 14.5. This algorithm is a generalization of - the \emph{Excess Scaling}, hence it shouldn't be too difficult to adapt - the linear MCF already implemented to the general convex cost case. - - \subsection*{June 20 -- July 4} - Benchmark the canditate solutions. - Optimize the code to find a good tradeoff of runtimes and probability of - success. - Parallelization of the time critical routines. - - \subsection*{July 4 -- July 18} - Implementation of a small C-lightning plugin to demonstrate the use of the - library. - - \subsection*{July 18 -- August 1} - Document the library and buffer time for unfinished issues. - - \subsection*{August 1 -- August 15} - Preparation of the final report and draft the research article. - - \subsection*{August 15 -- August 22} - Submission of the final project report. diff --git a/subprojects/MinCostFlow/documentation.md b/subprojects/MinCostFlow/documentation.md deleted file mode 100644 index e143c6b..0000000 --- a/subprojects/MinCostFlow/documentation.md +++ /dev/null @@ -1,78 +0,0 @@ -Graph representation -=== -``` - class digraph; -``` - -Path solvers -=== - -Finds a path in a directed graph with un-weighted edges. -``` - class pathSearch_BFS; -``` - -Finds a path in a directed graph with un-weighted edges. -This label-relabel search algorithm has meaninful state. -Figure 7.6 of Ahuja 93. -``` - class pathSearch_labeling; -``` - -Pseudo-polynomial generic path optimization. -``` - class shortestPath_FIFO; -``` - -Bellman-Ford path optimization -``` - class shortestPath_BellmanFord; -``` - -Dijkstra path optimization. -``` - class shortestPath_Dijkstra; -``` - -Max-flow -=== - -Generic augmenting path algorithm, template on the path finder algorithm. -Ahuja figure 6.12. -``` - template - class maxflow_augmenting_path; -``` - -Capacity scaling algorithm template on the path finder. -Ahuja figure 7.3. -``` - template - class maxflow_scaling; -``` - -Preflow-push algorithm. -Ahuja figure 7.12. -``` - class maxflow_preflow; -``` - -Min-Cost-Flow -=== -Min-cost-max-flow based on Edmonds-Karp augmenting path algorithm, it greedily -searches for the smallest cost routes. The path optimizer could be any shortest -path algorithm that deals with negative weights. -``` - template - class mincostflow_EdmondsKarp; -``` - -Min-cost-max-flow Primal Dual algorithm. -It uses a potential function to set the edges costs to zero and it pushes flow -along cost-zero edges. -It is a template on a path optimizer and a maxflow engines. -The path optimizer could be any shortest path algorithm including Dijkstra. -``` - template - class mincostflow_PrimalDual; -``` diff --git a/subprojects/MinCostFlow/examples/kattis-maxflow.cpp b/subprojects/MinCostFlow/examples/kattis-maxflow.cpp deleted file mode 100644 index ea5d066..0000000 --- a/subprojects/MinCostFlow/examples/kattis-maxflow.cpp +++ /dev/null @@ -1,72 +0,0 @@ -// https://open.kattis.com/problems/maxflow - -#include -#include -#include -#include - -using value_type = int; -// typedef ln::maxflow_augmenting_path maxflow_t; // 1.04s -// typedef ln::maxflow_augmenting_path maxflow_t; // 0.04s -// typedef ln::maxflow_scaling maxflow_t; // 0.42s -typedef ln::maxflow_scaling maxflow_t; // 0.07s -// typedef ln::maxflow_preflow maxflow_t; // 0.53s - -int main() -{ - int N,M,S,T; - std::cin >> N >> M >> S >> T; - - - - ln::digraph Graph; - std::vector capacity; - - for(int e=0;e>a>>b>>c; - auto [arc,arc2] = Graph.add_arc(a,b,e); - - if(capacity.size() 0 ? 1 : 0); - } - - - std::cout << N << ' ' << max_flow << ' ' << M_count << '\n'; - - for(int i=0;i -#include - -typedef long long value_type; - -// typedef ln::maxflow_augmenting_path maxflow_t; -typedef ln::maxflow_augmenting_path maxflow_t; -// typedef ln::maxflow_scaling maxflow_t; -// typedef ln::maxflow_scaling maxflow_t; -// typedef ln::maxflow_preflow maxflow_t; - -//typedef ln::mincostflow_EdmondsKarp< -// value_type,ln::shortestPath_FIFO> mincostflow_t; // 2.18s -// typedef ln::mincostflow_EdmondsKarp< -// value_type,ln::shortestPath_BellmanFord> mincostflow_t; // TLE -// typedef ln::mincostflow_EdmondsKarp mincostflow_t; // expected failure - -// typedef ln::mincostflow_PrimalDual< -// ln::shortestPath_FIFO,maxflow_t> mincostflow_t; // 2.54s, 3.02s, 2.75s, 3.18s, TLE -// typedef ln::mincostflow_PrimalDual< -// ln::shortestPath_BellmanFord,maxflow_t> mincostflow_t; // TLE -//typedef ln::mincostflow_PrimalDual< -// ln::shortestPath_Dijkstra,maxflow_t> mincostflow_t; -// 1.34s, 1.76s, 1.52s, 2.42s, TLE - -//typedef ln::mincostflow_capacityScaling< -// ln::shortestPath_Dijkstra,maxflow_t> mincostflow_t; -// TLE -typedef ln::mincostflow_costScaling< - maxflow_t> mincostflow_t; //0.37s, 0.06s, 0.10s, 0.07s, 0.15s - -int main() -{ - int N,M,S,T; - std::cin >> N >> M >> S >> T; - - ln::digraph G; - std::vector capacity; - std::vector weight; - - G.add_node(S); - G.add_node(T); - - for(int e=0;e>a>>b>>c>>w; - auto [arc,arc2] = G.add_arc(a,b,e); - - capacity.resize(G.max_num_arcs()); - weight.resize(G.max_num_arcs()); - - capacity[arc] = c; - capacity[arc2]=0; - - weight[arc] = w; - weight[arc2]=-w; - } - - mincostflow_t f; - - auto Flow = f.solve(G,G.get_node(S),G.get_node(T),weight,capacity); - long long Cost = 0; - - for(int e=0;e(weight[arc])*f.flow_at(G,arc,capacity); - } - - std::cout << Flow << " " << Cost << '\n'; - - return 0; -} - diff --git a/subprojects/MinCostFlow/examples/kattis-shortestpath1.cpp b/subprojects/MinCostFlow/examples/kattis-shortestpath1.cpp deleted file mode 100644 index 7da53aa..0000000 --- a/subprojects/MinCostFlow/examples/kattis-shortestpath1.cpp +++ /dev/null @@ -1,64 +0,0 @@ -// https://open.kattis.com/problems/shortestpath1 - -#include "mincostflow/graph.hpp" -#include "mincostflow/shortestpath.hpp" -#include - -using length_type = int; -typedef ln::shortestPath_Dijkstra pathsolver_t; // 0.28s -//typedef ln::shortestPath_FIFO pathsolver_t; // 1.14s -//typedef ln::shortestPath_BellmanFord pathsolver_t; // >3.0s - -int main() -{ - - while(1) - { - int N_vertex,N_edges,Q,S; - std::cin>>N_vertex>>N_edges>>Q>>S; - if(N_vertex==0)break; - - ln::digraph Graph; - std::vector weights(N_edges); - pathsolver_t solver; - - - for(int e=0;e> a>>b>>w; - auto [arc,arc2] = Graph.add_arc(a,b,e); - - if(weights.size()>v; - auto pos = Graph.get_node(v); - if(Graph.is_valid(pos) && solver.distance.at(pos) - -#include -#include -#include -#include -#include -#include -#include -#include - - -namespace ln -{ - - class digraph_types - { - public: - - // a lot more efficient if we can internally identify arcs and nodes by their unique - // position in the buffer array - typedef std::size_t pos_type; - static constexpr pos_type NONE = std::numeric_limits::max(); - - struct node_pos_t - { - pos_type x{NONE}; - bool operator<(const node_pos_t& that)const - { - return x < that.x; - } - bool operator==(const node_pos_t& that)const - { - return x==that.x; - } - - operator pos_type() const - { - return x; - } - node_pos_t& operator++() - { - ++x; - return *this; - } - node_pos_t operator++(int) - { - ++x; - return node_pos_t{x-1}; - } - }; - struct arc_pos_t - { - pos_type x{NONE}; - bool operator<(const arc_pos_t& that)const - { - return x < that.x; - } - bool operator==(const arc_pos_t& that)const - { - return x==that.x; - } - operator pos_type() const - { - return x; - } - arc_pos_t& operator++() - { - ++x; - return *this; - } - arc_pos_t operator++(int) - { - ++x; - return arc_pos_t{x-1}; - } - }; - - struct arc_data_t - { - node_pos_t a{NONE},b{NONE}; - arc_pos_t dual{NONE}; - }; - - struct node_data_t - { - std::vector out_arcs,in_arcs; - - void rm_arc(arc_pos_t arc) - { - auto rm_arc_vec=[arc](std::vector& V) - { - if (auto ptr = std::find(V.begin(),V.end(),arc); - ptr!=V.end()) - { - *ptr = V.back(); - V.pop_back(); - } - }; - rm_arc_vec(in_arcs); - rm_arc_vec(out_arcs); - } - void add_in_arc(arc_pos_t arc) - { - in_arcs.push_back(arc); - } - void add_out_arc(arc_pos_t arc) - { - out_arcs.push_back(arc); - } - }; - }; - - - // TODO: template on custom allocator - template - class digraph : public digraph_types - /* - Represents: a directed graph with dual arcs to simulate the residual network. - This data structure represents only the topological information. - Nodes and arcs have fixed positions. - */ - /* - ideally I would like to: - - digraph g; - - g.add_arc(node_a,node_b,arc_ab); // adds also the dual - - g::node_id_t n_1 = g.source_node(arc_x); - g::node_id_t n_2 = g.dest_node(arc_x); - - g::node_id_t n = g.add_node(node_factory); // generates a new node - - - - // an optimized interface allows to pass additional structure as arrays - g.max_nodes(); // size of the node array - g.max_arcs() ; // size of arc array - - for() - */ - { - vectorized_map my_arcs; - vectorized_map my_nodes; - - std::unordered_map arcs_htable; - std::vector arcs_ids; - std::vector arcs_ids_flag; - - std::unordered_map nodes_htable; - std::vector nodes_ids; - std::vector nodes_ids_flag; - - public: - - const auto& arcs()const - { - return my_arcs; - } - const auto& nodes()const - { - return my_nodes; - } - - bool is_valid(arc_pos_t arc)const - { - return my_arcs.is_valid(arc); - } - bool is_valid(node_pos_t node)const - { - return my_nodes.is_valid(node); - } - bool has_id(node_pos_t node)const - { - return nodes_ids_flag.at(node); - } - bool has_id(arc_pos_t arc)const - { - return arcs_ids_flag.at(arc); - } - - auto arc_ends(arc_pos_t arc)const - { - return std::pair{ - my_arcs.at(arc).a, - my_arcs.at(arc).b - }; - } - arc_pos_t arc_dual(arc_pos_t arc)const - { - return my_arcs.at(arc).dual; - } - - auto arc_ends_nocheck(arc_pos_t arc)const noexcept - { - return std::pair{ - my_arcs[arc].a, - my_arcs[arc].b - }; - } - arc_pos_t arc_dual_nocheck(arc_pos_t arc)const noexcept - { - return my_arcs[arc].dual; - } - - void erase(arc_pos_t arc) - { - if(! is_valid(arc)) - return; - - auto [a,b] = arc_ends(arc); - - my_nodes.at(a).rm_arc(arc); - my_nodes.at(b).rm_arc(arc); - - if(has_id(arc)) - { - auto id = arcs_ids.at(arc); - arcs_htable.erase(id); - } - - my_arcs.erase(arc); - arcs_ids.resize(my_arcs.capacity()); - arcs_ids_flag.resize(my_arcs.capacity()); - } - void erase(node_pos_t node) - { - if(! is_valid(node)) - return; - std::vector ls_arcs; - std::copy(my_nodes.at(node).in_arcs.begin(), - my_nodes.at(node).in_arcs.end(), - std::back_inserter(ls_arcs)); - std::copy(my_nodes.at(node).out_arcs.begin(), - my_nodes.at(node).out_arcs.end(), - std::back_inserter(ls_arcs)); - - // first remove all incoming and outgoin arcs - for(auto arc: ls_arcs) - erase(arc); - - if(has_id(node)) - { - auto id = nodes_ids.at(node); - nodes_htable.erase(id); - } - my_nodes.erase(node); - nodes_ids.resize(my_nodes.capacity()); - nodes_ids_flag.resize(my_nodes.capacity()); - } - node_pos_t new_node() - { - node_pos_t node = my_nodes.insert(node_data_t{}); - - nodes_ids.resize(my_nodes.capacity()); - nodes_ids_flag.resize(my_nodes.capacity()); - - nodes_ids_flag.at(node)=false; - - return node; - } - - const auto& out_arcs(node_pos_t node)const - { - if(!is_valid(node)) - throw std::runtime_error( - "digraph::out_arcs invalid node"); - return my_nodes.at(node).out_arcs; - } - const auto& in_arcs(node_pos_t node)const - { - if(!is_valid(node)) - throw std::runtime_error( - "digraph::in_arcs invalid node"); - return my_nodes.at(node).in_arcs; - } - - arc_pos_t new_arc(node_pos_t a, node_pos_t b) - { - if(!is_valid(a) || !is_valid(b)) - throw std::runtime_error("digraph::new_arc add a new arc with invalid end nodes"); - - arc_pos_t arc = my_arcs.insert(arc_data_t{a,b,NONE}); - arcs_ids.resize(my_arcs.capacity()); - arcs_ids_flag.resize(my_arcs.capacity()); - - arcs_ids_flag.at(arc)=false; - - my_nodes.at(a).add_out_arc(arc); - my_nodes.at(b).add_in_arc(arc); - return arc; - } - void set_dual(arc_pos_t arc1, arc_pos_t arc2) - { - if(!is_valid(arc1) || !is_valid(arc2)) - throw std::runtime_error("digraph::set_dual invalid arcs"); - - my_arcs.at(arc1).dual = arc2; - my_arcs.at(arc2).dual = arc1; - } - - auto max_num_arcs()const - { - return my_arcs.capacity(); - } - auto num_arcs()const - { - return my_arcs.size(); - } - auto max_num_nodes()const - { - return my_nodes.capacity(); - } - auto num_nodes()const - { - return my_nodes.size(); - } - - // translation - node_id_t get_node_id(node_pos_t node)const - { - if(!is_valid(node)) - throw std::runtime_error("digraph::get_node_id invalid node"); - if(!has_id(node)) - throw std::runtime_error("digraph::get_node_id node without id"); - return nodes_ids.at(node); - } - arc_id_t get_arc_id(arc_pos_t arc)const - { - if(!is_valid(arc)) - throw std::runtime_error("digraph::get_arc_id invalid arc"); - if(!has_id(arc)) - throw std::runtime_error("digraph::get_arc_id arc without id"); - return arcs_ids.at(arc); - } - node_pos_t get_node(node_id_t id)const - { - node_pos_t node{NONE}; - if(auto ptr = nodes_htable.find(id); - ptr != nodes_htable.end()) - { - node = ptr->second; - } - return node; - } - arc_pos_t get_arc(arc_id_t id)const - { - arc_pos_t arc{NONE}; - if(auto ptr = arcs_htable.find(id); - ptr != arcs_htable.end()) - { - arc = ptr->second; - } - return arc; - } - node_pos_t add_node(node_id_t id) - { - auto node = get_node(id); - if(!is_valid(node)) - { - node = new_node(); - nodes_ids.at(node) = id; - nodes_ids_flag.at(node) = true; - nodes_htable[id] = node; - } - return node; - } - std::pair add_arc(node_id_t a, node_id_t b, arc_id_t id) - { - auto n_a = add_node(a); - auto n_b = add_node(b); - - if(auto arc = get_arc(id); is_valid(arc)) - throw std::runtime_error("digraph::add_arc arc id already exists"); - - auto arc1 = new_arc(n_a,n_b); - auto arc2 = new_arc(n_b,n_a); - set_dual(arc1,arc2); - - arcs_ids.at(arc1) = id; - arcs_ids_flag.at(arc1) = true; - arcs_htable[id] = arc1; - - return {arc1,arc2}; - } - void remove_node(node_id_t id) - { - auto node = get_node(id); - if(!is_valid(node)) - return; - nodes_htable.erase(id); - erase(node); - } - void remove_arc(arc_id_t id) - { - auto arc = get_arc(id); - if(!is_valid(arc)) - return; - arcs_htable.erase(id); - - auto arc2 = arc_dual(arc); - erase(arc); - erase(arc2); - } - - digraph() - {} - }; -} diff --git a/subprojects/MinCostFlow/include/mincostflow/maxflow.hpp b/subprojects/MinCostFlow/include/mincostflow/maxflow.hpp deleted file mode 100644 index 0f9da7c..0000000 --- a/subprojects/MinCostFlow/include/mincostflow/maxflow.hpp +++ /dev/null @@ -1,289 +0,0 @@ -#pragma once - -#include - -namespace ln -{ - - template - class maxflow_base : public digraph_types - { - public: - using value_type = T; - static constexpr value_type INFINITY = std::numeric_limits::max(); - - template - value_type flow_at( - const graph_t& g, - const arc_pos_t e, - const std::vector& capacity) - { - auto e2 = g.arc_dual(e); - return capacity.at(e2.x); - } - }; - - template - class maxflow_augmenting_path : public maxflow_base - { - public: - using base_type = maxflow_base; - using value_type = typename base_type::value_type; - using node_pos_t = typename base_type::node_pos_t; - using arc_pos_t = typename base_type::arc_pos_t; - using base_type::flow_at; - using base_type::INFINITY; - - - - template - value_type solve( - const graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - std::vector& capacity, - condition_t valid_arc) - { - value_type sent=0; - path_solver_type path_solver; - - while(1) - { - bool found = path_solver.solve( - g, - Source,Dest, - [valid_arc,&capacity](arc_pos_t e) - { - return capacity.at(e)>0 && valid_arc(e); - }); - - - if(!found) - break; - - auto path = path_solver.get_path(g,Dest); - - value_type k = INFINITY; - for(auto e : path) - { - k = std::min(k,capacity.at(e)); - } - - for(auto e: path) - { - capacity.at(e) -= k; - capacity.at(g.arc_dual(e)) += k; - } - - sent += k; - } - return sent; - } - - maxflow_augmenting_path() - {} - }; - - template - class maxflow_scaling : public maxflow_base - { - public: - using base_type = maxflow_base; - using value_type = typename base_type::value_type; - using node_pos_t = typename base_type::node_pos_t; - using arc_pos_t = typename base_type::arc_pos_t; - using base_type::flow_at; - using base_type::INFINITY; - - template - value_type solve( - const graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - std::vector& residual_cap, - condition_t valid_arc) - // augmenting path - { - value_type sent=0; - path_solver_type search_algo; - - value_type cap_flow = 1; - for(auto e : g.out_arcs(Source)) - cap_flow = std::max(cap_flow,residual_cap.at(e)); - - cap_flow = lower_bound_power2(cap_flow); - - // int cycle=0; - for(;cap_flow>0;) - { - // cycle++; - // std::cerr << "augmenting path cycle: " << cycle << '\n'; - // std::cerr << "flow sent: " << sent << '\n'; - // std::cerr << "cap flow: " << cap_flow << '\n'; - - bool found = search_algo.solve( - g, - Source,Dest, - // edge is valid if - [this,valid_arc,cap_flow,&residual_cap](arc_pos_t e) - { - return residual_cap.at(e)>=cap_flow && valid_arc(e); - }); - - if(! found) - { - cap_flow/=2; - // std::cerr << "path not found!\n"; - search_algo.reset(); - continue; - } - - auto path = search_algo.get_path(g,Dest); - - // std::cerr << "path found!\n"; - - for(auto e: path) - { - residual_cap[e] -= cap_flow; - residual_cap[g.arc_dual(e)] += cap_flow; - } - - sent += cap_flow; - } - return sent; - } - - maxflow_scaling() - {} - }; - - - template - class maxflow_preflow : public maxflow_base, public distance_structure - { - public: - using base_type = maxflow_base; - using value_type = typename base_type::value_type; - using node_pos_t = typename base_type::node_pos_t; - using arc_pos_t = typename base_type::arc_pos_t; - using base_type::flow_at; - using base_type::INFINITY; - - static constexpr auto flow_INFINITY = base_type::INFINITY; - static constexpr auto dist_INFINITY = distance_structure::INFINITY; - - std::vector excess; - - template - void initialize_distance( - const graph_t& g, - const node_pos_t Dest, - condition_t valid_arc) - { - distance_structure::init(g); - - distance.at(Dest)=0; - - std::queue q; - q.push(Dest); - - while(!q.empty()) - { - auto n = q.front(); - q.pop(); - - for(auto e: g.in_arcs(n)) - if( valid_arc(e) ) - { - // assert b==n - auto [a,b] = g.arc_ends(e); - int dnew = distance[b] + 1; - - if(distance[a]==dist_INFINITY) - { - distance[a] = dnew; - q.push(a); - } - } - } - } - public: - - template - value_type solve( - const graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - std::vector& residual_cap, - condition_t valid_arc) - { - excess.resize(g.max_num_nodes()); - std::fill(excess.begin(),excess.end(),0); - - initialize_distance(g,Dest,valid_arc); - std::queue q; - - auto push = [&](arc_pos_t e) - { - auto [a,b] = g.arc_ends(e); - const auto delta = std::min(excess[a],residual_cap.at(e)); - residual_cap.at(e) -= delta; - residual_cap.at(g.arc_dual(e)) += delta; - - assert(delta>=0); - - excess.at(a) -= delta; - excess.at(b) += delta; - - if(delta>0 && excess.at(b)==delta) - q.push(b); - }; - - auto relabel = [&](node_pos_t v) - { - int hmin = dist_INFINITY; - for(auto e : g.out_arcs(v)) - if(valid_arc(e) && residual_cap.at(e)>0) - hmin = std::min(hmin,distance.at(g.arc_ends(e).second)); - if(hmin0) - { - auto b = g.arc_ends(e).second; - if(distance[a]== distance[b]+1) - push(e); - } - - if(excess.at(a)==0) - break; - - relabel(a); - } - }; - - excess.at(Source) = flow_INFINITY; - distance.at(Source) = g.num_nodes(); - - for(auto e : g.out_arcs(Source)) - if(valid_arc(e)) - push(e); - - while(!q.empty()) - { - auto node = q.front(); - q.pop(); - - if(node!=Dest && node!=Source) - discharge(node); - } - return excess.at(Dest); - } - maxflow_preflow() - {} - - }; -} diff --git a/subprojects/MinCostFlow/include/mincostflow/meson.build b/subprojects/MinCostFlow/include/mincostflow/meson.build deleted file mode 100644 index 42a02b3..0000000 --- a/subprojects/MinCostFlow/include/mincostflow/meson.build +++ /dev/null @@ -1 +0,0 @@ -my_headers = files(['graph.hpp','shortestpath.hpp','maxflow.hpp','mincostflow.hpp']) diff --git a/subprojects/MinCostFlow/include/mincostflow/mincostflow.hpp b/subprojects/MinCostFlow/include/mincostflow/mincostflow.hpp deleted file mode 100644 index afc1530..0000000 --- a/subprojects/MinCostFlow/include/mincostflow/mincostflow.hpp +++ /dev/null @@ -1,479 +0,0 @@ -#pragma once - -#include -#include - -namespace ln -{ - template - class mincostflow_EdmondsKarp : public maxflow_base - { - public: - using base_type = maxflow_base; - using value_type = typename base_type::value_type; - using node_pos_t = typename base_type::node_pos_t; - using arc_pos_t = typename base_type::arc_pos_t; - using base_type::flow_at; - using base_type::INFINITY; - - template - value_type solve( - const graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - const std::vector& weight, - std::vector& residual_cap - ) - // augmenting path - { - value_type sent =0 ; - path_optimizer_type path_opt; - - while(true) - { - path_opt.solve( - g, - Source, - weight, - // edge is valid if - [&residual_cap](arc_pos_t e){ - return residual_cap.at(e)>0; - }); - - if(! path_opt.is_reacheable(Dest)) - break; - - auto path = path_opt.get_path(g,Dest); - - value_type k = INFINITY; - for(auto e : path) - { - k = std::min(k,residual_cap.at(e)); - } - - for(auto e: path) - { - residual_cap[e] -= k; - residual_cap[g.arc_dual(e)] += k; - } - - sent += k; - } - return sent; - } - - mincostflow_EdmondsKarp() - {} - }; - - - template - class mincostflow_PrimalDual : public maxflow_type - { - public: - using base_type = maxflow_type; - using value_type = typename base_type::value_type; - using node_pos_t = typename base_type::node_pos_t; - using arc_pos_t = typename base_type::arc_pos_t; - using base_type::flow_at; - using base_type::INFINITY; - - - template - value_type solve( - const graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - const std::vector& weight, - std::vector& residual_cap - ) - { - std::vector reduced_weight = weight; - - value_type sent =0 ; - path_optimizer_type path_opt; - - while(true) - { - path_opt.solve( - g, - Source, - reduced_weight, - // edge is valid if - [&residual_cap](arc_pos_t e) -> bool - { - return residual_cap.at(e)>0; - }); - - if(! path_opt.is_reacheable(Dest)) - break; - - const auto& distance{path_opt.distance}; - - for(auto e : g.arcs()) - { - - auto [a,b] = g.arc_ends(e); - if(distance[a]bool - { - return reduced_weight[e]==0; - }); - - sent += F; - } - return sent; - } - - mincostflow_PrimalDual() - {} - }; - - template - class mincostflow_capacityScaling : public maxflow_type - { - public: - using base_type = maxflow_type; - using value_type = typename base_type::value_type; - using node_pos_t = typename base_type::node_pos_t; - using arc_pos_t = typename base_type::arc_pos_t; - using base_type::flow_at; - using base_type::INFINITY; - - template - value_type solve( - graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - const std::vector& weight, - std::vector& residual_cap) - { - - std::vector reduced_weight = weight; - value_type maxflow{0}; - - // find the max-flow-anycost - maxflow = maxflow_type::solve( - g,Source,Dest, - residual_cap, - [](arc_pos_t)->bool{return true;}); - - value_type cap_flow = lower_bound_power2(maxflow); - - std::vector excess(g.max_num_nodes(),0); - - std::vector weight_ex = weight; - - auto update_reduced_costs = - [&](const std::vector& potential) - { - for(auto e : g.arcs()) - { - auto [src,dst] = g.arc_ends(e); - auto p_src = potential.at(src), p_dst = potential.at(dst); - - p_src = p_src == INFINITY ? 0 : p_src; - p_dst = p_dst == INFINITY ? 0 : p_dst; - - weight_ex.at(e) += p_src - p_dst; - } - }; - - auto push_flow = - [&](arc_pos_t e,value_type delta) - { - // std::cerr << " push flow at " << e << " delta = " << delta << "\n"; - auto [src,dst] = g.arc_ends(e); - - residual_cap[e]-=delta; - residual_cap[g.arc_dual(e)]+=delta; - - excess.at(src) -= delta; - excess.at(dst) += delta; - - // std::cerr << "push " << delta << " over " << e << '\n'; - }; - - // auto report = - // [&]() - // { - // std::cerr << "residual cap + mod. costs\n"; - // for(auto e : g.arcs()) - // { - // std::cerr << " " << e << " -> " << residual_cap[e] << " " << weight_ex[e] << "\n"; - // } - // std::cerr << "potential + excess\n"; - // for(auto v : g.nodes()) - // { - // std::cerr << " " << v << " -> " << excess[v] << "\n"; - // } - // }; - // - // std::cerr << " maxflow = " << maxflow << "\n"; - - // int cycle=0; - for(;cap_flow>0;cap_flow/=2) - { - // cycle++; - // std::cerr << "cycle " << cycle << " cap_flow = " << cap_flow << '\n'; - // report(); - - // saturate edges with negative cost - for(auto e : g.arcs()) - while(residual_cap.at(e)>=cap_flow && weight_ex.at(e)<0) - { - push_flow(e,cap_flow); - } - - path_optimizer_type path_opt; - - // build S and T - std::set Sset,Tset; - for(auto v : g.nodes()) - { - if(excess.at(v)>=cap_flow) - Sset.insert(v); - if(excess.at(v)<=-cap_flow) - Tset.insert(v); - } - - const auto multi_source_node = g.new_node(); - excess.resize(g.max_num_nodes()); - excess.at(multi_source_node) = 0; - const Scope_guard rm_node = [&](){ g.erase(multi_source_node);}; - - - - for(auto v : Sset) - { - auto arc1 = g.new_arc(multi_source_node,v); - auto arc2 = g.new_arc(v,multi_source_node); - - g.set_dual(arc1,arc2); - - weight_ex.resize(g.max_num_arcs()); - residual_cap.resize(g.max_num_arcs()); - - weight_ex.at(arc1) = 0; - residual_cap.at(arc1) = excess.at(v); - - weight_ex.at(arc2) = 0; - residual_cap.at(arc2) = 0; - - excess.at(multi_source_node) += excess.at(v); - excess.at(v) = 0; - } - - const Scope_guard restore_excess = [&]() - { - for(auto e : g.out_arcs(multi_source_node)) - { - auto [src,dst] = g.arc_ends(e); - excess.at(dst) = residual_cap.at(e); - } - }; - - while(!Sset.empty() && !Tset.empty()) - { - path_opt.solve( - g, multi_source_node, - weight_ex, - [cap_flow,&residual_cap](arc_pos_t e)->bool - { - return residual_cap.at(e)>=cap_flow; - } - ); - - const auto& distance{path_opt.distance}; - - auto it = std::find_if(Tset.begin(),Tset.end(), - [&](node_pos_t v)->bool { - return distance.at(v) " << distance[v]<<"\n"; - // } - - update_reduced_costs(distance); - - auto path = path_opt.get_path(g,dst); - for(auto e: path) - { - // auto [src,dst] = g.arc_ends(e); - push_flow(e,cap_flow); - } - - if(excess.at(dst)>-cap_flow) - Tset.erase(dst); - } - } - - return maxflow; - } - - public: - mincostflow_capacityScaling() - {} - }; - - template - class mincostflow_costScaling : public maxflow_type - { - public: - using base_type = maxflow_type; - using value_type = typename base_type::value_type; - using node_pos_t = typename base_type::node_pos_t; - using arc_pos_t = typename base_type::arc_pos_t; - using base_type::flow_at; - using base_type::INFINITY; - - template - value_type solve( - graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - const std::vector& weight, - std::vector& residual_cap) - { - value_type maxflow{0}; - - // find the max-flow-anycost - maxflow = maxflow_type::solve( - g,Source,Dest, - residual_cap, - [](arc_pos_t)->bool{return true;}); - - // return maxflow; - - std::vector reduced_weight = weight; - std::vector potential(g.max_num_nodes(),0); - std::vector excess(g.max_num_nodes(),0); - - auto relabel = - [&](const node_pos_t x,value_type eps) - { - // std::cerr << " relabel " << x << " delta = " << eps << "\n"; - - potential[x] -= eps; - for(auto e : g.out_arcs(x)) - { - reduced_weight[e] -= eps; - } - for(auto e : g.in_arcs(x)) - { - reduced_weight[e] += eps; - } - }; - - auto push_flow = - [&](arc_pos_t e,value_type delta) - { - auto [src,dst] = g.arc_ends(e); - // std::cerr << " push flow at " << e - // << " (" << src << "," << dst << ")" - // << " delta = " << delta << "\n"; - - residual_cap[e]-=delta; - residual_cap[g.arc_dual(e)]+=delta; - - excess[src] -= delta; - excess[dst] += delta; - - // std::cerr << "push " << delta << " over " << e << '\n'; - }; - - value_type maxC = 0; - const int N = g.num_nodes(); - for(auto e : g.arcs()) - { - reduced_weight[e] *= N; - maxC = std::max(maxC,reduced_weight[e]); - } - maxC = lower_bound_power2(maxC); - - - // int cycle=0; - for(;maxC>0;maxC/=2){ - // cycle++; - // std::cerr << "cycle " << cycle << " eps = " << maxC << '\n'; - - // improve - for(auto e : g.arcs()) - { - if(reduced_weight[e]<0 && residual_cap[e]>0) - { - push_flow(e,residual_cap[e]); - } - // this also does the trick of making the flow = 0 for arcs with - // reduced_weight>0, - } - std::set active; - for(auto n : g.nodes()) - if(excess[n]>0) - { - active.insert(n); - // std::cerr << n << " is added to active\n"; - } - // int ct=0; - while(!active.empty()) - { - // ct++; - // if(ct>20)exit(1); - - auto t = *active.begin(); - - // std::cerr << t << " is active\n"; - - bool pushed = false; - - for(auto e : g.out_arcs(t)) - { - auto rw = reduced_weight[e]; - auto rc = residual_cap[e]; - if(rw<0 && rw>=-maxC && rc>0) - { - pushed = true; - auto [a,b] = g.arc_ends(e); - auto d = std::min(excess[a],rc); - - push_flow(e,d); - - if(excess[a]<=0) - active.erase(a); - if(excess[b]>0) - active.insert(b); - - break; - } - } - - if(!pushed) - relabel(t,maxC); - } - } - // std::cerr << "done\n"; - return maxflow; - } - - public: - mincostflow_costScaling() - {} - }; - -} diff --git a/subprojects/MinCostFlow/include/mincostflow/scope_guard.hpp b/subprojects/MinCostFlow/include/mincostflow/scope_guard.hpp deleted file mode 100644 index 3624bd8..0000000 --- a/subprojects/MinCostFlow/include/mincostflow/scope_guard.hpp +++ /dev/null @@ -1,45 +0,0 @@ -#pragma once - -/* - credits to - https://stackoverflow.com/questions/31365013/what-is-scopeguard-in-c -*/ - -#include // std::function -#include // std::move - -namespace ln -{ - class Non_copyable - { - private: - auto operator=( Non_copyable const& ) -> Non_copyable& = delete; - Non_copyable( Non_copyable const& ) = delete; - public: - auto operator=( Non_copyable&& ) -> Non_copyable& = default; - Non_copyable() = default; - Non_copyable( Non_copyable&& ) = default; - }; - - class Scope_guard - : public Non_copyable - { - private: - std::function cleanup_; - - public: - friend - void dismiss( Scope_guard& g ) { g.cleanup_ = []{}; } - - ~Scope_guard() { cleanup_(); } - - template< class Func > - Scope_guard( Func const& cleanup ) - : cleanup_( cleanup ) - {} - - Scope_guard( Scope_guard&& other ) - : cleanup_( std::move( other.cleanup_ ) ) - { dismiss( other ); } - }; -} diff --git a/subprojects/MinCostFlow/include/mincostflow/shortestpath.hpp b/subprojects/MinCostFlow/include/mincostflow/shortestpath.hpp deleted file mode 100644 index 127726b..0000000 --- a/subprojects/MinCostFlow/include/mincostflow/shortestpath.hpp +++ /dev/null @@ -1,518 +0,0 @@ -#pragma once - -#include - -#include -#include -#include -#include -#include -#include - - -namespace ln -{ - inline long long int lower_bound_power2(long long int n) - { - if(n<=2) return n; - while(n != (n & -n)) - n -= (n & -n); - return n; - } - - - class parent_structure : public digraph_types - { - public: - std::vector parent; - - bool has_parent(node_pos_t x)const - { - return parent.at(x)!=arc_pos_t{NONE}; - } - bool is_reacheable(node_pos_t x)const - { - return has_parent(x); - } - - template - void init(const graph_t& g) - { - parent.resize(g.max_num_nodes()); - std::fill(parent.begin(),parent.end(),arc_pos_t{NONE}); - } - - - template - auto get_path(const graph_t& g,node_pos_t last)const - { - std::vector path; - while(1) - { - auto e = parent.at(last); - if(!g.is_valid(e)) - break; - - path.push_back(e); - auto [a,b] = g.arc_ends(e); - last = a; - - } - return path; - } - }; - - template - class distance_structure - { - public: - using value_type = T; - static constexpr value_type INFINITY = std::numeric_limits::max(); - std::vector distance; - - template - void init(const graph_t& g) - { - distance.resize(g.max_num_nodes()); - std::fill(distance.begin(),distance.end(),INFINITY); - } - }; - - - class pathSearch_BFS : public parent_structure, public distance_structure - /* - Represents: weigthless path using BFS - Invariant: - - User interface: - Complexity: |E|+|V| - */ - { - public: - using value_type = int; - using parent_structure::parent; - using parent_structure::init; - using distance_structure::distance; - using distance_structure::init; - using distance_structure::INFINITY; - - pathSearch_BFS() - { - } - - void reset() - {} - - template - bool solve ( - const graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - condition_t valid_arc) - // each call resets the state of the tree and distances - // O(|E|+|V|) - { - bool found = false; - - if(!g.is_valid(Source)) - throw std::runtime_error( - "pathSearch_BFS::solve source node is not valid"); - - if(!g.is_valid(Dest)) - throw std::runtime_error( - "pathSearch_BFS::solve destination node is not valid"); - - parent_structure::init(g); - distance_structure::init(g); - - distance.at(Source) = 0; - - std::queue q; - q.push(Source); - - while(!q.empty()) - { - auto node = q.front(); - q.pop(); - - if(node==Dest) - { - found = true; - break; - } - for(auto e: g.out_arcs(node)) - if( valid_arc(e) ) - { - auto [a,b] = g.arc_ends(e); - - if(distance.at(b)==INFINITY) - { - distance.at(b) = distance.at(a)+1; - parent.at(b) = e; - q.push(b); - } - } - } - return found; - } - }; - - class pathSearch_labeling : public parent_structure, - public distance_structure - /* - Represents: shortest path with labeling - Invariant: - - User interface: - Complexity: - */ - { - public: - using value_type = distance_structure::value_type; - using parent_structure::parent; - using parent_structure::init; - using distance_structure::distance; - using distance_structure::init; - using distance_structure::INFINITY; - - node_pos_t last_source{NONE},last_dest{NONE}; - std::vector dist_freq; - - template - void initialize ( - const graph_t &g, - condition_t valid_arc) - { - parent_structure::init(g); - distance_structure::init(g); - - dist_freq.resize(g.num_nodes()+1); - std::fill(dist_freq.begin(),dist_freq.end(),0); - - std::queue q; - distance.at(last_dest)=0; - - // TODO: write a general purpose BFS label solver - q.push(last_dest); - - while(!q.empty()) - { - auto n = q.front(); - q.pop(); - - for(auto e: g.in_arcs(n)) - if( valid_arc(e) ) - { - auto [a,b] = g.arc_ends(e); - value_type dnew = distance[b] + 1; - - if(distance[a]==INFINITY) - { - distance[a] = dnew; - dist_freq.at(dnew)++; - q.push(a); - } - } - } - } - - - public: - - pathSearch_labeling() - {} - - - void reset() - { - last_source = last_dest = node_pos_t{NONE}; - } - template - bool solve ( - const graph_t& g, - const node_pos_t Source, const node_pos_t Dest, - condition_t valid_arc) - { - if(last_source!=Source || last_dest!=Dest) - { - last_source = Source; - last_dest = Dest; - initialize(g,valid_arc); - } - - parent_structure::init(g); - - - for(auto current = Source; - distance.at(Source) - class shortestPath_FIFO : public parent_structure, public distance_structure - /* - Represents: shortest path label-correcting FIFO - Invariant: - - User interface: - Complexity: pseudo-polynomial - */ - { - public: - using value_type = T; - using parent_structure::parent; - using parent_structure::init; - using distance_structure::distance; - using distance_structure::init; - using distance_structure::INFINITY; - - shortestPath_FIFO() - {} - - template - void solve( - const graph_t& g, - const node_pos_t Source, - const std::vector& weight, - condition_t valid_arc) - // shortest path FIFO - // each call resets the state of the tree and distances - // O( pseudo-polynomial ) - { - parent_structure::init(g); - distance_structure::init(g); - - if(!g.is_valid(Source)) - throw std::runtime_error( - "shortestPath_FIFO::solve source node is not valid"); - - if(weight.size() q; - q.push(Source); - distance.at(Source)=0; - - while(!q.empty()) - { - auto node = q.front(); - q.pop(); - - for(auto e: g.out_arcs(node)) - if( valid_arc(e) ) - { - auto [a,b] = g.arc_ends(e); - const value_type dnew = distance.at(a)+weight.at(e); - - if(distance.at(b)>dnew) - { - distance.at(b) = dnew; - parent.at(b) = e; - q.push(b); - } - } - } - } - }; - - template - class shortestPath_BellmanFord : public parent_structure, public distance_structure - /* - Represents: shortest path using Bellman-Ford - Invariant: - - User interface: - Complexity: |V| |E| - */ - { - public: - using value_type = T; - using parent_structure::parent; - using parent_structure::init; - using distance_structure::distance; - using distance_structure::init; - using distance_structure::INFINITY; - - shortestPath_BellmanFord() - {} - - template - void solve ( - const graph_t& g, - const node_pos_t Source, - const std::vector& weight, - condition_t valid_arc) - // shortest path Bellman-Ford - // each call resets the state of the tree and distances - // O(|V||E|) - { - parent_structure::init(g); - distance_structure::init(g); - - if(!g.is_valid(Source)) - throw std::runtime_error( - "shortestPath_BellmanFord::solve source node is not valid"); - - if(weight.size()dnew) - { - distance[b]=dnew; - parent.at(b) = e; - updates = true; - } - } - } - if(! updates) - break; - } - // TODO: check for negative cycles - } - }; - - template - class shortestPath_Dijkstra : public parent_structure, public distance_structure - /* - Represents: shortest path with weights using Dijkstra - Invariant: - - User interface: - Complexity: |E| + |V| log |V| - */ - { - public: - using value_type = T; - using parent_structure::parent; - using parent_structure::init; - using distance_structure::distance; - using distance_structure::init; - using distance_structure::INFINITY; - - shortestPath_Dijkstra() - {} - - template - void solve ( - const graph_t& g, - const node_pos_t Source, - const std::vector& weight, - condition_t valid_arc) - // Dijkstra algorithm - // precondition: doesnt work with negative weights! - // O( |E|+|V| log |V| ) - { - parent_structure::init(g); - distance_structure::init(g); - - if(!g.is_valid(Source)) - throw std::runtime_error( - "shortestPath_Dijkstra::solve source node is not valid"); - - if(weight.size() visited(g.max_num_nodes(),false); - - distance.at(Source) = 0; - std::priority_queue< - std::pair, - std::vector< std::pair >, - std::greater > - > q; - q.push( {0,Source} ); - - while(!q.empty()) - { - const auto [dist,node] = q.top(); - q.pop(); - - if(visited.at(node)) - continue; - - visited[node]=true; - - for(auto e: g.out_arcs(node)) - if( valid_arc(e) ) - { - auto [a,b] = g.arc_ends(e); - - if(weight.at(e)<0) - throw std::runtime_error( - "shortestPath_Dijkstra::solve found a negative edge"); - - value_type dnew = dist + weight.at(e); - if(distance.at(b)>dnew) - { - distance[b] = dnew; - parent[b] = e; - q.push({dnew,b}); - } - } - } - } - }; -} diff --git a/subprojects/MinCostFlow/include/mincostflow/vectorized_map.hpp b/subprojects/MinCostFlow/include/mincostflow/vectorized_map.hpp deleted file mode 100644 index 6ec9071..0000000 --- a/subprojects/MinCostFlow/include/mincostflow/vectorized_map.hpp +++ /dev/null @@ -1,259 +0,0 @@ -#pragma once - -#include -#include -#include -#include -#include -#include -#include - -namespace ln -{ - template - class vectorized_map - /* - a map-like data structure - that holds elements in a vector, hence access is O(1) - new elements are added to the smallest available slot and the key is returned, - we don't choose the key associated to an element, but once the key is set it will - remain valid until the element is removed. - iterators point to key values, not data, - though data can be accessed through the key - */ - { - using size_t = std::size_t; - - static constexpr index_t MAX_IDX = std::numeric_limits::max(); - std::vector< bool > valid_flag; - std::vector< data_t > data; - std::set< index_t > free_slots; - - void check_invariants()const - { - assert(valid_flag.size()==data.size()); - } - - void free_space() - { - while(!data.empty() && !valid_flag.back()) - // eliminate unused elements from the back of the buffer - { - auto ptr = index_t{data.size()-1}; - assert(free_slots.find(ptr)!=free_slots.end()); - - free_slots.erase(ptr); - valid_flag.pop_back(); - data.pop_back(); - } - check_invariants(); - } - - class base_iterator - { - protected: - const vectorized_map& const_ref; - index_t pos; - - void check_invariants()const - { - bool is_valid = const_ref.is_valid(pos); - bool is_infty = size_t(pos)==MAX_IDX; - - assert(is_valid ^ is_infty); - } - - - void next_valid() - { - while(size_t(pos)=const_ref.capacity()) - pos = index_t{MAX_IDX}; - } - - public: - base_iterator(const vectorized_map& c, index_t x): - const_ref{c}, - pos{x} - {} - - }; - - class iterator : public base_iterator - { - using base_iterator::pos; - using base_iterator::next_valid; - using base_iterator::check_invariants; - using base_iterator::const_ref; - - public: - iterator(vectorized_map& c, index_t x): - base_iterator{c} - { - next_valid(); - check_invariants(); - } - - bool operator==(const iterator& that)const - { - return pos == that.pos; - } - bool operator != (const iterator& that)const - { - return !(pos==that.pos); - } - index_t operator * ()const - { - return pos; - } - - iterator& operator++() - { - if(!const_ref.is_valid(pos)) - return *this; - ++pos; - next_valid(); - check_invariants(); - return *this; - } - }; - class const_iterator : public base_iterator - { - using base_iterator::pos; - using base_iterator::next_valid; - using base_iterator::check_invariants; - using base_iterator::const_ref; - - public: - - const_iterator(const vectorized_map& c, index_t x): - base_iterator{c,x} - { - next_valid(); - check_invariants(); - } - - bool operator==(const const_iterator& that)const - { - return pos == that.pos; - } - bool operator != (const const_iterator& that)const - { - return !(pos==that.pos); - } - index_t operator * ()const - { - return pos; - } - - const_iterator& operator++() - { - if(!const_ref.is_valid(pos)) - return *this; - ++pos; - next_valid(); - check_invariants(); - return *this; - } - }; - - public: - bool is_valid(index_t pos)const - { - return size_t(pos) -#include - -#define CHECK(cond,mesg) \ - if(!(cond)) throw std::runtime_error(mesg); - -template -void test_case(const std::vector>& arcs, - const std::vector & capacity, - const int source, - const int sink, - const std::vector& sol) -{ - ln::digraph graph; - pathsolver_t solver; - - graph.add_node(source); - graph.add_node(sink); - - std::vector res_cap; - for(auto i=0UL;i -void test() -{ - { // case 1 - int source=0; - int sink = 1; - std::vector flow{1,0,0,0,0,0}; - - std::vector< std::pair > arc{{0,1},{0,2},{1,3},{1,2},{1,0},{3,1}}; - std::vector capacity{1,9,5,1,7,4}; - - test_case(arc,capacity,source,sink,flow); - } - { // case 2 - int source=0; - int sink = 1; - std::vector flow{1,2,0,1,2}; - - std::vector< std::pair > arc{{0,2},{0,3}, {3,2},{2,1},{3,1}}; - std::vector capacity{1,2,2,2,2}; - - test_case(arc,capacity,source,sink,flow); - } -} - -int main() -{ - using value_type = int; - - try - { - test< ln::maxflow_augmenting_path >(); - test< ln::maxflow_augmenting_path >(); - test< ln::maxflow_preflow >(); - test< ln::maxflow_scaling >(); - test< ln::maxflow_scaling >(); - }catch(...) - { - return 1; - } - return 0; -} - - diff --git a/subprojects/MinCostFlow/test/meson.build b/subprojects/MinCostFlow/test/meson.build deleted file mode 100644 index 89d1e07..0000000 --- a/subprojects/MinCostFlow/test/meson.build +++ /dev/null @@ -1,15 +0,0 @@ -test('Vectorized Map', - executable('vec-map','vec-map.cpp', - dependencies: dep_mincostflow)) - -test('Shortest Path', - executable('short-path','short-path.cpp', - dependencies: dep_mincostflow)) - -test('Max Flow', - executable('max-flow','max-flow.cpp', - dependencies: dep_mincostflow)) - -test('Min Cost Max Flow', - executable('min-cost-flow','min-cost-flow.cpp', - dependencies: dep_mincostflow)) diff --git a/subprojects/MinCostFlow/test/min-cost-flow.cpp b/subprojects/MinCostFlow/test/min-cost-flow.cpp deleted file mode 100644 index bd8f3fb..0000000 --- a/subprojects/MinCostFlow/test/min-cost-flow.cpp +++ /dev/null @@ -1,274 +0,0 @@ -#include "mincostflow/graph.hpp" -#include "mincostflow/maxflow.hpp" -#include "mincostflow/mincostflow.hpp" -#include -#include - -#define CHECK(cond,mesg) \ - if(!(cond)) throw std::runtime_error(mesg); - -template -void test_case(const std::vector>& arcs, - const std::vector & capacity, - const std::vector & weight, - const int source, - const int sink, - const std::vector& sol) -{ - ln::digraph graph; - pathsolver_t solver; - - graph.add_node(source); - graph.add_node(sink); - - std::vector res_cap; - std::vector res_cost; - for(auto i=0UL;i -void test() -{ - { // case 1 - int source=0; - int sink = 1; - std::vector flow{1,0,0,0,0,0}; - - std::vector< std::pair > arc{{0,1},{0,2},{1,3},{1,2},{1,0},{3,1}}; - std::vector capacity{1,9,5,1,7,4}; - std::vector weight{1,1,1,1,1,1}; - - test_case(arc,capacity,weight,source,sink,flow); - } - { // case 2 - int source=0; - int sink = 1; - std::vector flow{1,2,0,1,2}; - - std::vector< std::pair > arc{{0,2},{0,3}, {3,2},{2,1},{3,1}}; - std::vector capacity{1,2,2,2,2}; - std::vector weight{1,1,1,1,1,1}; - - test_case(arc,capacity,weight,source,sink,flow); - } - { - // case 3 - int source=0; - int sink=1; - std::vector flow{2,5,2,0,0}; - std::vector< std::pair > arc{{0,2},{0,1},{2,1},{1,3},{0,3}}; - std::vector capacity{2,5,7,8,6}; - std::vector weight{1,3,2,2,6}; - - test_case(arc,capacity,weight,source,sink,flow); - } - { - // case 4 - int source = 0; - int sink = 1; - std::vector flow{0,4,1,0,0,1,1,0}; - std::vector< std::pair > arc{{0,2},{0,1},{0,3},{1,3},{2,3},{2,1},{3,2},{3,0}}; - std::vector capacity{2,4,3,3,3,1,1,4}; - std::vector weight{2,3,1,0,2,0,0,4}; - - test_case(arc,capacity,weight,source,sink,flow); - } - { - // case 5 - int source = 0; - int sink = 1; - std::vector flow{1,1,0,0,1,2}; - std::vector< std::pair > arc{{0,3},{0,2},{1,2},{1,0},{2,3},{3,1}}; - std::vector capacity{2,1,1,1,4,2}; - std::vector weight{4,1,0,1,2,0}; - - test_case(arc,capacity,weight,source,sink,flow); - - } -} - -int main() -{ - using value_type = int; - - try - { - test< ln::mincostflow_EdmondsKarp< - value_type, - ln::shortestPath_FIFO > >(); - test< ln::mincostflow_EdmondsKarp< - value_type, - ln::shortestPath_BellmanFord > >(); - - // expected failure - // test< ln::mincostflow_EdmondsKarp< - // value_type, - // ln::shortestPath_Dijkstra > >(); - - - test< ln::mincostflow_PrimalDual< - ln::shortestPath_FIFO, - ln::maxflow_augmenting_path >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_FIFO, - ln::maxflow_augmenting_path >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_FIFO, - ln::maxflow_scaling >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_FIFO, - ln::maxflow_scaling >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_FIFO, - ln::maxflow_preflow >>(); - - - test< ln::mincostflow_PrimalDual< - ln::shortestPath_BellmanFord, - ln::maxflow_augmenting_path >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_BellmanFord, - ln::maxflow_augmenting_path >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_BellmanFord, - ln::maxflow_scaling >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_BellmanFord, - ln::maxflow_scaling >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_BellmanFord, - ln::maxflow_preflow >>(); - - - test< ln::mincostflow_PrimalDual< - ln::shortestPath_Dijkstra, - ln::maxflow_augmenting_path >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_Dijkstra, - ln::maxflow_augmenting_path >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_Dijkstra, - ln::maxflow_scaling >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_Dijkstra, - ln::maxflow_scaling >>(); - test< ln::mincostflow_PrimalDual< - ln::shortestPath_Dijkstra, - ln::maxflow_preflow >>(); - - - - - - test< ln::mincostflow_capacityScaling< - ln::shortestPath_FIFO, - ln::maxflow_augmenting_path - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_FIFO, - ln::maxflow_augmenting_path - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_FIFO, - ln::maxflow_scaling - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_FIFO, - ln::maxflow_scaling - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_FIFO, - ln::maxflow_preflow - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_BellmanFord, - ln::maxflow_augmenting_path - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_BellmanFord, - ln::maxflow_augmenting_path - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_BellmanFord, - ln::maxflow_scaling - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_BellmanFord, - ln::maxflow_scaling - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_BellmanFord, - ln::maxflow_preflow - >> (); - - test< ln::mincostflow_capacityScaling< - ln::shortestPath_Dijkstra, - ln::maxflow_augmenting_path - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_Dijkstra, - ln::maxflow_augmenting_path - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_Dijkstra, - ln::maxflow_scaling - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_Dijkstra, - ln::maxflow_scaling - >> (); - test< ln::mincostflow_capacityScaling< - ln::shortestPath_Dijkstra, - ln::maxflow_preflow - >> (); - - test< ln::mincostflow_costScaling< - ln::maxflow_augmenting_path - >> (); - test< ln::mincostflow_costScaling< - ln::maxflow_augmenting_path - >> (); - test< ln::mincostflow_costScaling< - ln::maxflow_scaling - >> (); - test< ln::mincostflow_costScaling< - ln::maxflow_scaling - >> (); - test< ln::mincostflow_costScaling< - ln::maxflow_preflow - >> (); - }catch(std::exception& e) - { - std::cerr << e.what() << '\n'; - return 1; - } - return 0; -} - - diff --git a/subprojects/MinCostFlow/test/short-path.cpp b/subprojects/MinCostFlow/test/short-path.cpp deleted file mode 100644 index 4650435..0000000 --- a/subprojects/MinCostFlow/test/short-path.cpp +++ /dev/null @@ -1,86 +0,0 @@ -#include "mincostflow/graph.hpp" -#include "mincostflow/shortestpath.hpp" -#include -#include - -#define CHECK(cond,mesg) \ - if(!(cond)) throw std::runtime_error(mesg); - -template -void test_case(const std::vector>& arcs, - const std::vector & length, - const int source, - const std::vector& sol) -{ - ln::digraph graph; - pathsolver_t solver; - - for(auto i=0UL;i weights; - for(auto i=0UL;i -void test() -{ - { // case 1 - int source=0; - std::vector distance{0,1,2,6}; - - std::vector< std::pair > arc{{0,1},{0,2},{1,3},{1,2},{1,0},{3,1}}; - std::vector length{1,9,5,1,7,4}; - - test_case(arc,length,source,distance); - } - { // case 2 - int source=0; - std::vector distance{0,4,11,9}; - - std::vector< std::pair > arc{{0,1},{1,3},{1,0},{1,2},{2,1},{3,2}}; - std::vector length{4,5,4,7,7,3}; - - test_case(arc,length,source,distance); - } -} - -int main() -{ - using length_type = int; - - try - { - test< ln::shortestPath_Dijkstra >(); - test< ln::shortestPath_FIFO >(); - test< ln::shortestPath_BellmanFord >(); - }catch(...) - { - return 1; - } - return 0; -} - diff --git a/subprojects/MinCostFlow/test/vec-map.cpp b/subprojects/MinCostFlow/test/vec-map.cpp deleted file mode 100644 index b50e920..0000000 --- a/subprojects/MinCostFlow/test/vec-map.cpp +++ /dev/null @@ -1,41 +0,0 @@ -#include - -#define CHECK(cond) \ - if(!(cond)) return 1; - - -int main() -{ - ln::vectorized_map vm; - - // initial state - CHECK(vm.capacity()==0 && vm.size()==0); - - // add three elements - vm.insert(1); - vm.insert(2); - vm.insert(3); - CHECK(vm.capacity()==3 && vm.size()==3); - - // remove the first - vm.erase(0); - CHECK(vm.capacity()==3 && vm.size()==2); - - // remove twice, no effect - vm.erase(0); - CHECK(vm.capacity()==3 && vm.size()==2); - - // add a new one, to the first empty slot - vm.insert(11); - CHECK(vm.capacity()==3 && vm.size()==3); - - // remove a non-existent index, no effect - vm.erase(4); - CHECK(vm.capacity()==3 && vm.size()==3); - - // remove from the tail, effect capacity changes - vm.erase(1); - vm.erase(2); - CHECK(vm.capacity()==1 && vm.size()==1); - return 0; -} From 8cc8d2a4450bd8a7f40882d586fd262287d5b7cc Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Mon, 15 Aug 2022 11:01:19 +0200 Subject: [PATCH 07/20] added mincostflow submodule --- .gitmodules | 3 +++ subprojects/MinCostFlow | 1 + 2 files changed, 4 insertions(+) create mode 100644 .gitmodules create mode 160000 subprojects/MinCostFlow diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..a7e90b7 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "subprojects/MinCostFlow"] + path = subprojects/MinCostFlow + url = https://github.com/Lagrang3/mincostflow.git diff --git a/subprojects/MinCostFlow b/subprojects/MinCostFlow new file mode 160000 index 0000000..076959c --- /dev/null +++ b/subprojects/MinCostFlow @@ -0,0 +1 @@ +Subproject commit 076959cc37c7d5bf25f7d72de34502bad050b833 From d8e809055d0f0aa9aad24f360df22756bf3ccf88 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Mon, 15 Aug 2022 11:55:35 +0200 Subject: [PATCH 08/20] added meson to main project --- meson.build | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 meson.build diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..9d5d994 --- /dev/null +++ b/meson.build @@ -0,0 +1,11 @@ +project('pickhardtpayments','cpp','c', + default_options : ['cpp_std=c++17', + 'warning_level=3', + 'optimization=3'], + version: '0.0.0') + + +libMCF_proj = subproject('MinCostFlow') +libMCF_dep = libMCF_proj.get_variable('dep_mincostflow') +libMCF = libMCF_proj.get_variable('libmincostflow') + From fa00c92317eabe618f595d96a91eefd210ba0c34 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Mon, 15 Aug 2022 11:56:33 +0200 Subject: [PATCH 09/20] update MinCostFlow submodule --- subprojects/MinCostFlow | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subprojects/MinCostFlow b/subprojects/MinCostFlow index 076959c..c6519b6 160000 --- a/subprojects/MinCostFlow +++ b/subprojects/MinCostFlow @@ -1 +1 @@ -Subproject commit 076959cc37c7d5bf25f7d72de34502bad050b833 +Subproject commit c6519b6ec0e5b6f5150a6f284bd6739321efac2e From e615272dd482929363b6318e99058ac39ab8e7f4 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 18 Aug 2022 08:37:00 +0200 Subject: [PATCH 10/20] update the default MinCostFlow --- subprojects/MinCostFlow | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subprojects/MinCostFlow b/subprojects/MinCostFlow index c6519b6..1d476a0 160000 --- a/subprojects/MinCostFlow +++ b/subprojects/MinCostFlow @@ -1 +1 @@ -Subproject commit c6519b6ec0e5b6f5150a6f284bd6739321efac2e +Subproject commit 1d476a00c296df8d594358f2ff416e6e1c8a6e67 From 349172619ab98406c4f274b042874947de6329a9 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 18 Aug 2022 13:08:11 +0200 Subject: [PATCH 11/20] only channels that changed their cost/capacity are updated this commit is makes sense only after this PR https://github.com/Lagrang3/mincostflow/pull/14 is merged. --- .../SyncSimulatedPaymentSession.py | 141 +++++++++++------- 1 file changed, 86 insertions(+), 55 deletions(-) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index 6e836d7..edd065c 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -1,5 +1,6 @@ from .UncertaintyNetwork import UncertaintyNetwork from .OracleLightningNetwork import OracleLightningNetwork +from .UncertaintyChannel import DEFAULT_N from MinCostFlow import MCFNetwork @@ -31,23 +32,76 @@ def __init__(self, self._oracle = oracle self._uncertainty_network = uncertainty_network self._prune_network = prune_network - self._prepare_integer_indices_for_nodes() - def _prepare_integer_indices_for_nodes(self): - """ - necessary for the OR-lib by google and the min cost flow solver - - - let's initialize the look up tables for node_ids to integers from [0,...,#number of nodes] - this is necessary because of the API of the Google Operations Research min cost flow solver - """ - self._mcf_id = {} - self._node_key = {} - for k, node_id in enumerate(self._uncertainty_network.network.nodes()): - self._mcf_id[node_id] = k - self._node_key[k] = node_id + def _mcf_demands(self,src,dest,amt: int = 1): + # add amount to sending node + self._min_cost_flow.SetNodeSupply( + src, int(amt)) # /QUANTIZATION)) - def _prepare_mcf_solver(self, src, dest, amt: int = 1, mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): + # add -amount to recipient nods + self._min_cost_flow.SetNodeSupply( + dest, -int(amt)) # /QUANTIZATION)) + + def _mcf_channel_encode_part(self,channel,part: int=0): + direction = 0 + if channel.src>channel.dest: + direction = 1 + return direction + part*2 + + def _mcf_good_channel(self,channel,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): + # ignore channels with too large base fee + if channel.base_fee > base_fee: + return False + # FIXME: Remove Magic Number for pruning + # Prune channels away thay have too low success probability! This is a huge runtime boost + # However the pruning would be much better to work on quantiles of normalized cost + # So as soon as we have better Scaling, Centralization and feature engineering we can + # probably have a more focused pruning + if self._prune_network and channel.success_probability(250_000) < 0.9: + return False + return True + + def _mcf_delete_channel(self,channel): + for part in range(DEFAULT_N): + index = self._min_cost_flow.RemoveArc(channel.short_channel_id, + self._mcf_channel_encode_part(channel,part)) + self._arc_to_channel.pop(index) + + def _mcf_update_channel(self,channel,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): + if not self._mcf_good_channel(channel,mu,base_fee): + self._mcf_delete_channel(channel) + return + # QUANTIZATION): + for part,(capacity, cost) in enumerate(channel.get_piecewise_linearized_costs(mu=mu)): + self._min_cost_flow.UpdateArc(channel.short_channel_id, + self._mcf_channel_encode_part(channel,part), + capacity, + cost) + + def _mcf_update_used_channels(self,payments,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): + for attempt in payments.values(): + path = attempt['path'] + for channel in path: + self._mcf_update_channel(channel,mu,base_fee) + + def _mcf_new_arc(self,channel,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): + if not self._mcf_good_channel(channel,mu,base_fee): + return + cnt = 0 + # QUANTIZATION): + for part,(capacity, cost) in enumerate(channel.get_piecewise_linearized_costs(mu=mu)): + index = self._min_cost_flow.AddArc(channel.src, + channel.dest, + channel.short_channel_id, + self._mcf_channel_encode_part(channel,part), + capacity, + cost) + self._arc_to_channel[index] = (channel.src, channel.dest, channel, 0) + if self._prune_network and cnt > 1: + break + cnt += 1 + + def _prepare_mcf_solver(self,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): """ computes the uncertainty network given our prior belief and prepares the min cost flow solver @@ -60,45 +114,12 @@ def _prepare_mcf_solver(self, src, dest, amt: int = 1, mu: int = 100_000_000, ba self._arc_to_channel = {} for s, d, channel in self._uncertainty_network.network.edges(data="channel"): - # ignore channels with too large base fee - if channel.base_fee > base_fee: - continue - # FIXME: Remove Magic Number for pruning - # Prune channels away thay have too low success probability! This is a huge runtime boost - # However the pruning would be much better to work on quantiles of normalized cost - # So as soon as we have better Scaling, Centralization and feature engineering we can - # probably have a more focused pruning - if self._prune_network and channel.success_probability(250_000) < 0.9: - continue - cnt = 0 - # QUANTIZATION): - direction = 0 - if s>d: - direction = 1 - for part,(capacity, cost) in enumerate(channel.get_piecewise_linearized_costs(mu=mu)): - index = self._min_cost_flow.AddArc(s, - d, - channel.short_channel_id, - direction + part*2, - capacity, - cost) - self._arc_to_channel[index] = (s, d, channel, 0) - if self._prune_network and cnt > 1: - break - cnt += 1 + self._mcf_new_arc(channel,mu,base_fee) # Add node supply to 0 for all nodes for i in self._uncertainty_network.network.nodes(): self._min_cost_flow.SetNodeSupply(i, 0) - # add amount to sending node - self._min_cost_flow.SetNodeSupply( - src, int(amt)) # /QUANTIZATION)) - - # add -amount to recipient nods - self._min_cost_flow.SetNodeSupply( - dest, -int(amt)) # /QUANTIZATION)) - def _next_hop(self, path): """ generator to iterate through edges indext by node id of paths @@ -148,11 +169,15 @@ def _disect_flow_to_paths(self, s, d): # first collect all linearized edges which are assigned a non zero flow put them into a networkx graph G = nx.MultiDiGraph() - for i in range(self._min_cost_flow.NumArcs()): - flow = self._min_cost_flow.Flow(i) # *QUANTIZATION - if flow == 0: - continue - + + #for i in range(self._min_cost_flow.NumArcs()): + + index_list, flow_list = self._min_cost_flow.FlowArray() + for i,flow in zip(index_list,flow_list): + #flow = self._min_cost_flow.Flow(i) # *QUANTIZATION + #if flow == 0: + # continue + # print(i,flow) src, dest, channel, _ = self._arc_to_channel[i] if G.has_edge(src, dest): if channel.short_channel_id in G[src][dest]: @@ -198,7 +223,9 @@ def _generate_candidate_paths(self, src, dest, amt, mu: int = 100_000_000, base: """ # First we prepare the min cost flow by getting arcs from the uncertainty network - self._prepare_mcf_solver(src, dest, amt, mu, base) + # self._prepare_mcf_solver(mu, base) + + self._mcf_demands(src,dest,amt) start = time.time() #print("solving mcf...") @@ -347,6 +374,8 @@ def pickhardt_pay(self, src, dest, amt, mu=1, base=0): number_number_of_onions = 0 total_number_failed_paths = 0 + self._prepare_mcf_solver(mu,base) + # This is the main payment loop. It is currently blocking and synchronous but may be # implemented in a concurrent way. Also we stop after 10 rounds which is pretty arbitrary # a better stop criteria would be if we compute infeasable flows or if the probabilities @@ -364,6 +393,8 @@ def pickhardt_pay(self, src, dest, amt, mu=1, base=0): # matke attempts and update our information about the UncertaintyNetwork self._attempt_payments(payments) + + self._mcf_update_used_channels(payments,mu,base) # run some simple statistics and depict them amt, paid_fees, num_paths, number_failed_paths = self._evaluate_attempts( From 656c054931d926d0005d6a5d89211cef2a33750b Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 18 Aug 2022 20:54:20 +0200 Subject: [PATCH 12/20] testing the simpleMCFNetwork solver --- pickhardtpayments/SyncSimulatedPaymentSession.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index edd065c..dc7fb47 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -2,7 +2,7 @@ from .OracleLightningNetwork import OracleLightningNetwork from .UncertaintyChannel import DEFAULT_N -from MinCostFlow import MCFNetwork +from MinCostFlow import MCFNetwork, simpleMCFNetwork from typing import List @@ -110,7 +110,9 @@ def _prepare_mcf_solver(self,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE returns the instantiated min_cost_flow object from the google OR-lib that contains the piecewise linearized problem """ - self._min_cost_flow = MCFNetwork() + start = time.time() + # self._min_cost_flow = MCFNetwork() + self._min_cost_flow = simpleMCFNetwork() self._arc_to_channel = {} for s, d, channel in self._uncertainty_network.network.edges(data="channel"): @@ -119,6 +121,8 @@ def _prepare_mcf_solver(self,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE # Add node supply to 0 for all nodes for i in self._uncertainty_network.network.nodes(): self._min_cost_flow.SetNodeSupply(i, 0) + end = time.time() + return end-start def _next_hop(self, path): """ @@ -374,7 +378,8 @@ def pickhardt_pay(self, src, dest, amt, mu=1, base=0): number_number_of_onions = 0 total_number_failed_paths = 0 - self._prepare_mcf_solver(mu,base) + total_mcf_time = 0 + time_prepare_mcf = self._prepare_mcf_solver(mu,base) # This is the main payment loop. It is currently blocking and synchronous but may be # implemented in a concurrent way. Also we stop after 10 rounds which is pretty arbitrary @@ -387,7 +392,7 @@ def pickhardt_pay(self, src, dest, amt, mu=1, base=0): # transfer to a min cost flow problem and rund the solver paths, runtime = self._generate_candidate_paths( src, dest, amt, mu, base) - + total_mcf_time += runtime # compute some statistics about candidate paths payments = self._estimate_payment_statistics(paths) @@ -415,6 +420,8 @@ def pickhardt_pay(self, src, dest, amt, mu=1, base=0): print("Number of failed onions: ", total_number_failed_paths) print("Failure rate: {:4.2f}% ".format( total_number_failed_paths*100./number_number_of_onions)) + print("runtime for graph initialization: {:4.3f} sec".format(time_prepare_mcf)) + print("total runtime for MCF solve: {:4.3f} sec".format(total_mcf_time)) print("total runtime (including inefficient memory managment): {:4.3f} sec".format( end-start)) print("Learnt entropy: {:5.2f} bits".format(entropy_start-entropy_end)) From 822d1aa635809d0521a1031173d657dea9c4d4a1 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 18 Aug 2022 21:58:29 +0200 Subject: [PATCH 13/20] rollback to the MCFNetwork --- pickhardtpayments/SyncSimulatedPaymentSession.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index dc7fb47..278cd5b 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -111,8 +111,8 @@ def _prepare_mcf_solver(self,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE returns the instantiated min_cost_flow object from the google OR-lib that contains the piecewise linearized problem """ start = time.time() - # self._min_cost_flow = MCFNetwork() - self._min_cost_flow = simpleMCFNetwork() + self._min_cost_flow = MCFNetwork() + # self._min_cost_flow = simpleMCFNetwork() self._arc_to_channel = {} for s, d, channel in self._uncertainty_network.network.edges(data="channel"): From 1f135f32258d727e2f2dc3caa724137d0184294b Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Fri, 19 Aug 2022 10:23:17 +0200 Subject: [PATCH 14/20] new naming convention for MCF solvers --- pickhardtpayments/SyncSimulatedPaymentSession.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index 278cd5b..013edf0 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -2,7 +2,7 @@ from .OracleLightningNetwork import OracleLightningNetwork from .UncertaintyChannel import DEFAULT_N -from MinCostFlow import MCFNetwork, simpleMCFNetwork +from MinCostFlow import MCFAugmentingPaths, MCFCostScaling from typing import List @@ -111,8 +111,8 @@ def _prepare_mcf_solver(self,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE returns the instantiated min_cost_flow object from the google OR-lib that contains the piecewise linearized problem """ start = time.time() - self._min_cost_flow = MCFNetwork() - # self._min_cost_flow = simpleMCFNetwork() + self._min_cost_flow = MCFAugmentingPaths() + # self._min_cost_flow = MCFCostScaling() self._arc_to_channel = {} for s, d, channel in self._uncertainty_network.network.edges(data="channel"): From 49f67e5f04e9ec9f55f96fe4f23dd79e32879f7d Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Fri, 19 Aug 2022 10:26:39 +0200 Subject: [PATCH 15/20] tracking changes in MinCostFlow master branch --- subprojects/MinCostFlow | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subprojects/MinCostFlow b/subprojects/MinCostFlow index 1d476a0..ee7965a 160000 --- a/subprojects/MinCostFlow +++ b/subprojects/MinCostFlow @@ -1 +1 @@ -Subproject commit 1d476a00c296df8d594358f2ff416e6e1c8a6e67 +Subproject commit ee7965af55d3afc01f96694f4d46e633652f4cbb From fb9b70461edb0b6caccd02f42146bdf1ba1ac515 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 25 Aug 2022 09:44:13 +0200 Subject: [PATCH 16/20] uploading the files I used for benchmark --- benchmark/analize.py | 73 +++++++++++ benchmark/benchmark.py | 83 ++++++++++++ .../SyncSimulatedPaymentSession.py | 121 ++++++++++-------- 3 files changed, 222 insertions(+), 55 deletions(-) create mode 100644 benchmark/analize.py create mode 100644 benchmark/benchmark.py diff --git a/benchmark/analize.py b/benchmark/analize.py new file mode 100644 index 0000000..fd413cf --- /dev/null +++ b/benchmark/analize.py @@ -0,0 +1,73 @@ +import json +import matplotlib.pyplot as plt +import numpy + +def open_json_data(filename): + with open(filename,'r') as f: + return json.load(f) + +ortool_data = open_json_data('ortools.json') +apath_data = open_json_data('augmenting_path.json') +cscaling_data = open_json_data('cost_scaling.json') + +def plot_success_rate(data): + x = [ int(i) for i in data.keys() ] + y = [ data[i]['nsuccess']/(data[i]['nsuccess']+data[i]['nfails']) for i in data.keys() ] + fig = plt.figure() + ax = fig.subplots() + ax.grid() + ax.set_title('Success rate') + ax.set_xlabel('amount (sats)') + ax.set_ylabel('success rate') + ax.loglog(x,y) + fig.savefig('plot_success.png') + plt.show() + +def plot_time(data): + fig=plt.figure() + ax=fig.subplots() + ax.grid() + ax.set_title('MCF timing') + ax.set_xlabel('amount (sats)') + ax.set_ylabel('time (ms)') + for k in data: + D = data[k] + xser = [ int(i) for i in D.keys() ] + tser = [ numpy.mean(D[i]['time_mcf'])*1000 for i in D.keys() ] + tmin = [ numpy.min(D[i]['time_mcf'])*1000 for i in D.keys() ] + tmax = [ numpy.max(D[i]['time_mcf'])*1000 for i in D.keys() ] + ax.fill_between(xser,tmin,tmax,alpha = 0.2) + ax.loglog(xser,tser,label=k) + ax.legend() + fig.savefig('plot_time.png') + plt.show() + #x = [ int(i) for i in data['ortools'].keys() ] + +def plot_totaltime(data): + fig=plt.figure() + ax=fig.subplots() + ax.grid() + ax.set_title('Total Payment Time') + ax.set_xlabel('amount (sats)') + ax.set_ylabel('time (ms)') + for k in data: + D = data[k] + xser = [ int(i) for i in D.keys() ] + tser = [ numpy.mean(D[i]['time_total'])*1000 for i in D.keys() ] + tmin = [ numpy.min(D[i]['time_total'])*1000 for i in D.keys() ] + tmax = [ numpy.max(D[i]['time_total'])*1000 for i in D.keys() ] + ax.fill_between(xser,tmin,tmax,alpha = 0.2) + ax.loglog(xser,tser,label=k) + ax.legend() + fig.savefig('plot_totaltime.png') + plt.show() + + +#print(ortool_data) +#print(apath_data) +#print(cscaling_data) + +plot_success_rate(ortool_data) +in_data = {"ortools" : ortool_data, "cost_scaling" : cscaling_data, "augmenting_path" : apath_data} +plot_time(in_data) +plot_totaltime(in_data) diff --git a/benchmark/benchmark.py b/benchmark/benchmark.py new file mode 100644 index 0000000..3284389 --- /dev/null +++ b/benchmark/benchmark.py @@ -0,0 +1,83 @@ +import random +import json + +#prefix='cost_scaling' +prefix='augmenting_path' +#prefix='ortools' +random.seed(64) + +from pickhardtpayments.ChannelGraph import ChannelGraph +from pickhardtpayments.UncertaintyNetwork import UncertaintyNetwork +from pickhardtpayments.OracleLightningNetwork import OracleLightningNetwork +from pickhardtpayments.SyncSimulatedPaymentSession import SyncSimulatedPaymentSession + + +#we first need to import the chanenl graph from c-lightning jsondump +#you can get your own data set via: +# $: lightning-cli listchannels > listchannels20220412.json +# alternatively you can go to https://ln.rene-pickhardt.de to find a data dump +channel_graph = ChannelGraph("listchannels20220412.json") + +uncertainty_network = UncertaintyNetwork(channel_graph) +oracle_lightning_network = OracleLightningNetwork(channel_graph) +#we create ourselves a payment session which in this case operates by sending out the onions +#sequentially +payment_session = SyncSimulatedPaymentSession(oracle_lightning_network, + uncertainty_network, + mu = 0, + base = 0, + prune_network=False) + +#we need to make sure we forget all learnt information on the Uncertainty Nework + +net = uncertainty_network._channel_graph +nodes = net.nodes() +funds = {} +for n in nodes: + funds[n] = 0 + +for a,b,d in net.edges(data=True): + cap = d['channel'].capacity + funds[a] += cap + +# exit(0) + +def random_node_pairs(n_pairs,amt): + list_of_pairs = [] + for i in range(n_pairs): + while(True): + a,b = random.sample(nodes,2) + if min(funds[a],funds[b])>=amt: + if oracle_lightning_network.theoretical_maximum_payable_amount(a,b,base_fee=0)>=amt: + list_of_pairs.append((a,b)) + break + return list_of_pairs + +#we run the simulation of pickhardt payments and track all the results + +n_pairs = 10 +amount_list = [ 2**i for i in range(24,25)] +stat = {} + +for a in amount_list: + stat[a] = { 'time_mcf': [], 'time_total': [], 'nfails': 0 , 'nsuccess': 0 } + +# random_node_pairs(n_pairs) + +with open(prefix+'.log','w') as flog: + for amt in amount_list: + for (A,B) in random_node_pairs(n_pairs,amt): + print("trying a payment of",amt,"sats from",A,"to",B) + print("trying a payment of",amt,"sats from",A,"to",B, file=flog) + payment_session.forget_information() + try: + time_mcf,time_total = payment_session.pickhardt_pay(A,B, amt,log_out=flog) + stat[amt]['time_mcf'].append( time_mcf ) + stat[amt]['time_total'].append( time_total ) + stat[amt]['nsuccess'] += 1 + except: + stat[amt]['nfails'] +=1 + continue + +with open(prefix+'.json','w') as f: + json.dump(stat,f) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index 013edf0..dd5920b 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -2,13 +2,13 @@ from .OracleLightningNetwork import OracleLightningNetwork from .UncertaintyChannel import DEFAULT_N -from MinCostFlow import MCFAugmentingPaths, MCFCostScaling +from MinCostFlow import MCFNetwork from typing import List import time import networkx as nx - +import sys DEFAULT_BASE_THRESHOLD = 0 @@ -28,10 +28,14 @@ class SyncSimulatedPaymentSession(): def __init__(self, oracle: OracleLightningNetwork, uncertainty_network: UncertaintyNetwork, + mu = 1, base = 0, prune_network: bool = True): self._oracle = oracle self._uncertainty_network = uncertainty_network self._prune_network = prune_network + self._mu = mu + self._base_fee_threshold = base + self._prepare_mcf_solver() def _mcf_demands(self,src,dest,amt: int = 1): # add amount to sending node @@ -48,9 +52,9 @@ def _mcf_channel_encode_part(self,channel,part: int=0): direction = 1 return direction + part*2 - def _mcf_good_channel(self,channel,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): + def _mcf_good_channel(self,channel): # ignore channels with too large base fee - if channel.base_fee > base_fee: + if channel.base_fee > self._base_fee_threshold: return False # FIXME: Remove Magic Number for pruning # Prune channels away thay have too low success probability! This is a huge runtime boost @@ -67,29 +71,29 @@ def _mcf_delete_channel(self,channel): self._mcf_channel_encode_part(channel,part)) self._arc_to_channel.pop(index) - def _mcf_update_channel(self,channel,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): - if not self._mcf_good_channel(channel,mu,base_fee): + def _mcf_update_channel(self,channel): + if not self._mcf_good_channel(channel): self._mcf_delete_channel(channel) return # QUANTIZATION): - for part,(capacity, cost) in enumerate(channel.get_piecewise_linearized_costs(mu=mu)): + for part,(capacity, cost) in enumerate(channel.get_piecewise_linearized_costs(mu=self._mu)): self._min_cost_flow.UpdateArc(channel.short_channel_id, self._mcf_channel_encode_part(channel,part), capacity, cost) - def _mcf_update_used_channels(self,payments,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): + def _mcf_update_used_channels(self,payments): for attempt in payments.values(): path = attempt['path'] for channel in path: - self._mcf_update_channel(channel,mu,base_fee) + self._mcf_update_channel(channel) - def _mcf_new_arc(self,channel,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): - if not self._mcf_good_channel(channel,mu,base_fee): + def _mcf_new_arc(self,channel): + if not self._mcf_good_channel(channel): return cnt = 0 # QUANTIZATION): - for part,(capacity, cost) in enumerate(channel.get_piecewise_linearized_costs(mu=mu)): + for part,(capacity, cost) in enumerate(channel.get_piecewise_linearized_costs(mu=self._mu)): index = self._min_cost_flow.AddArc(channel.src, channel.dest, channel.short_channel_id, @@ -101,7 +105,7 @@ def _mcf_new_arc(self,channel,mu: int = 100_000_000, base_fee: int = DEFAULT_BAS break cnt += 1 - def _prepare_mcf_solver(self,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE_THRESHOLD): + def _prepare_mcf_solver(self): """ computes the uncertainty network given our prior belief and prepares the min cost flow solver @@ -111,12 +115,11 @@ def _prepare_mcf_solver(self,mu: int = 100_000_000, base_fee: int = DEFAULT_BASE returns the instantiated min_cost_flow object from the google OR-lib that contains the piecewise linearized problem """ start = time.time() - self._min_cost_flow = MCFAugmentingPaths() - # self._min_cost_flow = MCFCostScaling() + self._min_cost_flow = MCFNetwork() self._arc_to_channel = {} for s, d, channel in self._uncertainty_network.network.edges(data="channel"): - self._mcf_new_arc(channel,mu,base_fee) + self._mcf_new_arc(channel) # Add node supply to 0 for all nodes for i in self._uncertainty_network.network.nodes(): @@ -214,7 +217,7 @@ def _disect_flow_to_paths(self, s, d): G.remove_edge(src, dest, key=channel.short_channel_id) return channel_paths - def _generate_candidate_paths(self, src, dest, amt, mu: int = 100_000_000, base: int = DEFAULT_BASE_THRESHOLD): + def _generate_candidate_paths(self, src, dest, amt): """ computes the optimal payment split to deliver `amt` from `src` to `dest` and updates our belief about the liquidity @@ -234,7 +237,9 @@ def _generate_candidate_paths(self, src, dest, amt, mu: int = 100_000_000, base: start = time.time() #print("solving mcf...") try: - self._min_cost_flow.Solve() + #self._min_cost_flow.Solve() + self._min_cost_flow._Solve_by_AugmentingPaths() + #self._min_cost_flow._Solve_by_CostScaling() except: raise BaseException('There was an issue with the min cost flow input.') @@ -284,7 +289,7 @@ def _attempt_payments(self, payments): self._uncertainty_network.allocate_amount_on_path( attempt["path"], attempt["amount"]) - def _evaluate_attempts(self, payments): + def _evaluate_attempts(self, payments, log_out = sys.stdout): """ helper function to collect statistics about attempts and print them @@ -296,11 +301,11 @@ def _evaluate_attempts(self, payments): number_failed_paths = 0 expected_sats_to_deliver = 0 amt = 0 - print("\nStatistics about {} candidate onions:\n".format(len(payments))) + print("\nStatistics about {} candidate onions:\n".format(len(payments)), file=log_out) has_failed_attempt = False - print("successful attempts:") - print("--------------------") + print("successful attempts:", file=log_out) + print("--------------------", file=log_out) for attempt in payments.values(): success = attempt["success"] if success == False: @@ -314,12 +319,12 @@ def _evaluate_attempts(self, payments): total_fees += fee expected_sats_to_deliver += probability * amount print(" p = {:6.2f}% amt: {:9} sats hops: {} ppm: {:5}".format( - probability*100, amount, len(path), int(fee*1000_000/amount))) + probability*100, amount, len(path), int(fee*1000_000/amount)), file=log_out) paid_fees += fee if has_failed_attempt: - print("\nfailed attempts:") - print("----------------") + print("\nfailed attempts:", file=log_out) + print("----------------", file=log_out) for attempt in payments.values(): success = attempt["success"] if success: @@ -332,30 +337,34 @@ def _evaluate_attempts(self, payments): total_fees += fee expected_sats_to_deliver += probability * amount print(" p = {:6.2f}% amt: {:9} sats hops: {} ppm: {:5} ".format( - probability*100, amount, len(path), int(fee*1000_000/amount))) + probability*100, amount, len(path), int(fee*1000_000/amount)), file=log_out) number_failed_paths += 1 residual_amt += amount - print("\nAttempt Summary:") - print("=================") - print("\nTried to deliver {:10} sats".format(amt)) + print("\nAttempt Summary:", file=log_out) + print("=================", file=log_out) + print("\nTried to deliver {:10} sats".format(amt), file=log_out) fraction = expected_sats_to_deliver*100./amt print("expected to deliver {:10} sats \t({:4.2f}%)".format( - int(expected_sats_to_deliver), fraction)) + int(expected_sats_to_deliver), fraction), file=log_out) fraction = (amt-residual_amt)*100./(amt) print("actually deliverd {:10} sats \t({:4.2f}%)".format( - amt-residual_amt, fraction)) + amt-residual_amt, fraction), file=log_out) print("deviation: {:4.2f}".format( - (amt-residual_amt)/(expected_sats_to_deliver+1))) - print("planned_fee: {:8.3f} sat".format(total_fees)) - print("paid fees: {:8.3f} sat".format(paid_fees)) + (amt-residual_amt)/(expected_sats_to_deliver+1)), file=log_out) + print("planned_fee: {:8.3f} sat".format(total_fees), file=log_out) + print("paid fees: {:8.3f} sat".format(paid_fees), file=log_out) return residual_amt, paid_fees, len(payments), number_failed_paths def forget_information(self): """ forgets all the information in the UncertaintyNetwork that is a member of the PaymentSession """ + start = time.time() self._uncertainty_network.reset_uncertainty_network() + self._min_cost_flow.Forget() + end = time.time() + return end-start def activate_network_wide_uncertainty_reduction(self, n): """ @@ -364,7 +373,7 @@ def activate_network_wide_uncertainty_reduction(self, n): self._uncertainty_network.activate_network_wide_uncertainty_reduction( n, self._oracle) - def pickhardt_pay(self, src, dest, amt, mu=1, base=0): + def pickhardt_pay(self, src, dest, amt, log_out = sys.stdout): """ conduct one experiment! might need to call oracle.reset_uncertainty_network() first I could not put it here as some experiments require sharing of liqudity information @@ -379,19 +388,20 @@ def pickhardt_pay(self, src, dest, amt, mu=1, base=0): total_number_failed_paths = 0 total_mcf_time = 0 - time_prepare_mcf = self._prepare_mcf_solver(mu,base) + # time_prepare_mcf = self._prepare_mcf_solver(mu,base) + time_prepare_mcf = self.forget_information() # This is the main payment loop. It is currently blocking and synchronous but may be # implemented in a concurrent way. Also we stop after 10 rounds which is pretty arbitrary # a better stop criteria would be if we compute infeasable flows or if the probabilities # are to low or residual amounts decrease to slowly while amt > 0 and cnt < 10: - print("Round number: ", cnt+1) - print("Try to deliver", amt, "satoshi:") + print("Round number: ", cnt+1, file=log_out) + print("Try to deliver", amt, "satoshi:", file=log_out) # transfer to a min cost flow problem and rund the solver paths, runtime = self._generate_candidate_paths( - src, dest, amt, mu, base) + src, dest, amt) total_mcf_time += runtime # compute some statistics about candidate paths payments = self._estimate_payment_statistics(paths) @@ -399,13 +409,13 @@ def pickhardt_pay(self, src, dest, amt, mu=1, base=0): # matke attempts and update our information about the UncertaintyNetwork self._attempt_payments(payments) - self._mcf_update_used_channels(payments,mu,base) + self._mcf_update_used_channels(payments) # run some simple statistics and depict them amt, paid_fees, num_paths, number_failed_paths = self._evaluate_attempts( - payments) - print("Runtime of flow computation: {:4.2f} sec ".format(runtime)) - print("\n================================================================\n") + payments, log_out) + print("Runtime of flow computation: {:4.2f} sec ".format(runtime), file=log_out) + print("\n================================================================\n", file=log_out) number_number_of_onions += num_paths total_number_failed_paths += number_failed_paths @@ -413,18 +423,19 @@ def pickhardt_pay(self, src, dest, amt, mu=1, base=0): cnt += 1 end = time.time() entropy_end = self._uncertainty_network.entropy() - print("SUMMARY:") - print("========") - print("Rounds of mcf-computations: ", cnt) - print("Number of onions sent: ", number_number_of_onions) - print("Number of failed onions: ", total_number_failed_paths) + print("SUMMARY:", file=log_out) + print("========", file=log_out) + print("Rounds of mcf-computations: ", cnt, file=log_out) + print("Number of onions sent: ", number_number_of_onions, file=log_out) + print("Number of failed onions: ", total_number_failed_paths, file=log_out) print("Failure rate: {:4.2f}% ".format( - total_number_failed_paths*100./number_number_of_onions)) - print("runtime for graph initialization: {:4.3f} sec".format(time_prepare_mcf)) - print("total runtime for MCF solve: {:4.3f} sec".format(total_mcf_time)) + total_number_failed_paths*100./number_number_of_onions), file=log_out) + print("runtime for graph initialization: {:4.3f} sec".format(time_prepare_mcf), file=log_out) + print("total runtime for MCF solve: {:4.3f} sec".format(total_mcf_time), file=log_out) print("total runtime (including inefficient memory managment): {:4.3f} sec".format( - end-start)) - print("Learnt entropy: {:5.2f} bits".format(entropy_start-entropy_end)) + end-start), file=log_out) + print("Learnt entropy: {:5.2f} bits".format(entropy_start-entropy_end), file=log_out) print("Fees for successfull delivery: {:8.3f} sat --> {} ppm".format( - total_fees, int(total_fees*1000*1000/full_amt))) - print("used mu:", mu) + total_fees, int(total_fees*1000*1000/full_amt)), file=log_out) + print("used mu:", self._mu, file=log_out) + return total_mcf_time,end-start From c5315d2772763a9c0d9dd02dfb83a261c44ccbbc Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 25 Aug 2022 12:44:13 +0200 Subject: [PATCH 17/20] show error code on MCFNetwork failure --- pickhardtpayments/SyncSimulatedPaymentSession.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index dd5920b..2d41678 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -237,11 +237,12 @@ def _generate_candidate_paths(self, src, dest, amt): start = time.time() #print("solving mcf...") try: - #self._min_cost_flow.Solve() - self._min_cost_flow._Solve_by_AugmentingPaths() + self._min_cost_flow.Solve() + #self._min_cost_flow._Solve_by_AugmentingPaths() #self._min_cost_flow._Solve_by_CostScaling() except: - raise BaseException('There was an issue with the min cost flow input.') + raise BaseException('There was an issue with the min cost flow. Error code %d' % + self._min_cost_flow.error_code) paths = self._disect_flow_to_paths(src, dest) end = time.time() From 01ae6b8be6842b82056c21d8e05c270fe57df467 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 25 Aug 2022 12:51:40 +0200 Subject: [PATCH 18/20] updated submodule --- subprojects/MinCostFlow | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subprojects/MinCostFlow b/subprojects/MinCostFlow index ee7965a..6c1925e 160000 --- a/subprojects/MinCostFlow +++ b/subprojects/MinCostFlow @@ -1 +1 @@ -Subproject commit ee7965af55d3afc01f96694f4d46e633652f4cbb +Subproject commit 6c1925e661ccf92bf82b58838d8bb869333ef291 From b1e28d5d0683828a9c52709029f2e7d35c343b9f Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 25 Aug 2022 15:18:51 +0200 Subject: [PATCH 19/20] we forgot to update the channels in the MCF solver --- pickhardtpayments/SyncSimulatedPaymentSession.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pickhardtpayments/SyncSimulatedPaymentSession.py b/pickhardtpayments/SyncSimulatedPaymentSession.py index 443417b..21b8a64 100644 --- a/pickhardtpayments/SyncSimulatedPaymentSession.py +++ b/pickhardtpayments/SyncSimulatedPaymentSession.py @@ -108,10 +108,9 @@ def _mcf_update_channel(self,channel): capacity, cost) - def _mcf_update_used_channels(self,payments): - for attempt in payments.values(): - path = attempt['path'] - for channel in path: + def _mcf_update_used_channels(self, attempts: List[Attempt]): + for attempt in attempts: + for channel in attempt.path: self._mcf_update_channel(channel) def _mcf_new_arc(self,channel): @@ -419,6 +418,7 @@ def pickhardt_pay(self, src, dest, amt, log_out = sys.stdout): # make attempts, try to send onion and register if success or not # update our information about the UncertaintyNetwork self._attempt_payments(sub_payment.attempts) + self._mcf_update_used_channels(sub_payment.attempts) # run some simple statistics and depict them amt, paid_fees, num_paths, number_failed_paths = self._evaluate_attempts( From f94879c0410946b90b059eee660bb79689b43033 Mon Sep 17 00:00:00 2001 From: Lagrang3 Date: Thu, 25 Aug 2022 15:24:25 +0200 Subject: [PATCH 20/20] adapted the basicexample.py --- examples/basicexample.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/basicexample.py b/examples/basicexample.py index dcf6b7b..bf1bf0b 100644 --- a/examples/basicexample.py +++ b/examples/basicexample.py @@ -16,6 +16,7 @@ #sequentially payment_session = SyncSimulatedPaymentSession(oracle_lightning_network, uncertainty_network, + mu=0,base=0, prune_network=False) #we need to make sure we forget all learnt information on the Uncertainty Nework @@ -29,4 +30,4 @@ C_OTTO = "027ce055380348d7812d2ae7745701c9f93e70c1adeb2657f053f91df4f2843c71" tested_amount = 10_000_000 #10 million sats -payment_session.pickhardt_pay(RENE, C_OTTO, tested_amount, mu=0, base=0) +payment_session.pickhardt_pay(RENE, C_OTTO, tested_amount)