diff --git a/Cargo.lock b/Cargo.lock index d1b5e53..2926032 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -13,6 +13,15 @@ dependencies = [ "version_check", ] +[[package]] +name = "aho-corasick" +version = "0.7.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e37cfd5e7657ada45f742d6e99ca5788580b5c529dc78faf11ece6dc702656f" +dependencies = [ + "memchr", +] + [[package]] name = "ansi_term" version = "0.12.1" @@ -28,6 +37,24 @@ version = "1.0.44" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "61604a8f862e1d5c3229fdd78f8b02c68dcf73a4c4b05fd636d12240aaa242c1" +[[package]] +name = "approx" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3f2a05fd1bd10b2527e20a2cd32d8873d115b8b39fe219ee25f42a8aca6ba278" +dependencies = [ + "num-traits", +] + +[[package]] +name = "approx" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "072df7202e63b127ab55acfe16ce97013d5b97bf160489336d3f1840fd78e99e" +dependencies = [ + "num-traits", +] + [[package]] name = "arrayvec" version = "0.5.2" @@ -252,6 +279,12 @@ version = "3.7.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d9df67f7bf9ef8498769f994239c45613ef0c5899415fb58e9add412d2c1a538" +[[package]] +name = "bytemuck" +version = "1.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72957246c41db82b8ef88a5486143830adeb8227ef9837740bdec67724cf2c5b" + [[package]] name = "byteorder" version = "1.4.3" @@ -464,6 +497,16 @@ dependencies = [ "subtle", ] +[[package]] +name = "ctor" +version = "0.1.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ccc0a48a9b826acdf4028595adc9db92caea352f7af011a3034acd172a52a0aa" +dependencies = [ + "quote", + "syn", +] + [[package]] name = "daemon" version = "0.1.0" @@ -478,6 +521,12 @@ dependencies = [ "futures", "hex", "hkdf", + "itertools", + "nalgebra", + "ndarray", + "ndarray_einsum_beta", + "num", + "pretty_assertions", "rand 0.6.5", "reqwest", "rocket", @@ -572,6 +621,12 @@ dependencies = [ "syn", ] +[[package]] +name = "diff" +version = "0.1.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0e25ea47919b1560c4e3b7fe0aaab9becf5b84a10325ddf7db0f0ba5e1026499" + [[package]] name = "digest" version = "0.9.0" @@ -854,9 +909,9 @@ checksum = "9b919933a397b79c37e33b77bb2aa3dc8eb6e165ad809e58ff75bc7db2e34574" [[package]] name = "h2" -version = "0.3.4" +version = "0.3.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d7f3675cfef6a30c8031cf9e6493ebdc3bb3272a3fea3923c4210d1830e6a472" +checksum = "6c06815895acec637cd6ed6e9662c935b866d20a106f8361892893a7d9234964" dependencies = [ "bytes", "fnv", @@ -1087,9 +1142,9 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" [[package]] name = "libc" -version = "0.2.102" +version = "0.2.103" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a2a5ac8f984bfcf3a823267e5fde638acc3325f6496633a5da6bb6eb2171e103" +checksum = "dd8f7255a17a627354f321ef0055d63b898c6fb27eff628af4d1b66b7331edf6" [[package]] name = "libsqlite3-sys" @@ -1160,6 +1215,15 @@ version = "0.1.9" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a3e378b66a060d48947b590737b30a1be76706c8dd7b8ba0f2fe3989c68a853f" +[[package]] +name = "matrixmultiply" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5a8a15b776d9dfaecd44b03c5828c2199cddff5247215858aac14624f8d6b741" +dependencies = [ + "rawpointer", +] + [[package]] name = "memchr" version = "2.4.1" @@ -1174,9 +1238,9 @@ checksum = "2a60c7ce501c71e03a9c9c0d35b861413ae925bd979cc7a4e30d060069aaac8d" [[package]] name = "minimal-lexical" -version = "0.1.3" +version = "0.1.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0c835948974f68e0bd58636fc6c5b1fbff7b297e3046f11b3b3c18bbac012c6d" +checksum = "9c64630dcdd71f1a64c435f54885086a0de5d6a12d104d69b165fb7d5286d677" [[package]] name = "miniscript" @@ -1230,6 +1294,47 @@ dependencies = [ "version_check", ] +[[package]] +name = "nalgebra" +version = "0.29.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d506eb7e08d6329505faa8a3a00a5dcc6de9f76e0c77e4b75763ae3c770831ff" +dependencies = [ + "approx 0.5.0", + "matrixmultiply", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "ndarray" +version = "0.15.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08e854964160a323e65baa19a0b1a027f76d590faba01f05c0cbc3187221a8c9" +dependencies = [ + "approx 0.4.0", + "matrixmultiply", + "num-complex", + "num-integer", + "num-traits", + "rawpointer", +] + +[[package]] +name = "ndarray_einsum_beta" +version = "0.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "668b3abeae3e0637740340e0e32a9bf9308380e146ea6797950f9ff16e88d88a" +dependencies = [ + "lazy_static", + "ndarray", + "num-traits", + "regex", +] + [[package]] name = "nom" version = "7.0.0" @@ -1250,6 +1355,40 @@ dependencies = [ "winapi 0.3.9", ] +[[package]] +name = "num" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "43db66d1170d347f9a065114077f7dccb00c1b9478c89384490a3425279a4606" +dependencies = [ + "num-bigint", + "num-complex", + "num-integer", + "num-iter", + "num-rational", + "num-traits", +] + +[[package]] +name = "num-bigint" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "74e768dff5fb39a41b3bcd30bb25cf989706c90d028d1ad71971987aa309d535" +dependencies = [ + "autocfg 1.0.1", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-complex" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "26873667bbbb7c5182d4a37c1add32cdf09f841af72da53318fdb81543c15085" +dependencies = [ + "num-traits", +] + [[package]] name = "num-integer" version = "0.1.44" @@ -1260,6 +1399,29 @@ dependencies = [ "num-traits", ] +[[package]] +name = "num-iter" +version = "0.1.42" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b2021c8337a54d21aca0d59a92577a029af9431cb59b909b03252b9c164fad59" +dependencies = [ + "autocfg 1.0.1", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d41702bd167c2df5520b384281bc111a4b5efcf7fbc4c9c222c815b07e0a6a6a" +dependencies = [ + "autocfg 1.0.1", + "num-bigint", + "num-integer", + "num-traits", +] + [[package]] name = "num-traits" version = "0.2.14" @@ -1303,6 +1465,15 @@ version = "3.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "6acbef58a60fe69ab50510a55bc8cdd4d6cf2283d27ad338f54cb52747a9cf2d" +[[package]] +name = "output_vt100" +version = "0.1.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53cdc5b785b7a58c5aad8216b3dfa114df64b0b06ae6e1501cef91df2fbdf8f9" +dependencies = [ + "winapi 0.3.9", +] + [[package]] name = "parking_lot" version = "0.11.2" @@ -1328,6 +1499,12 @@ dependencies = [ "winapi 0.3.9", ] +[[package]] +name = "paste" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "acbf547ad0c65e31259204bd90935776d1c693cec2f4ff7abb7a1bbbd40dfe58" + [[package]] name = "pear" version = "0.2.3" @@ -1407,6 +1584,18 @@ version = "0.2.10" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ac74c624d6b2d21f425f752262f42188365d7b8ff1aff74c82e45136510a4857" +[[package]] +name = "pretty_assertions" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec0cfe1b2403f172ba0f234e500906ee0a3e493fb81092dac23ebefe129301cc" +dependencies = [ + "ansi_term", + "ctor", + "diff", + "output_vt100", +] + [[package]] name = "proc-macro-error" version = "1.0.4" @@ -1693,6 +1882,12 @@ dependencies = [ "rand_core 0.6.3", ] +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + [[package]] name = "rdrand" version = "0.4.0" @@ -1737,6 +1932,8 @@ version = "1.5.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d07a8629359eb56f1e2fb1652bb04212c072a87ba68546a04065d525673ac461" dependencies = [ + "aho-corasick", + "memchr", "regex-syntax", ] @@ -2044,6 +2241,15 @@ version = "1.0.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "71d301d4193d031abdd79ff7e3dd721168a9572ef3fe51a1517aba235bd8f86e" +[[package]] +name = "safe_arch" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "794821e4ccb0d9f979512f9c1973480123f9bd62a90d74ab0f9426fcf8f4a529" +dependencies = [ + "bytemuck", +] + [[package]] name = "same-file" version = "1.0.6" @@ -2269,6 +2475,19 @@ dependencies = [ "libc", ] +[[package]] +name = "simba" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0b7840f121a46d63066ee7a99fc81dcabbc6105e437cae43528cea199b5a05f" +dependencies = [ + "approx 0.5.0", + "num-complex", + "num-traits", + "paste", + "wide", +] + [[package]] name = "slab" version = "0.4.4" @@ -2277,9 +2496,9 @@ checksum = "c307a32c1c5c437f38c7fd45d753050587732ba8628319fbdf12a7e289ccc590" [[package]] name = "smallvec" -version = "1.6.1" +version = "1.7.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fe0f37c9e8f3c5a4a66ad655a93c74daac4ad00c441533bf5c6e7990bb42604e" +checksum = "1ecab6c735a6bb4139c0caafd0cc3635748bbb3acf4550e8138122099251f309" [[package]] name = "socket2" @@ -2477,9 +2696,9 @@ checksum = "6bdef32e8150c2a081110b42772ffe7d7c9032b606bc226c8260fd97e0976601" [[package]] name = "syn" -version = "1.0.77" +version = "1.0.78" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5239bc68e0fef57495900cfea4e8dc75596d9a319d7e16b1e0a440d24e6fe0a0" +checksum = "a4eac2e6c19f5c3abc0c229bea31ff0b9b091c7b14990e8924b92902a303a0c0" dependencies = [ "proc-macro2", "quote", @@ -2601,9 +2820,9 @@ dependencies = [ [[package]] name = "tokio-macros" -version = "1.3.0" +version = "1.4.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "54473be61f4ebe4efd09cec9bd5d16fa51d70ea0192213d754d2d500457db110" +checksum = "154794c8f499c2619acd19e839294703e9e32e7630ef5f46ea80d4ef0fbee5eb" dependencies = [ "proc-macro2", "quote", @@ -3051,14 +3270,24 @@ dependencies = [ [[package]] name = "whoami" -version = "1.1.4" +version = "1.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cabfe22aa4936611957e0b5ad9ed0472ac52b2bfb9aedac4a3f3a91a03bd1ff0" +checksum = "483a59fee1a93fec90eb08bc2eb4315ef10f4ebc478b3a5fadc969819cb66117" dependencies = [ "wasm-bindgen", "web-sys", ] +[[package]] +name = "wide" +version = "0.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2f2548e954f6619da26c140d020e99e59a2ca872a11f1e6250b829e8c96c893" +dependencies = [ + "bytemuck", + "safe_arch", +] + [[package]] name = "winapi" version = "0.2.8" diff --git a/cfd_protocol/src/interval.rs b/cfd_protocol/src/interval.rs index f0e41c5..30bdab5 100644 --- a/cfd_protocol/src/interval.rs +++ b/cfd_protocol/src/interval.rs @@ -14,7 +14,7 @@ const MAX_DIGITS: usize = 20; const BASE: usize = 2; /// Binary representation of a price interval. -#[derive(Clone, Debug)] +#[derive(Clone, Debug, PartialEq)] pub struct Digits(BitVec); impl Digits { diff --git a/cfd_protocol/src/protocol.rs b/cfd_protocol/src/protocol.rs index 16a643d..6c8a7fc 100644 --- a/cfd_protocol/src/protocol.rs +++ b/cfd_protocol/src/protocol.rs @@ -385,7 +385,7 @@ impl PartialEq for Announcement { } } -#[derive(Debug, Clone)] +#[derive(Debug, Clone, PartialEq)] pub struct Payout { digits: interval::Digits, maker_amount: Amount, diff --git a/daemon/Cargo.toml b/daemon/Cargo.toml index 7b259f1..6f37d7e 100644 --- a/daemon/Cargo.toml +++ b/daemon/Cargo.toml @@ -14,6 +14,11 @@ clap = "3.0.0-beta.4" futures = { version = "0.3", default-features = false } hex = "0.4" hkdf = "0.11" +itertools = "0.10" +nalgebra = { version = "0.29", default-features = false, features = ["std"] } +ndarray = "0.15.3" +ndarray_einsum_beta = "0.7.0" +num = "0.4.0" rand = "0.6" reqwest = { version = "0.11", default-features = false, features = ["json", "rustls-tls-webpki-roots"] } rocket = { version = "0.5.0-rc.1", features = ["json"] } @@ -47,5 +52,6 @@ name = "maker" path = "src/maker.rs" [dev-dependencies] +pretty_assertions = "1" tempfile = "3" time = { version = "0.3", features = ["std"] } diff --git a/daemon/src/payout_curve.rs b/daemon/src/payout_curve.rs index ea86576..f5e2b98 100644 --- a/daemon/src/payout_curve.rs +++ b/daemon/src/payout_curve.rs @@ -1,30 +1,744 @@ +use std::fmt; + use crate::model::{Leverage, Usd}; -use anyhow::Result; +use crate::payout_curve::curve::Curve; +use anyhow::{Context, Result}; use bdk::bitcoin; -use cfd_protocol::interval::MAX_PRICE_DEC; use cfd_protocol::{generate_payouts, Payout}; +use itertools::Itertools; +use ndarray::prelude::*; +use num::{FromPrimitive, ToPrimitive}; +use rust_decimal::Decimal; + +mod basis; +mod basis_eval; +mod compat; +mod csr_tools; +mod curve; +mod curve_factory; +mod splineobject; +mod utils; + +/// function to generate an iterator of values, heuristically viewed as: +/// +/// `[left_price_boundary, right_price_boundary], maker_payout_value` +/// +/// with units +/// +/// `[Usd, Usd], bitcoin::Amount` +/// +/// A key item to note is that although the POC logic has been to imposed +/// that maker goes short every time, there is no reason to make the math +/// have this imposition as well. As such, the `long_position` parameter +/// is used to indicate which party (Maker or Taker) has the long position, +/// and everything else is handled internally. +/// +/// As well, the POC has also demanded that the Maker always has unity +/// leverage, hence why the ability to to specify this amount has been +/// omitted from the parameters. Internally, it is hard-coded to unity +/// in the call to PayoutCurve::new(), so this behaviour can be changed in +/// the future trivially. +/// +/// ### Paramters +/// +/// * price: BTC-USD exchange rate used to create CFD contract +/// * quantity: Interger number of one-dollar USD contracts contained in the +/// CFD; expressed as a Usd amount +/// * leverage: Leveraging used by the taker +/// +/// ### Returns +/// +/// The list of [`Payout`]s for the given price, quantity and leverage. +pub fn calculate(price: Usd, quantity: Usd, leverage: Leverage) -> Result> { + let payouts = calculate_payout_parameters(price, quantity, leverage)? + .into_iter() + .map(PayoutParameter::into_payouts) + .flatten_ok() + .collect::>>()?; + + Ok(payouts) +} + +const CONTRACT_VALUE: f64 = 1.; +const N_PAYOUTS: usize = 200; +const SHORT_LEVERAGE: usize = 1; -pub fn calculate( +/// Internal calculate function for the payout curve. +/// +/// To ease testing, we write our tests against this function because it has a more human-friendly +/// output. The design goal here is that the the above `calculate` function is as thin as possible. +fn calculate_payout_parameters( price: Usd, - _quantity: Usd, - maker_payin: bitcoin::Amount, - (taker_payin, _leverage): (bitcoin::Amount, Leverage), -) -> Result> { - let dollars = price.try_into_u64()?; - let payouts = vec![ - generate_payouts( - 0..=(dollars - 10), - maker_payin + taker_payin, - bitcoin::Amount::ZERO, - )?, - generate_payouts((dollars - 9)..=(dollars + 10), maker_payin, taker_payin)?, + quantity: Usd, + long_leverage: Leverage, +) -> Result> { + let initial_rate = price + .try_into_u64() + .context("Cannot convert price to u64")? as f64; + let quantity = quantity + .try_into_u64() + .context("Cannot convert quantity to u64")? as usize; + + let payout_curve = PayoutCurve::new( + initial_rate as f64, + long_leverage.0 as usize, + SHORT_LEVERAGE, + quantity, + CONTRACT_VALUE, + None, + )?; + + let payout_parameters = payout_curve + .generate_payout_scheme(N_PAYOUTS)? + .rows() + .into_iter() + .map(|row| { + let left_bound = row[0] as u64; + let right_bound = row[1] as u64; + let long_amount = row[2]; + + let short_amount = to_sats(payout_curve.total_value - long_amount)?; + let long_amount = to_sats(long_amount)?; + + Ok(PayoutParameter { + left_bound, + right_bound, + long_amount, + short_amount, + }) + }) + .collect::>>()?; + + Ok(payout_parameters) +} + +#[derive(PartialEq)] +struct PayoutParameter { + left_bound: u64, + right_bound: u64, + long_amount: u64, + short_amount: u64, +} + +impl PayoutParameter { + fn into_payouts(self) -> Result> { generate_payouts( - (dollars + 11)..=MAX_PRICE_DEC, - bitcoin::Amount::ZERO, - maker_payin + taker_payin, - )?, - ] - .concat(); + self.left_bound..=self.right_bound, + bitcoin::Amount::from_sat(self.short_amount), + bitcoin::Amount::from_sat(self.long_amount), + ) + } +} - Ok(payouts) +impl fmt::Debug for PayoutParameter { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!( + f, + "payout({}..={}, {}, {})", + self.left_bound, self.right_bound, self.short_amount, self.long_amount + ) + } +} + +/// Converts a float with any precision to a [`bitcoin::Amount`]. +fn to_sats(btc: f64) -> Result { + let sats_per_btc = Decimal::from(100_000_000); + + let btc = Decimal::from_f64(btc).context("Cannot create decimal from float")?; + let sats = btc * sats_per_btc; + let sats = sats.to_u64().context("Cannot fit sats into u64")?; + + Ok(sats) +} + +#[derive(thiserror::Error, Debug)] +pub enum Error { + #[error("failed to init CSR object--is the specified shape correct?")] + #[allow(clippy::upper_case_acronyms)] + CannotInitCSR, + #[error("matrix must be square")] + MatrixMustBeSquare, + #[error("evaluation outside parametric domain")] + InvalidDomain, + #[error("einsum error--array size mismatch?")] + Einsum, + #[error("no operand string found")] + NoEinsumOperatorString, + #[error("cannot connect periodic curves")] + CannotConnectPeriodicCurves, + #[error("degree must be strictly positive")] + DegreeMustBePositive, + #[error("all parameter arrays must have the same length if not using a tensor grid")] + InvalidDerivative, + #[error("Rational derivative not implemented for order sum(d) > 1")] + DerivativeNotImplemented, + #[error("requested segmentation is too coarse for this curve")] + InvalidSegmentation, + #[error("concatonation error")] + NdArray { + #[from] + source: ndarray::ShapeError, + }, + #[error(transparent)] + NotOneDimensional { + #[from] + source: compat::NotOneDimensional, + }, +} + +#[derive(Clone, Debug)] +struct PayoutCurve { + curve: Curve, + has_upper_limit: bool, + lower_corner: f64, + upper_corner: f64, + total_value: f64, +} + +impl PayoutCurve { + fn new( + initial_rate: f64, + leverage_long: usize, + leverage_short: usize, + n_contracts: usize, + contract_value: f64, + tolerance: Option, + ) -> Result { + let tolerance = tolerance.unwrap_or(1e-6); + let bounds = cutoffs(initial_rate, leverage_long, leverage_short); + let total_value = pool_value( + initial_rate, + n_contracts, + contract_value, + leverage_long, + leverage_short, + ); + let mut curve = curve_factory::line((0., 0.), (bounds.0, 0.), false)?; + + let payout = + create_long_payout_function(initial_rate, n_contracts, contract_value, leverage_long); + let variable_payout = + curve_factory::fit(payout, bounds.0, bounds.1, Some(tolerance), None)?; + curve.append(variable_payout)?; + + let upper_corner; + if bounds.2 { + let upper_liquidation = curve_factory::line( + (bounds.1, total_value), + (4. * initial_rate, total_value), + false, + )?; + curve.append(upper_liquidation)?; + upper_corner = bounds.1; + } else { + upper_corner = curve.spline.bases[0].end(); + } + + Ok(PayoutCurve { + curve, + has_upper_limit: bounds.2, + lower_corner: bounds.0, + upper_corner, + total_value, + }) + } + + pub fn generate_payout_scheme(&self, n_segments: usize) -> Result, Error> { + let n_min; + if self.has_upper_limit { + n_min = 3; + } else { + n_min = 2; + } + + if n_segments < n_min { + return Result::Err(Error::InvalidSegmentation); + } + + let t; + if self.has_upper_limit { + t = self.build_sampling_vector_upper_bounded(n_segments); + } else { + t = self.build_sampling_vector_upper_unbounded(n_segments) + } + + let mut z_arr = self.curve.evaluate(&mut &[t][..])?; + if self.has_upper_limit { + self.modify_samples_bounded(&mut z_arr); + } else { + self.modify_samples_unbounded(&mut z_arr); + } + self.generate_segments(&mut z_arr); + + Ok(z_arr) + } + + fn build_sampling_vector_upper_bounded(&self, n_segs: usize) -> Array1 { + let knots = &self.curve.spline.knots(0, None).unwrap()[0]; + let klen = knots.len(); + let n_64 = (n_segs + 1) as f64; + let d = knots[klen - 2] - knots[1]; + let delta_0 = d / (2. * (n_64 - 5.)); + let delta_1 = d * (n_64 - 6.) / ((n_64 - 5.) * (n_64 - 4.)); + + let mut vec = Vec::::with_capacity(n_segs + 2); + for i in 0..n_segs + 2 { + if i == 0 { + vec.push(self.curve.spline.bases[0].start()); + } else if i == 1 { + vec.push(knots[1]); + } else if i == 2 { + vec.push(knots[1] + delta_0); + } else if i == n_segs - 1 { + vec.push(knots[klen - 2] - delta_0); + } else if i == n_segs { + vec.push(knots[klen - 2]); + } else if i == n_segs + 1 { + vec.push(self.curve.spline.bases[0].end()); + } else { + let c = (i - 2) as f64; + vec.push(knots[1] + delta_0 + c * delta_1); + } + } + Array1::::from_vec(vec) + } + + fn build_sampling_vector_upper_unbounded(&self, n_segs: usize) -> Array1 { + let knots = &self.curve.spline.knots(0, None).unwrap()[0]; + let klen = knots.len(); + let n_64 = (n_segs + 1) as f64; + let d = knots[klen - 1] - knots[1]; + let delta = d / (n_64 - 1_f64); + let delta_x = d / (2. * (n_64 - 1_f64)); + let delta_y = 3. * d / (2. * (n_64 - 1_f64)); + + let mut vec = Vec::::with_capacity(n_segs + 2); + for i in 0..n_segs + 2 { + if i == 0 { + vec.push(self.curve.spline.bases[0].start()); + } else if i == 1 { + vec.push(knots[1]); + } else if i == 2 { + vec.push(knots[1] + delta_x); + } else if i == n_segs { + vec.push(knots[klen - 1] - delta_y); + } else if i == n_segs + 1 { + vec.push(knots[klen - 1]); + } else { + let c = (i - 2) as f64; + vec.push(knots[1] + delta_x + c * delta); + } + } + Array1::::from_vec(vec) + } + + fn modify_samples_bounded(&self, arr: &mut Array2) { + let n = arr.shape()[0]; + let capacity = 2 * (n - 2); + let mut vec = Vec::::with_capacity(2 * capacity); + for (i, e) in arr.slice(s![.., 0]).iter().enumerate() { + if i < 2 || i > n - 3 { + vec.push(*e); + } else if i == 2 { + vec.push(arr[[i - 1, 0]]); + vec.push(arr[[i, 1]]); + vec.push((*e + arr[[i + 1, 0]]) / 2.); + } else if i == n - 3 { + vec.push((arr[[i - 1, 0]] + *e) / 2.); + vec.push(arr[[i, 1]]); + vec.push(arr[[i + 1, 0]]); + } else { + vec.push((arr[[i - 1, 0]] + *e) / 2.); + vec.push(arr[[i, 1]]); + vec.push((*e + arr[[i + 1, 0]]) / 2.); + } + vec.push(arr[[i, 1]]); + } + + *arr = Array2::::from_shape_vec((capacity, 2), vec).unwrap(); + } + + fn modify_samples_unbounded(&self, arr: &mut Array2) { + let n = arr.shape()[0]; + let capacity = 2 * (n - 1); + let mut vec = Vec::::with_capacity(2 * capacity); + for (i, e) in arr.slice(s![.., 0]).iter().enumerate() { + if i < 2 { + vec.push(*e); + } else if i == 2 { + vec.push(arr[[i - 1, 0]]); + vec.push(arr[[i, 1]]); + vec.push((*e + arr[[i + 1, 0]]) / 2.); + } else if i == n - 1 { + vec.push((arr[[i - 1, 0]] + *e) / 2.); + vec.push(arr[[i, 1]]); + vec.push(arr[[i, 0]]); + } else { + vec.push((arr[[i - 1, 0]] + *e) / 2.); + vec.push(arr[[i, 1]]); + vec.push((*e + arr[[i + 1, 0]]) / 2.); + } + vec.push(arr[[i, 1]]); + } + + *arr = Array2::::from_shape_vec((capacity, 2), vec).unwrap(); + } + + /// this should only be used on an array `arr` that has been + /// processed by self.modify_samples_* first, otherwise the results + /// will be jibberish. + fn generate_segments(&self, arr: &mut Array2) { + let capacity = 3 * arr.shape()[0] / 2; + let mut vec = Vec::::with_capacity(capacity); + for (i, e) in arr.slice(s![.., 0]).iter().enumerate() { + if i == 0 { + vec.push(e.floor()); + } else if i % 2 == 1 { + vec.push(e.round()); + vec.push(arr[[i, 1]]); + } else { + vec.push(e.round() + 1_f64); + } + } + + *arr = Array2::::from_shape_vec((capacity / 3, 3), vec).unwrap(); + } +} + +fn cutoffs(initial_rate: f64, leverage_long: usize, leverage_short: usize) -> (f64, f64, bool) { + let ll_64 = leverage_long as f64; + let ls_64 = leverage_short as f64; + let a = initial_rate * ll_64 / (ll_64 + 1_f64); + if leverage_short == 1 { + let b = 2. * initial_rate; + return (a, b, false); + } + let b = initial_rate * ls_64 / (ls_64 - 1_f64); + + (a, b, true) +} + +fn pool_value( + initial_rate: f64, + n_contracts: usize, + contract_value: f64, + leverage_long: usize, + leverage_short: usize, +) -> f64 { + let ll_64 = leverage_long as f64; + let ls_64 = leverage_short as f64; + let n_64 = n_contracts as f64; + + (n_64 * contract_value / initial_rate) * (1_f64 / ll_64 + 1_f64 / ls_64) +} + +fn create_long_payout_function( + initial_rate: f64, + n_contracts: usize, + contract_value: f64, + leverage_long: usize, +) -> impl Fn(&Array1) -> Array2 { + let n_64 = n_contracts as f64; + let ll_64 = leverage_long as f64; + + move |t: &Array1| { + let mut vec = Vec::::with_capacity(2 * t.len()); + for e in t.iter() { + let eval = (n_64 * contract_value) + * (1_f64 / (initial_rate * ll_64) + (1_f64 / initial_rate - 1_f64 / e)); + vec.push(*e); + vec.push(eval); + } + + Array2::::from_shape_vec((t.len(), 2), vec).unwrap() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use rust_decimal_macros::dec; + use std::ops::RangeInclusive; + + #[test] + fn test_bounded() { + let initial_rate = 40000.0; + let leverage_long = 5; + let leverage_short = 2; + let n_contracts = 200; + let contract_value = 100.; + + let payout = PayoutCurve::new( + initial_rate, + leverage_long, + leverage_short, + n_contracts, + contract_value, + None, + ) + .unwrap(); + + let z = payout.generate_payout_scheme(5000).unwrap(); + + assert!(z.shape()[0] == 5000); + } + + #[test] + fn test_unbounded() { + let initial_rate = 40000.0; + let leverage_long = 5; + let leverage_short = 1; + let n_contracts = 200; + let contract_value = 100.; + + let payout = PayoutCurve::new( + initial_rate, + leverage_long, + leverage_short, + n_contracts, + contract_value, + None, + ) + .unwrap(); + + let z = payout.generate_payout_scheme(5000).unwrap(); + + // out-by-one error expected at this point in time + assert!(z.shape()[0] == 5001); + } + + #[test] + fn calculate_snapshot() { + let actual_payouts = + calculate_payout_parameters(Usd(dec!(54000.00)), Usd(dec!(3500.00)), Leverage(5)) + .unwrap(); + + let expected_payouts = vec![ + payout(0..=45000, 7777777, 0), + payout(45001..=45315, 7750759, 27018), + payout(45316..=45630, 7697244, 80533), + payout(45631..=45945, 7644417, 133359), + payout(45946..=46260, 7592270, 185507), + payout(46261..=46575, 7540793, 236984), + payout(46576..=46890, 7489978, 287799), + payout(46891..=47205, 7439816, 337961), + payout(47206..=47520, 7390298, 387479), + payout(47521..=47835, 7341415, 436362), + payout(47836..=48150, 7293159, 484618), + payout(48151..=48465, 7245520, 532257), + payout(48466..=48780, 7198490, 579287), + payout(48781..=49095, 7152060, 625717), + payout(49096..=49410, 7106222, 671555), + payout(49411..=49725, 7060965, 716812), + payout(49726..=50040, 7016282, 761494), + payout(50041..=50355, 6972164, 805612), + payout(50356..=50670, 6928602, 849174), + payout(50671..=50985, 6885587, 892189), + payout(50986..=51300, 6843111, 934666), + payout(51301..=51615, 6801163, 976613), + payout(51616..=51930, 6759737, 1018040), + payout(51931..=52245, 6718822, 1058955), + payout(52246..=52560, 6678410, 1099367), + payout(52561..=52875, 6638493, 1139284), + payout(52876..=53190, 6599060, 1178716), + payout(53191..=53505, 6560105, 1217672), + payout(53506..=53820, 6521617, 1256160), + payout(53821..=54135, 6483588, 1294189), + payout(54136..=54450, 6446009, 1331768), + payout(54451..=54765, 6408872, 1368905), + payout(54766..=55080, 6372166, 1405610), + payout(55081..=55395, 6335885, 1441892), + payout(55396..=55710, 6300018, 1477758), + payout(55711..=56025, 6264558, 1513219), + payout(56026..=56340, 6229494, 1548282), + payout(56341..=56655, 6194820, 1582957), + payout(56656..=56970, 6160524, 1617253), + payout(56971..=57285, 6126599, 1651177), + payout(57286..=57600, 6093037, 1684740), + payout(57601..=57915, 6059827, 1717949), + payout(57916..=58230, 6026965, 1750812), + payout(58231..=58545, 5994445, 1783332), + payout(58546..=58860, 5962264, 1815512), + payout(58861..=59175, 5930419, 1847358), + payout(59176..=59490, 5898905, 1878872), + payout(59491..=59805, 5867718, 1910059), + payout(59806..=60120, 5836855, 1940922), + payout(60121..=60435, 5806311, 1971465), + payout(60436..=60750, 5776084, 2001693), + payout(60751..=61065, 5746168, 2031608), + payout(61066..=61380, 5716561, 2061216), + payout(61381..=61695, 5687258, 2090519), + payout(61696..=62010, 5658255, 2119522), + payout(62011..=62325, 5629549, 2148228), + payout(62326..=62640, 5601135, 2176642), + payout(62641..=62955, 5573010, 2204767), + payout(62956..=63270, 5545170, 2232607), + payout(63271..=63585, 5517611, 2260165), + payout(63586..=63900, 5490330, 2287447), + payout(63901..=64215, 5463321, 2314455), + payout(64216..=64530, 5436583, 2341194), + payout(64531..=64845, 5410109, 2367667), + payout(64846..=65160, 5383898, 2393879), + payout(65161..=65475, 5357944, 2419833), + payout(65476..=65790, 5332245, 2445532), + payout(65791..=66105, 5306795, 2470982), + payout(66106..=66420, 5281592, 2496185), + payout(66421..=66735, 5256631, 2521146), + payout(66736..=67050, 5231909, 2545868), + payout(67051..=67365, 5207421, 2570356), + payout(67366..=67680, 5183164, 2594612), + payout(67681..=67995, 5159135, 2618642), + payout(67996..=68310, 5135328, 2642449), + payout(68311..=68625, 5111740, 2666037), + payout(68626..=68940, 5088368, 2689409), + payout(68941..=69255, 5065207, 2712569), + payout(69256..=69570, 5042254, 2735523), + payout(69571..=69885, 5019505, 2758272), + payout(69886..=70200, 4996955, 2780821), + payout(70201..=70515, 4974602, 2803175), + payout(70516..=70830, 4952442, 2825335), + payout(70831..=71145, 4930473, 2847304), + payout(71146..=71460, 4908694, 2869083), + payout(71461..=71775, 4887102, 2890675), + payout(71776..=72090, 4865695, 2912081), + payout(72091..=72405, 4844473, 2933304), + payout(72406..=72720, 4823433, 2954344), + payout(72721..=73035, 4802573, 2975204), + payout(73036..=73350, 4781891, 2995886), + payout(73351..=73665, 4761385, 3016391), + payout(73666..=73980, 4741054, 3036722), + payout(73981..=74295, 4720896, 3056881), + payout(74296..=74610, 4700909, 3076868), + payout(74611..=74925, 4681090, 3096686), + payout(74926..=75240, 4661439, 3116338), + payout(75241..=75555, 4641953, 3135824), + payout(75556..=75870, 4622630, 3155146), + payout(75871..=76185, 4603469, 3174307), + payout(76186..=76500, 4584468, 3193309), + payout(76501..=76815, 4565624, 3212153), + payout(76816..=77130, 4546937, 3230840), + payout(77131..=77445, 4528403, 3249374), + payout(77446..=77760, 4510022, 3267755), + payout(77761..=78075, 4491791, 3285986), + payout(78076..=78390, 4473708, 3304068), + payout(78391..=78705, 4455773, 3322004), + payout(78706..=79020, 4437982, 3339795), + payout(79021..=79335, 4420333, 3357443), + payout(79336..=79650, 4402827, 3374950), + payout(79651..=79965, 4385459, 3392318), + payout(79966..=80280, 4368228, 3409548), + payout(80281..=80595, 4351133, 3426643), + payout(80596..=80910, 4334172, 3443605), + payout(80911..=81225, 4317343, 3460434), + payout(81226..=81540, 4300643, 3477134), + payout(81541..=81855, 4284071, 3493705), + payout(81856..=82170, 4267626, 3510151), + payout(82171..=82485, 4251305, 3526472), + payout(82486..=82800, 4235107, 3542670), + payout(82801..=83115, 4219029, 3558748), + payout(83116..=83430, 4203070, 3574707), + payout(83431..=83745, 4187229, 3590547), + payout(83746..=84060, 4171506, 3606271), + payout(84061..=84375, 4155899, 3621878), + payout(84376..=84690, 4140406, 3637371), + payout(84691..=85005, 4125028, 3652749), + payout(85006..=85320, 4109763, 3668014), + payout(85321..=85635, 4094610, 3683167), + payout(85636..=85950, 4079567, 3698209), + payout(85951..=86265, 4064635, 3713142), + payout(86266..=86580, 4049812, 3727965), + payout(86581..=86895, 4035096, 3742680), + payout(86896..=87210, 4020488, 3757289), + payout(87211..=87525, 4005985, 3771792), + payout(87526..=87840, 3991587, 3786189), + payout(87841..=88155, 3977293, 3800484), + payout(88156..=88470, 3963102, 3814675), + payout(88471..=88785, 3949013, 3828764), + payout(88786..=89100, 3935024, 3842753), + payout(89101..=89415, 3921135, 3856642), + payout(89416..=89730, 3907344, 3870432), + payout(89731..=90045, 3893652, 3884125), + payout(90046..=90360, 3880056, 3897721), + payout(90361..=90675, 3866555, 3911221), + payout(90676..=90990, 3853150, 3924627), + payout(90991..=91305, 3839837, 3937940), + payout(91306..=91620, 3826618, 3951159), + payout(91621..=91935, 3813489, 3964287), + payout(91936..=92250, 3800452, 3977325), + payout(92251..=92565, 3787504, 3990273), + payout(92566..=92880, 3774644, 4003133), + payout(92881..=93195, 3761872, 4015905), + payout(93196..=93510, 3749186, 4028591), + payout(93511..=93825, 3736585, 4041192), + payout(93826..=94140, 3724069, 4053708), + payout(94141..=94455, 3711636, 4066141), + payout(94456..=94770, 3699286, 4078491), + payout(94771..=95085, 3687016, 4090760), + payout(95086..=95400, 3674827, 4102949), + payout(95401..=95715, 3662718, 4115059), + payout(95716..=96030, 3650686, 4127091), + payout(96031..=96345, 3638733, 4139044), + payout(96346..=96660, 3626856, 4150920), + payout(96661..=96975, 3615057, 4162720), + payout(96976..=97290, 3603333, 4174444), + payout(97291..=97605, 3591684, 4186093), + payout(97606..=97920, 3580110, 4197666), + payout(97921..=98235, 3568611, 4209166), + payout(98236..=98550, 3557184, 4220593), + payout(98551..=98865, 3545831, 4231946), + payout(98866..=99180, 3534550, 4243227), + payout(99181..=99495, 3523340, 4254437), + payout(99496..=99810, 3512201, 4265576), + payout(99811..=100125, 3501133, 4276644), + payout(100126..=100440, 3490134, 4287643), + payout(100441..=100755, 3479205, 4298572), + payout(100756..=101070, 3468344, 4309433), + payout(101071..=101385, 3457551, 4320226), + payout(101386..=101700, 3446825, 4330952), + payout(101701..=102015, 3436165, 4341611), + payout(102016..=102330, 3425572, 4352205), + payout(102331..=102645, 3415044, 4362732), + payout(102646..=102960, 3404581, 4373196), + payout(102961..=103275, 3394182, 4383594), + payout(103276..=103590, 3383847, 4393930), + payout(103591..=103905, 3373574, 4404202), + payout(103906..=104220, 3363364, 4414412), + payout(104221..=104535, 3353216, 4424561), + payout(104536..=104850, 3343128, 4434648), + payout(104851..=105165, 3333101, 4444675), + payout(105166..=105480, 3323134, 4454643), + payout(105481..=105795, 3313226, 4464551), + payout(105796..=106110, 3303377, 4474400), + payout(106111..=106425, 3293585, 4484192), + payout(106426..=106740, 3283851, 4493926), + payout(106741..=107055, 3274174, 4503603), + payout(107056..=107370, 3264552, 4513225), + payout(107371..=107764, 3254986, 4522790), + payout(107765..=108000, 3240740, 4537037), + ]; + + pretty_assertions::assert_eq!(actual_payouts, expected_payouts); + } + + #[test] + fn verfiy_tails() { + let actual_payouts = + calculate_payout_parameters(Usd(dec!(54000.00)), Usd(dec!(3500.00)), Leverage(5)) + .unwrap(); + + let lower_tail = payout(0..=45000, 7777777, 0); + let upper_tail = payout(107765..=108000, 3240740, 4537037); + + pretty_assertions::assert_eq!(actual_payouts.first().unwrap(), &lower_tail); + pretty_assertions::assert_eq!(actual_payouts.last().unwrap(), &upper_tail); + } + + fn payout(range: RangeInclusive, short: u64, long: u64) -> PayoutParameter { + PayoutParameter { + left_bound: *range.start(), + right_bound: *range.end(), + long_amount: long, + short_amount: short, + } + } } diff --git a/daemon/src/payout_curve/README.md b/daemon/src/payout_curve/README.md new file mode 100644 index 0000000..c4db040 --- /dev/null +++ b/daemon/src/payout_curve/README.md @@ -0,0 +1,3 @@ +## Note + +This codebase is effectively a brute-force copy of the the [Splipy](https://github.com/SINTEF/Splipy) python lib, cherry-picked to do what we need while still remaining consistent with the source material. This is quite a bit of overkill for our purposes, so this will either be refined in the futre **or** a complete translation of the Splipy lib will be created as it's own crate. As such, this code should be considered transitional at best. diff --git a/daemon/src/payout_curve/basis.rs b/daemon/src/payout_curve/basis.rs new file mode 100644 index 0000000..21c4af1 --- /dev/null +++ b/daemon/src/payout_curve/basis.rs @@ -0,0 +1,208 @@ +use crate::payout_curve::basis_eval::*; +use crate::payout_curve::csr_tools::CSR; +use crate::payout_curve::utils::*; +use crate::payout_curve::Error; + +use core::cmp::max; +use ndarray::prelude::*; +use ndarray::{concatenate, s}; + +#[derive(Clone, Debug)] +pub struct BSplineBasis { + pub knots: Array1, + pub order: usize, + pub periodic: isize, + pub knot_tol: f64, +} + +impl BSplineBasis { + pub fn new( + order: Option, + knots: Option>, + periodic: Option, + ) -> Result { + let order = order.unwrap_or(2); + let periodic = periodic.unwrap_or(-1); + let knots = match knots { + Some(knots) => knots, + None => default_knot(order, periodic)?, + }; + let ktol = knot_tolerance(None); + + Ok(BSplineBasis { + order, + knots, + periodic, + knot_tol: ktol, + }) + } + + pub fn num_functions(&self) -> usize { + let p = (self.periodic + 1) as usize; + + self.knots.len() - self.order - p + } + + /// Start point of parametric domain. For open knot vectors, this is the + /// first knot. + pub fn start(&self) -> f64 { + self.knots[self.order - 1] + } + + /// End point of parametric domain. For open knot vectors, this is the + /// last knot. + pub fn end(&self) -> f64 { + self.knots[self.knots.len() - self.order] + } + + /// Fetch greville points, also known as knot averages + /// over entire knot vector: + /// .. math:: \\sum_{j=i+1}^{i+p-1} \\frac{t_j}{p-1} + pub fn greville(&self) -> Array1 { + let n = self.num_functions() as i32; + + (0..n).map(|idx| self.greville_single(idx)).collect() + } + + fn greville_single(&self, index: i32) -> f64 { + let p = self.order as i32; + let den = (self.order - 1) as f64; + + self.knots.slice(s![index + 1..index + p]).sum() / den + } + + /// Evaluate all basis functions in a given set of points. + /// ## parameters: + /// * t: The parametric coordinate(s) in which to evaluate + /// * d: Number of derivatives to compute + /// * from_right: true if evaluation should be done in the limit from above + /// ## returns: + /// * CSR (sparse) matrix N\[i,j\] of all basis functions j evaluated in all points j + pub fn evaluate(&self, t: &mut Array1, d: usize, from_right: bool) -> Result { + let basis = Basis::new( + self.order, + self.knots.clone().to_owned(), + Some(self.periodic), + Some(self.knot_tol), + ); + snap(t, &basis.knots, Some(basis.ktol)); + + if self.order <= d { + let csr = CSR::new( + Array1::::zeros(0), + Array1::::zeros(0), + Array1::::zeros(t.len() + 1), + (t.len(), self.num_functions()), + )?; + + return Ok(csr); + } + + let out = basis.evaluate(t, d, Some(from_right))?; + + Ok(out) + } + + /// Snap evaluation points to knots if they are sufficiently close + /// as given in by knot_tolerance. + /// + /// * t: The parametric coordinate(s) in which to evaluate + pub fn snap(&self, t: &mut Array1) { + snap(t, &self.knots, Some(self.knot_tol)) + } + + /// Create a knot vector with higher order. + /// + /// The continuity at the knots are kept unchanged by increasing their + /// multiplicities. + /// + /// ### parameters + /// * amount: relative (polynomial) degree to raise the basis function by + pub fn raise_order(&mut self, amount: usize) { + if amount > 0 { + let knot_spans_arr = self.knot_spans(true); + let knot_spans = knot_spans_arr.iter().collect::>(); + let temp = self.knots.clone(); + let mut knots = temp.iter().collect::>(); + + for _ in 0..amount { + knots.append(&mut knot_spans.clone()); + } + + let mut knots_vec = knots.iter().map(|e| **e).collect::>(); + knots_vec.sort_by(cmp_f64); + + let knots_arr = Array1::::from_vec(knots_vec); + + let new_knot; + if self.periodic > -1 { + let n_0 = bisect_left(&knots_arr, &self.start(), knots_arr.len()); + let n_1 = + knot_spans.len() - bisect_left(&knots_arr, &self.end(), knots_arr.len()) - 1; + let mut new_knot_vec = knots[n_0 * amount..n_1 * amount] + .iter() + .map(|e| **e) + .collect::>(); + new_knot_vec.sort_by(cmp_f64); + new_knot = Array1::::from_vec(new_knot_vec); + } else { + new_knot = knots_arr; + } + + self.order += amount; + self.knots = new_knot; + } + } + + /// Return the set of unique knots in the knot vector. + /// + /// ### parameters: + /// * include_ghosts: if knots outside start/end are to be included. These \ + /// knots are used by periodic basis. + /// + /// ### returns: + /// * 1-D array of unique knots + pub fn knot_spans(&self, include_ghosts: bool) -> Array1 { + let p = &self.order; + + // TODO: this is VERY sloppy! + let mut res: Vec = vec![]; + if include_ghosts { + res.push(self.knots[0]); + for elem in self.knots.slice(s![1..]).iter() { + if (elem - res[res.len() - 1]).abs() > self.knot_tol { + res.push(*elem); + } + } + } else { + res.push(self.knots[p - 1]); + let klen = self.knots.len(); + for elem in self.knots.slice(s![p - 1..klen - p + 1]).iter() { + if (elem - res[res.len() - 1]).abs() > self.knot_tol { + res.push(*elem); + } + } + } + + Array1::::from_vec(res) + } +} + +fn default_knot(order: usize, periodic: isize) -> Result, Error> { + let prd = max(periodic, -1); + let p = (prd + 1) as usize; + let mut knots = concatenate( + Axis(0), + &[ + Array1::::zeros(order).view(), + Array1::::ones(order).view(), + ], + )?; + + for i in 0..p { + knots[i] = -1.; + knots[2 * order - i - 1] = 2.; + } + + Ok(knots) +} diff --git a/daemon/src/payout_curve/basis_eval.rs b/daemon/src/payout_curve/basis_eval.rs new file mode 100644 index 0000000..a5ce8a1 --- /dev/null +++ b/daemon/src/payout_curve/basis_eval.rs @@ -0,0 +1,177 @@ +use crate::payout_curve::csr_tools::CSR; +use crate::payout_curve::utils::*; +use crate::payout_curve::Error; + +use ndarray::prelude::*; +use std::cmp::min; + +#[derive(Clone, Debug)] +pub struct Basis { + pub knots: Array1, + pub order: usize, + pub periodic: isize, + n_all: usize, + n: usize, + pub start: f64, + pub end: f64, + pub ktol: f64, +} + +impl Basis { + pub fn new( + order: usize, + knots: Array1, + periodic: Option, + knot_tol: Option, + ) -> Self { + let p = periodic.unwrap_or(-1); + let n_all = knots.len() - order; + let n = knots.len() - order - ((p + 1) as usize); + let start = knots[order - 1]; + let end = knots[n_all]; + + Basis { + order, + knots, + periodic: p, + n_all, + n, + start, + end, + ktol: knot_tolerance(knot_tol), + } + } + + /// Wrap periodic evaluation into domain + fn wrap_periodic(&self, t: &mut Array1, right: &bool) { + for i in 0..t.len() { + if t[i] < self.start || t[i] > self.end { + t[i] = (t[i] - self.start) % (self.end - self.start) + self.start; + } + if (t[i] - self.start).abs() < self.ktol && !right { + t[i] = self.end; + } + } + } + + pub fn evaluate( + &self, + t: &mut Array1, + d: usize, + from_right: Option, + ) -> Result { + let m = t.len(); + let mut right = from_right.unwrap_or(true); + + if self.periodic >= 0 { + self.wrap_periodic(t, &right); + } + + let mut store = Array1::::zeros(self.order); + let mut mu; + let mut idx_j; + let mut idx_k; + + let mut data = Array1::::zeros(m * self.order); + let mut indices = Array1::::zeros(m * self.order); + let indptr = + Array1::::from_vec((0..m * self.order + 1).step_by(self.order).collect()); + + for i in 0..m { + right = from_right.unwrap_or(true); + // Special-case the endpoint, so the user doesn't need to + if (t[i] - self.end).abs() < self.ktol { + right = false; + } + // Skip non-periodic evaluation points outside the domain + if t[i] < self.start + || t[i] > self.end + || ((t[i] - self.start).abs() < self.ktol && !right) + { + continue; + } + + // mu = index of last non-zero basis function + if right { + mu = bisect_right(&self.knots, &t[i], self.n_all + self.order); + } else { + mu = bisect_left(&self.knots, &t[i], self.n_all + self.order); + } + mu = min(mu, self.n_all); + + for k in 0..self.order - 1 { + store[k] = 0.; + } + + // the last entry is a dummy-zero which is never used + store[self.order - 1] = 1.; + + for q in 1..self.order - d { + idx_j = self.order - q - 1; + idx_k = mu - q - 1; + store[idx_j] += store[idx_j + 1] * (self.knots[idx_k + q + 1] - t[i]) + / (self.knots[idx_k + q + 1] - self.knots[idx_k + 1]); + + for j in self.order - q..self.order - 1 { + // 'i'-index in global knot vector (ref Hughes book pg.21) + let k = mu - self.order + j; + store[j] = + store[j] * (t[i] - self.knots[k]) / (self.knots[k + q] - self.knots[k]); + store[j] += store[j + 1] * (self.knots[k + q + 1] - t[i]) + / (self.knots[k + q + 1] - self.knots[k + 1]); + } + idx_j = self.order - 1; + idx_k = mu - 1; + store[idx_j] = store[idx_j] * (t[i] - self.knots[idx_k]) + / (self.knots[idx_k + q] - self.knots[idx_k]); + } + + for q in self.order - d..self.order { + for j in self.order - q - 1..self.order { + // 'i'-index in global knot vector (ref Hughes book pg.21) + idx_k = mu - self.order + j; + if j != self.order - q - 1 { + store[j] = + store[j] * (q as f64) / (self.knots[idx_k + q] - self.knots[idx_k]); + } + if j != self.order - 1 { + store[j] -= store[j + 1] * (q as f64) + / (self.knots[idx_k + q + 1] - self.knots[idx_k + 1]); + } + } + } + + for (j, k) in (i * self.order..(i + 1) * self.order).enumerate() { + data[k] = store[j]; + indices[k] = (mu - self.order + j) % self.n; + } + } + + let csr = CSR::new(data, indices, indptr, (m, self.n))?; + + Ok(csr) + } +} + +pub fn knot_tolerance(tol: Option) -> f64 { + tol.unwrap_or(1e-10) +} + +/// Snap evaluation points to knots if they are sufficiently close +/// as specified by self.ktol +/// +/// * t: The parametric coordinate(s) in which to evaluate +/// * knots: knot-vector +/// * knot_tol: default=1e-10 +pub fn snap(t: &mut Array1, knots: &Array1, knot_tol: Option) { + let ktol = knot_tolerance(knot_tol); + + for j in 0..t.len() { + let i = bisect_left(knots, &t[j], knots.len()); + if i < knots.len() && (knots[i] - t[j]).abs() < ktol { + t[j] = knots[i]; + } else if i > 0 && (knots[i - 1] - t[j]).abs() < ktol { + t[j] = knots[i - 1]; + } + } +} diff --git a/daemon/src/payout_curve/compat.rs b/daemon/src/payout_curve/compat.rs new file mode 100644 index 0000000..6306ec6 --- /dev/null +++ b/daemon/src/payout_curve/compat.rs @@ -0,0 +1,69 @@ +use nalgebra::{ComplexField, DMatrix, Dynamic, Scalar}; +use ndarray::{Array1, Array2}; +use std::fmt::Debug; + +pub trait ToNAlgebraMatrix { + fn to_nalgebra_matrix(&self) -> DMatrix; +} + +pub trait To1DArray { + fn to_1d_array(&self) -> Result, NotOneDimensional>; +} + +pub trait To2DArray { + fn to_2d_array(&self) -> Array2; +} + +impl ToNAlgebraMatrix + for Array1 +{ + fn to_nalgebra_matrix(&self) -> DMatrix { + DMatrix::from_row_slice_generic( + Dynamic::new(self.len()), + Dynamic::new(1), + self.to_vec().as_slice(), + ) + } +} + +impl ToNAlgebraMatrix + for Array2 +{ + fn to_nalgebra_matrix(&self) -> DMatrix { + let flattened = self.rows().into_iter().fold(Vec::new(), |mut acc, next| { + acc.extend(next.into_iter().cloned()); + + acc + }); + + DMatrix::from_row_slice_generic( + Dynamic::new(self.nrows()), + Dynamic::new(self.ncols()), + flattened.as_slice(), + ) + } +} + +impl To1DArray for DMatrix { + fn to_1d_array(&self) -> Result, NotOneDimensional> { + if self.ncols() != 1 { + return Err(NotOneDimensional); + } + + Ok(Array1::from_shape_fn(self.nrows(), |index| { + self.get(index).unwrap().clone() + })) + } +} + +impl To2DArray for DMatrix { + fn to_2d_array(&self) -> Array2 { + Array2::from_shape_fn((self.nrows(), self.ncols()), |indices| { + self.get(indices).unwrap().clone() + }) + } +} + +#[derive(Debug, thiserror::Error)] +#[error("The provided matrix is not one-dimensional and cannot be converted into a 1D array")] +pub struct NotOneDimensional; diff --git a/daemon/src/payout_curve/csr_tools.rs b/daemon/src/payout_curve/csr_tools.rs new file mode 100644 index 0000000..89c0270 --- /dev/null +++ b/daemon/src/payout_curve/csr_tools.rs @@ -0,0 +1,152 @@ +use crate::payout_curve::compat::{To1DArray, ToNAlgebraMatrix}; +use crate::payout_curve::Error; +use ndarray::prelude::*; +use std::ops::Mul; + +/// NOTE: +/// This struct is provided here in this form as nalgebra_sparse +/// is rather embryonic and incluldes (basically) no solvers +/// at present. As we only need to be able construct as CSR +/// matrix and perform multiplication with a (dense) vector, it +/// seemed to make more sense to define our own rudementary CSR struct +/// and avoid the bloat of using a crate that will introduce +/// breaking changes regularly, and we need to write our own +/// solver regardless. +#[derive(Clone, Debug, PartialEq)] +#[allow(clippy::upper_case_acronyms)] +pub struct CSR { + pub data: Array1, + pub indices: Array1, + pub indptr: Array1, + pub shape: (usize, usize), + pub nnz: usize, +} + +impl CSR { + pub fn new( + data: Array1, + indices: Array1, + indptr: Array1, + shape: (usize, usize), + ) -> Result { + let major_dim: isize = (indptr.len() as isize) - 1; + let nnz = &data.len(); + + if major_dim > 1 && shape.0 as isize == major_dim { + Result::Ok(CSR { + data, + indices, + indptr, + shape, + nnz: *nnz, + }) + } else { + Result::Err(Error::CannotInitCSR) + } + } + + // matrix version of `solve()`; useful for solving AX = B. Implementation + // is horrible. + pub fn matrix_solve(&self, b_arr: &Array2) -> Result, Error> { + let a_arr = self.todense(); + let ncols = b_arr.shape()[1]; + let mut temp = (0..ncols) + .rev() + .map(|e| { + let b = b_arr.slice(s![.., e]).to_owned(); + + let sol = lu_solve(&a_arr, &b).unwrap(); + Ok(sol.to_vec()) + }) + .collect::, Error>>()?; + + let nrows = temp[0].len(); + let mut raveled = Vec::with_capacity(nrows * temp.len()); + + for _ in 0..nrows { + for vec in &mut temp { + raveled.push(vec.pop().unwrap()); + } + } + + raveled.reverse(); + + let result = Array2::::from_shape_vec((nrows, ncols), raveled)?.to_owned(); + + Ok(result) + } + + pub fn todense(&self) -> Array2 { + let mut out = Array2::::zeros(self.shape); + for i in 0..self.shape.0 { + for j in self.indptr[i]..self.indptr[i + 1] { + out[[i, self.indices[j]]] += self.data[j]; + } + } + + out + } +} + +fn lu_solve(a: &Array2, b: &Array1) -> Result, Error> { + let a = a.to_nalgebra_matrix().lu(); + let b = b.to_nalgebra_matrix(); + + let x = a + .solve(&b) + .ok_or(Error::MatrixMustBeSquare)? + .to_1d_array()?; + + Ok(x) +} + +impl Mul<&Array1> for CSR { + type Output = Array1; + + fn mul(self, rhs: &Array1) -> Array1 { + let mut out = Array1::::zeros(self.shape.0); + for i in 0..self.shape.0 { + for j in self.indptr[i]..self.indptr[i + 1] { + out[i] += self.data[j] * rhs[self.indices[j]]; + } + } + + out + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_lu_solve() { + let a = Array2::::from(vec![[11., 12., 0.], [0., 22., 23.], [31., 0., 33.]]); + let b = Array1::::from_vec(vec![35., 113., 130.]); + let x_expected = Array1::::from_vec(vec![1., 2., 3.]); + + let x = lu_solve(&a, &b).unwrap(); + + for (x, expected) in x.into_iter().zip(x_expected) { + assert!( + (x - expected).abs() < f64::EPSILON * 10., + "{} {}", + x, + expected + ) + } + } + + #[test] + fn negative_csr_test_00() { + let a = CSR::new( + Array1::::zeros(0), + Array1::::zeros(0), + Array1::::zeros(11), + (1, 3), + ) + .unwrap_err(); + + assert!(matches!(a, Error::CannotInitCSR)); + } +} diff --git a/daemon/src/payout_curve/curve.rs b/daemon/src/payout_curve/curve.rs new file mode 100644 index 0000000..c6112eb --- /dev/null +++ b/daemon/src/payout_curve/curve.rs @@ -0,0 +1,453 @@ +use crate::payout_curve::basis::BSplineBasis; +use crate::payout_curve::splineobject::SplineObject; +use crate::payout_curve::utils::*; +use crate::payout_curve::Error; + +use ndarray::prelude::*; +use ndarray::s; +use std::cmp::max; + +fn default_basis() -> Result, Error> { + let out = vec![BSplineBasis::new(None, None, None)?]; + Ok(out) +} + +#[derive(Clone, Debug)] +pub struct Curve { + pub spline: SplineObject, +} + +/// Represents a curve: an object with a one-dimensional parameter space. +impl Curve { + /// Construct a curve with the given basis and control points. In theory, + /// the curve could be defined in some Euclidean space of dimension N, but + /// for the moment only curves in E^2 are supported via hard-coding. That is, + /// any valid basis set can be passed in when instantiating this object, + /// but only the first one is every considered in the methods provided. + /// + /// The default is to create a linear one-element mapping from (0,1) to the + /// unit interval. + /// + /// ### parameters + /// * basis: The underlying B-Spline basis + /// * controlpoints: An *n* × *d* matrix of control points + /// * rational: Whether the curve is rational (in which case the + /// control points are interpreted as pre-multiplied with the weight, + /// which is the last coordinate) + pub fn new( + bases: Option>, + controlpoints: Option>, + rational: Option, + ) -> Result { + let bases = bases.unwrap_or(default_basis()?); + let spline = SplineObject::new(bases, controlpoints, rational)?; + + Ok(Curve { spline }) + } + + /// Evaluate the object at given parametric values. + /// + /// ### parameters + /// * t: collection of parametric coordinates in which to evaluate + /// Realistically, this should actually be an Array1 object, but for + /// consistency with the underlying SplineObject methods the collection + /// is used instead. + /// + /// ### returns + /// * 2D array + // pub fn evaluate(&self, t: &mut &[Array1], tensor: bool) -> Result, Error> { + pub fn evaluate(&self, t: &mut &[Array1]) -> Result, Error> { + self.spline.validate_domain(t)?; + + let mut tx = t[0].clone().to_owned(); + let n_csr = self.spline.bases[0].evaluate(&mut tx, 0, true)?; + + // kludge... + let n0 = self.spline.controlpoints.shape()[0]; + let n1 = self.spline.controlpoints.shape()[1]; + let flat_controlpoints = self + .ravel(&self.spline.controlpoints) + .into_raw_vec() + .to_vec(); + + let arr = Array2::::from_shape_vec((n0, n1), flat_controlpoints)?; + let mut result = n_csr.todense().dot(&arr); + + // if the spline is rational, we apply the weights and omit the weights column + if self.spline.rational { + let wpos = result.shape()[1] - 1; + let weights = &&result.slice(s![.., wpos]); + let mut temp = Array2::::zeros((result.shape()[0], wpos)); + + for i in 0..self.spline.dimension { + { + let mut slice = temp.slice_mut(s![.., i]); + let update = result.slice(s![.., i]).to_owned() / weights.to_owned(); + slice.assign(&update.view()); + } + } + result = temp; + } + + Ok(result) + } + + /// Extend the curve by merging another curve to the end of it. + /// + /// The curves are glued together in a C0 fashion with enough repeated + /// knots. The function assumes that the end of this curve perfectly + /// matches the start of the input curve. + /// + /// Obviously, neither curve can be periodic, which enables us + /// to assume all knot vectors are open. + /// + /// ### parameters + /// * curve: Another curve + /// + /// ### returns + pub fn append(&mut self, othercurve: Curve) -> Result<(), Error> { + if self.spline.bases[0].periodic > -1 || othercurve.spline.bases[0].periodic > -1 { + return Result::Err(Error::CannotConnectPeriodicCurves); + }; + + let mut extending_curve = othercurve; + + // make sure both are in the same space, and (if needed) have rational weights + self.spline + .make_splines_compatible(&mut extending_curve.spline)?; + let p1 = self.spline.order(0)?[0]; + let p2 = extending_curve.spline.order(0)?[0]; + + if p1 < p2 { + self.raise_order(p2 - p1)?; + } else { + extending_curve.raise_order(p1 - p2)?; + } + + let p = max(p1, p2); + + let old_knot = self.spline.knots(0, Some(true))?[0].clone(); + let mut add_knot = extending_curve.spline.knots(0, Some(true))?[0].clone(); + add_knot -= add_knot[0]; + add_knot += old_knot[old_knot.len() - 1]; + + let mut new_knot = Array1::::zeros(add_knot.len() + old_knot.len() - p - 1); + { + let mut slice = new_knot.slice_mut(s![..old_knot.len() - 1]); + let update = old_knot.slice(s![..old_knot.len() - 1]); + slice.assign(&update.view()); + } + { + let mut slice = new_knot.slice_mut(s![old_knot.len() - 1..]); + let update = add_knot.slice(s![p..]); + slice.assign(&update.view()); + } + + let rational = self.spline.rational as usize; + let n1 = self.spline.controlpoints.shape()[0]; + let n2 = extending_curve.spline.controlpoints.shape()[0]; + let n3 = self.spline.dimension + rational; + let mut new_controlpoints = Array2::::zeros((n1 + n2 - 1, n3)); + { + let mut slice = new_controlpoints.slice_mut(s![..n1, ..]); + let update = self.spline.controlpoints.slice(s![.., ..]); + slice.assign(&update.view()); + } + { + let mut slice = new_controlpoints.slice_mut(s![n1.., ..]); + let update = extending_curve.spline.controlpoints.slice(s![1.., ..]); + slice.assign(&update.view()); + } + + self.spline.bases = vec![BSplineBasis::new(Some(p), Some(new_knot), None)?]; + self.spline.controlpoints = new_controlpoints.into_dyn().to_owned(); + + Ok(()) + } + + /// Raise the polynomial order of the curve. + /// + /// ### parameters + /// * amount: Number of times to raise the order + pub fn raise_order(&mut self, amount: usize) -> Result<(), Error> { + if amount == 0 { + return Ok(()); + } + + // work outside of self, copy back in at the end + let mut new_basis = self.spline.bases[0].clone(); + new_basis.raise_order(amount); + + // set up an interpolation problem. This is in projective space, + // so no problems for rational cases + let mut interpolation_pts_t = new_basis.greville(); + let n_old = self.spline.bases[0].evaluate(&mut interpolation_pts_t, 0, true)?; + let n_new = new_basis.evaluate(&mut interpolation_pts_t, 0, true)?; + + // Some kludge required to enable .dot(), which doesn't work on dynamic + // arrays. Surely a better way to do this, but this is quick and dirty + // and valid since we're in curve land + let n0 = self.spline.controlpoints.shape()[0]; + let n1 = self.spline.controlpoints.shape()[1]; + let flat_controlpoints = self + .ravel(&self.spline.controlpoints) + .into_raw_vec() + .to_vec(); + + let arr = Array2::::from_shape_vec((n0, n1), flat_controlpoints)?; + let interpolation_pts_x = n_old.todense().dot(&arr); + let res = n_new.matrix_solve(&interpolation_pts_x)?; + + self.spline.controlpoints = res.into_dyn().to_owned(); + self.spline.bases = vec![new_basis]; + + Ok(()) + } + + /// Computes the euclidian length of the curve in geometric space + /// + /// .. math:: \\int_{t_0}^{t_1}\\sqrt{x(t)^2 + y(t)^2 + z(t)^2} dt + /// + /// ### parameters + /// * t0: lower integration limit + /// * t1: upper integration limit + pub fn length(&self, t0: Option, t1: Option) -> Result { + let mut knots = &self.spline.knots(0, Some(false))?[0]; + + // keep only integration boundaries within given start (t0) and stop (t1) interval + let new_knots_0 = t0 + .map(|t0| { + let i = bisect_left(knots, &t0, knots.len()); + let mut vec = Vec::::with_capacity(&knots.to_vec()[i..].len() + 1); + vec.push(t0); + for elem in knots.to_vec()[i..].iter() { + vec.push(*elem); + } + Array1::::from_vec(vec) + }) + .unwrap_or_else(|| knots.to_owned()); + knots = &new_knots_0; + + let new_knots_1 = t1 + .map(|t1| { + let i = bisect_right(knots, &t1, knots.len()); + let mut vec = Vec::::with_capacity(&knots.to_vec()[..i].len() + 1); + for elem in knots.to_vec()[..i].iter() { + vec.push(*elem); + } + vec.push(t1); + Array1::::from_vec(vec) + }) + .unwrap_or_else(|| knots.to_owned()); + knots = &new_knots_1; + + let klen = knots.len(); + let gleg = GaussLegendreQuadrature::new(self.spline.order(0)?[0] + 1)?; + + let t = &knots.to_vec()[..klen - 1] + .iter() + .zip(knots.to_vec()[1..].iter()) + .map(|(t0, t1)| (&gleg.sample_points + 1.) / 2. * (t1 - t0) + *t0) + .collect::>(); + + let w = &knots.to_vec()[..klen - 1] + .iter() + .zip(knots.to_vec()[1..].iter()) + .map(|(t0, t1)| &gleg.weights / 2. * (t1 - t0)) + .collect::>(); + + let t_flat = self.flattened(&t[..]); + let w_flat = self.flattened(&w[..]); + + let dx = self.derivative(&t_flat, 1, Some(true))?; + let det_j = Array1::::from_vec( + dx.mapv(|e| e.powi(2)) + .sum_axis(Axis(1)) + .mapv(f64::sqrt) + .iter() + .copied() + .collect::>(), + ); + let out = det_j.dot(&w_flat); + + Ok(out) + } + + fn flattened(&self, vec_arr: &[Array1]) -> Array1 { + let alloc = vec_arr.iter().fold(0, |sum, e| sum + e.len()); + let mut vec_out = Vec::::with_capacity(alloc); + for arr in vec_arr.iter() { + for e in arr.to_vec().iter() { + vec_out.push(*e); + } + } + + Array1::::from_vec(vec_out) + } + + /// left here as a private method as it assumes C-contiguous ordering, + /// which is fine for where we use it here. + fn ravel(&self, arr: &ArrayD) -> Array1 { + let alloc = arr.shape().iter().product(); + let mut vec = Vec::::with_capacity(alloc); + for e in arr.iter() { + vec.push(*e) + } + + Array1::::from_vec(vec) + } + + /// Evaluate the derivative of the curve at the given parametric values. + /// + /// This function returns an *n* × *dim* array, where *n* is the number of + /// evaluation points, and *dim* is the physical dimension of the curve. + /// **At this point in time, only `dim == 1` works, owing to the provisional + /// constraints on the struct `Curve` itself.** + /// + /// ### parameters + /// * t: Parametric coordinates in which to evaluate + /// * d: Number of derivatives to compute + /// * from_right: Evaluation in the limit from right side + pub fn derivative( + &self, + t: &Array1, + d: usize, + from_right: Option, + ) -> Result, Error> { + let from_right = from_right.unwrap_or(true); + + if !self.spline.rational || d < 2 || d > 3 { + let mut tx = &vec![t.clone().to_owned()][..]; + let res = self.spline.derivative(&mut tx, &[d], &[from_right], true)?; + return Ok(res); + } + + // init rusult array + let mut res = Array2::::zeros((t.len(), self.spline.dimension)); + + // annoying fix to make the controlpoints not dynamic--implicit + // assumption of 2D curve only! + let n0 = self.spline.controlpoints.shape()[0]; + let n1 = self.spline.controlpoints.shape()[1]; + let flat_controlpoints = self + .ravel(&self.spline.controlpoints) + .into_raw_vec() + .to_vec(); + let static_controlpoints = Array2::::from_shape_vec((n0, n1), flat_controlpoints)?; + let mut t_eval = t.clone(); + let d2 = self.spline.bases[0] + .evaluate(&mut t_eval, 2, from_right)? + .todense() + .dot(&static_controlpoints); + let d1 = self.spline.bases[0] + .evaluate(&mut t_eval, 1, from_right)? + .todense() + .dot(&static_controlpoints); + let d0 = self.spline.bases[0] + .evaluate(&mut t_eval, 0, true)? + .todense() + .dot(&static_controlpoints); + let w0 = &d0.slice(s![.., d0.shape()[1] - 1]).to_owned(); + let w1 = &d1.slice(s![.., d1.shape()[1] - 1]).to_owned(); + let w2 = &d2.slice(s![.., d2.shape()[1] - 1]).to_owned(); + + if d == 2 { + let w0_cube = &w0.mapv(|e| e.powi(3)).to_owned(); + for i in 0..self.spline.dimension { + { + let update = &((d2.slice(s![.., i]).to_owned() * w0 * w0 + - 2. * w1 + * (d1.slice(s![.., i]).to_owned() * w0 + - d0.slice(s![.., i]).to_owned() * w1) + - d0.slice(s![.., i]).to_owned() * w2 * w1) + / w0_cube); + + let mut slice = res.slice_mut(s![.., i]); + slice.assign(update); + } + } + } + + if d == 3 { + let d3 = self.spline.bases[0] + .evaluate(&mut t_eval, 3, from_right)? + .todense() + .dot(&static_controlpoints); + let w3 = &d3.slice(s![.., d3.shape()[1] - 1]); + let w0_four = w0.mapv(|e| e.powi(6)); + for i in 0..self.spline.dimension { + { + let h0 = &(d1.slice(s![.., i]).to_owned() * w0 + - d0.slice(s![.., i]).to_owned() * w1); + let h1 = &(d2.slice(s![.., i]).to_owned() * w0 + - d0.slice(s![.., i]).to_owned() * w2); + let h2 = &(d3.slice(s![.., i]).to_owned() * w0 + + d2.slice(s![.., i]).to_owned() * w1 + - d1.slice(s![.., i]).to_owned() * w2 + - d0.slice(s![.., i]).to_owned() * w3); + let g0 = &(h1 * w0 - 2. * h0 * w1); + let g1 = &(h2 * w0 - 2. * h0 * w2 - h1 * w1); + + let update = (g1 * w0 - 3. * g0 * w1) / &w0_four; + let mut slice = res.slice_mut(s![.., i]); + slice.assign(&update); + } + } + } + + Ok(res.into_dyn().to_owned()) + } + + /// Computes the L2 (squared and per knot span) between this + /// curve and a target curve as well as the L_infinity error: + /// + /// .. math:: ||\\boldsymbol{x_h}(t)-\\boldsymbol{x}(t)||_{L^2(t_1,t_2)}^2 = \\int_{t_1}^{t_2} + /// |\\boldsymbol{x_h}(t)-\\boldsymbol{x}(t)|^2 dt, \\quad \\forall \\;\\text{knots}\\;t_1 < + /// t_2 + /// + /// .. math:: ||\\boldsymbol{x_h}(t)-\\boldsymbol{x}(t)||_{L^\\infty} = \\max_t + /// |\\boldsymbol{x_h}(t)-\\boldsymbol{x}(t)| + /// + /// ### parameters + /// * function target: callable function which takes as input a vector + /// of evaluation points t and gives as output a matrix x where + /// x[i,j] is component j evaluated at point t[i] + /// + /// ### returns + /// * L2 error per knot-span + pub fn error( + &self, + target: impl Fn(&Array1) -> Array2, + ) -> Result<(Array1, f64), Error> { + let knots = &self.spline.knots(0, Some(false))?[0]; + let n = self.spline.order(0)?[0]; + let gleg = GaussLegendreQuadrature::new(n + 1)?; + + let mut error_l2 = Vec::with_capacity(knots.len() - 1); + let mut error_linf = Vec::with_capacity(knots.len() - 1); + + for (t0, t1) in knots.to_vec()[..knots.len() - 1] + .iter() + .zip(&mut knots.to_vec()[1..].iter()) + { + let tg = (&gleg.sample_points + 1.) / 2. * (t1 - t0) + *t0; + let eval = vec![tg.clone()]; + let wg = &gleg.weights / 2. * (t1 - t0); + + let exact = target(&tg); + let error = self.evaluate(&mut &eval[..])? - exact; + let error_2 = &error.mapv(|e| e.powi(2)).sum_axis(Axis(1)); + let error_abs = &error_2.mapv(|e| e.sqrt()); + + let l2_val = error_2.dot(&wg); + let linf_val = error_abs.iter().copied().fold(f64::NEG_INFINITY, f64::max); + + error_l2.push(l2_val); + error_linf.push(linf_val); + } + + let out_inf = error_linf.iter().copied().fold(f64::NEG_INFINITY, f64::max); + + Ok((Array1::::from_vec(error_l2), out_inf)) + } +} diff --git a/daemon/src/payout_curve/curve_factory.rs b/daemon/src/payout_curve/curve_factory.rs new file mode 100644 index 0000000..fa57221 --- /dev/null +++ b/daemon/src/payout_curve/curve_factory.rs @@ -0,0 +1,159 @@ +use crate::payout_curve::basis::BSplineBasis; +use crate::payout_curve::curve::Curve; +use crate::payout_curve::utils::cmp_f64; +use crate::payout_curve::Error; + +use ndarray::prelude::*; + +/// Perform general spline interpolation on a provided basis. +/// +/// ### parameters +/// * x: Matrix *X\[i,j\]* of interpolation points *x_i* with components *j* +/// * basis: Basis on which to interpolate +/// * t: parametric values at interpolation points; defaults to +/// Greville points if not provided +/// +/// ### returns +/// * Interpolated curve +pub fn interpolate( + x: &Array2, + basis: &BSplineBasis, + t: Option>, +) -> Result { + let mut t = t.unwrap_or_else(|| basis.greville()); + let evals = basis.evaluate(&mut t, 0, true)?; + let controlpoints = evals.matrix_solve(x)?; + let out = Curve::new(Some(vec![basis.clone()]), Some(controlpoints), None)?; + + Ok(out) +} + +/// Computes an interpolation for a parametric curve up to a specified +/// tolerance. The method will iteratively refine parts where needed +/// resulting in a non-uniform knot vector with as optimized knot +/// locations as possible. +/// +/// ### parameters +/// * x: callable function `x: t --> (t, x(t))` which takes as input a vector +/// of evaluation points `t` and gives as output a matrix `x` where `x\[i,j\]` +/// is component `j` evaluated at point `t\[i\]` +/// * t0: start of parametric domain +/// * t1: end of parametric domain +/// * rtol: relative tolerance for stopping criterium. It is defined to be +/// `||e||_L2 / D`, where `D` is the length of the curve and `||e||_L2` is +/// the L2-error (see Curve.error) +/// * atol: absolute tolerance for stopping criterium. It is defined to be +/// the maximal distance between the curve approximation and the exact curve +/// +/// ### returns +/// Curve (NURBS) +pub fn fit( + x: impl Fn(&Array1) -> Array2, + t0: f64, + t1: f64, + rtol: Option, + atol: Option, +) -> Result { + let rtol = rtol.unwrap_or(1e-4); + let atol = atol.unwrap_or(0.0); + + let knot_vector = Array1::::from_vec(vec![t0, t0, t0, t0, t1, t1, t1, t1]); + let b = BSplineBasis::new(Some(4), Some(knot_vector), None)?; + let t = b.greville(); + let exact = &x(&t); + + let mut crv = interpolate(exact, &b, Some(t))?; + let err = crv.error(&x)?; + + // polynomial input (which can be exactly represented) only use one knot span + if err.1 < 1e-13 { + return Ok(crv); + } + + // for all other curves, start with 4 knot spans + let mut knot_vec = Vec::::with_capacity(12); + for _ in 0..4 { + knot_vec.push(t0) + } + for i in 0..4 { + let i_64 = (i + 1) as f64; + let val = i_64 / 5. * (t1 - t0) + t0; + knot_vec.push(val); + } + for _ in 0..4 { + knot_vec.push(t1) + } + let knot_vector = Array1::::from_vec(knot_vec.clone()); + let b = BSplineBasis::new(Some(4), Some(knot_vector), None)?; + let t = b.greville(); + let exact = &x(&t); + + crv = interpolate(exact, &b, Some(t))?; + let err = crv.error(&x)?; + let mut err_l2 = err.0; + let mut err_max = err.1; + + // this is technically false since we need the length of the target function *x* + // and not our approximation *crv*, but we don't have the derivative of *x*, so + // we can't compute it. This seems like a healthy compromise + let length = crv.length(None, None)?; + let mut target = (err_l2.sum() / length).sqrt(); + + // conv_order = 4 + // square_conv_order = 2 * conv_order + // scale = square_conv_order + 4 + let scale_64 = 12_f64; + + while target > rtol && err_max > atol { + let knot_span = &crv.spline.knots(0, Some(false))?[0]; + let target_error = (rtol * length).powi(2) / err_l2.len() as f64; + for i in 0..err_l2.len() { + // figure out how many new knots we require in this knot interval: + // if we converge with *scale* and want an error of *target_error* + // |e|^2 * (1/n)^scale = target_error^2 + let n = ((err_l2[i].ln() - target_error.ln()) / scale_64) + .exp() + .ceil() as usize; + + // add *n* new interior knots to this knot span + let new_knots = Array1::::linspace(knot_span[i], knot_span[i + 1], n + 1); + for e in new_knots.slice(s![1..new_knots.len() - 1]).iter() { + knot_vec.push(*e); + } + } + + // build new refined knot vector + knot_vec.sort_by(cmp_f64); + let knot_vector = Array1::::from_vec(knot_vec.clone()); + let b = BSplineBasis::new(Some(4), Some(knot_vector), None)?; + + // do interpolation and return result + let t = b.greville(); + let exact = &x(&t); + + crv = interpolate(exact, &b, Some(t))?; + let err = crv.error(&x)?; + err_l2 = err.0; + err_max = err.1; + target = err_l2.sum().sqrt() / length; + } + + Ok(crv) +} + +/// Create a line between two points. +/// +/// ### parameters +/// * a, b: start and end points (resp.) +/// * relative: whether `b` is relative to `a` or absolute +pub fn line(a: (f64, f64), b: (f64, f64), relative: bool) -> Result { + let vec; + if relative { + vec = vec![[a.0, a.1], [a.0 + b.0, a.1 + b.1]]; + } else { + vec = vec![[a.0, a.1], [b.0, b.1]]; + } + let controlpoints = Array2::::from(vec); + + Curve::new(None, Some(controlpoints), None) +} diff --git a/daemon/src/payout_curve/splineobject.rs b/daemon/src/payout_curve/splineobject.rs new file mode 100644 index 0000000..1b5507a --- /dev/null +++ b/daemon/src/payout_curve/splineobject.rs @@ -0,0 +1,494 @@ +use crate::payout_curve::basis::BSplineBasis; +use crate::payout_curve::csr_tools::CSR; +use crate::payout_curve::Error; +use itertools::Itertools; +use ndarray::prelude::*; +use ndarray::{concatenate, Order}; +use ndarray_einsum_beta::{einsum, tensordot}; +use std::collections::HashMap; + +#[derive(Clone, Debug)] +pub struct SplineObject { + pub bases: Vec, + pub controlpoints: ArrayD, + pub dimension: usize, + pub rational: bool, + pub pardim: usize, +} + +/// Master struct for spline objects with arbitrary dimensions. +/// +/// This class should be composed instead of used directly. +impl SplineObject { + pub fn new( + bases: Vec, + controlpoints: Option>, + rational: Option, + ) -> Result { + let mut controlpoints = match controlpoints { + Some(controlpoints) => controlpoints, + None => default_control_points(&bases)?, + }; + let rational = rational.unwrap_or(false); + + if controlpoints.slice(s![0, ..]).shape()[0] == 1 { + controlpoints = concatenate( + Axis(1), + &[ + controlpoints.view(), + Array1::::zeros(controlpoints.shape()[0]) + .insert_axis(Axis(1)) + .view(), + ], + )?; + } + + if rational { + controlpoints = concatenate( + Axis(1), + &[ + controlpoints.view(), + Array1::::ones(controlpoints.shape()[0]) + .insert_axis(Axis(1)) + .view(), + ], + )?; + } + + let dim = controlpoints.shape()[1] - (rational as usize); + let bases_shape = determine_shape(&bases)?; + let ncomps = dim + (rational as usize); + let cpts_shaped = reshaper(controlpoints, bases_shape, ncomps)?; + let pardim = cpts_shaped.shape().len() - 1; + + Ok(SplineObject { + bases, + controlpoints: cpts_shaped, + dimension: dim, + rational, + pardim, + }) + } + + /// Check whether the given evaluation parameters are valid + pub fn validate_domain(&self, t: &[Array1]) -> Result<(), Error> { + for (basis, params) in self.bases.iter().zip(t.to_owned().iter_mut()) { + if basis.periodic < 0 { + basis.snap(&mut *params); + let p_max = ¶ms.iter().copied().fold(f64::NEG_INFINITY, f64::max); + let p_min = ¶ms.iter().copied().fold(f64::INFINITY, f64::min); + if *p_min < basis.start() || basis.end() < *p_max { + return Result::Err(Error::InvalidDomain); + } + } + } + + Ok(()) + } + + fn tensor_evaluate(&self, eval_bases: &mut [CSR], tensor: bool) -> Result, Error> { + // KLUDGE! + // owing to the fact that the conventional ellipsis notation is not yet + // implemented for einsum, we use this workaround that should cover us. + // If not, just grow the maps as needed or address the issue: + // https://github.com/oracleofnj/einsum/issues/6 + let init_map: HashMap = [ + (2, "ij,jp->ip"), + (3, "ij,jpq->ipq"), + (4, "ij,jpqr->ipqr"), + (5, "ij,jpqrs->ipqrs"), + (6, "ij,jpqrst->ipqrst"), + ] + .iter() + .cloned() + .collect(); + + let iter_map: HashMap = [ + (3, "ij,ijp->ip"), + (4, "ij,ijpq->ipq"), + (5, "ij,ijpqr->ipqr"), + (6, "ij,ijpqrs->ipqrs"), + ] + .iter() + .cloned() + .collect(); + + let mut out; + if tensor { + eval_bases.reverse(); + let cpts = self.controlpoints.clone().to_owned(); + let idx = eval_bases.len() - 1; + + out = eval_bases.iter().fold(cpts, |e, tns| { + tensordot(&tns.todense(), &e, &[Axis(1)], &[Axis(idx)]) + }); + } else { + let mut pos = 0; + let mut key = self.bases.len() + 1; + let mut val = match init_map.get(&key) { + Some(val) => Ok(val), + _ => Result::Err(Error::NoEinsumOperatorString), + }?; + + out = einsum(val, &[&eval_bases[pos].todense(), &self.controlpoints]) + .map_err(|_| Error::Einsum)?; + + for _ in eval_bases.iter().skip(1) { + pos += 1; + val = match iter_map.get(&key) { + Some(val) => Ok(val), + _ => Result::Err(Error::NoEinsumOperatorString), + }?; + let temp = out.clone().to_owned(); + + out = + einsum(val, &[&eval_bases[pos].todense(), &temp]).map_err(|_| Error::Einsum)?; + + key -= 1; + } + } + // *** END KLUDGE **** + + Ok(out) + } + + /// Evaluate the derivative of the object at the given parametric values. + /// + /// If *tensor* is true, evaluation will take place on a tensor product + /// grid, i.e. it will return an *n1* × *n2* × ... × *dim* array, where + /// *ni* is the number of evaluation points in direction *i*, and *dim* is + /// the physical dimension of the object. + /// + /// If *tensor* is false, there must be an equal number *n* of evaluation + /// points in all directions, and the return value will be an *n* × *dim* + /// array. + /// + /// ### parameters + /// * t = [u,v,...]: Parametric coordinates in which to evaluate + /// * d: Order of derivative to compute, index corresponds to bases index + /// * from_right: Evaluation in the limit from above; index orresponds to bases index + /// * tensor: Whether to evaluate on a tensor product grid + pub fn derivative( + &self, + t: &mut &[Array1], + d: &[usize], + from_right: &[bool], + tensor: bool, + ) -> Result, Error> { + // check + let testlen = t.len(); + let ops = [t.len(), d.len(), from_right.len()] + .iter() + .all(|e| *e == testlen); + if !tensor && !ops { + return Result::Err(Error::InvalidDerivative); + } + + self.validate_domain(t)?; + + // Evaluate the derivatives of the corresponding bases at the corresponding points + // and build the result array + let mut evals = self + .bases + .iter() + .zip(t.iter().zip(d.iter().zip(from_right.iter()))) + .map(|(b, (t, (d, r)))| { + let mut tx = t.clone(); + let eval = b.evaluate(&mut tx, *d, *r)?; + Ok(eval) + }) + .collect::, Error>>()?; + let mut result = self.tensor_evaluate(&mut evals, tensor)?; + + // For rational curves, we need to use the quotient rule + // (n/W)' = (n' W - n W') / W^2 = n'/W - nW'/W^2 + // n'(i) = result[..., i] + // W'(i) = result[..., -1] + if self.rational { + if d.iter().sum::() > 1 { + return Result::Err(Error::DerivativeNotImplemented); + } + + let mut ns = self + .bases + .iter() + .zip(t.iter()) + .map(|(b, t)| { + let mut tx = t.clone(); + let eval = b.evaluate(&mut tx, 0, true)?; + Ok(eval) + }) + .collect::, Error>>()?; + let non_derivative = self.tensor_evaluate(&mut ns, tensor)?; + + let axis_w = non_derivative.shape().len() - 1; + let idx_w = non_derivative.shape()[axis_w] - 1; + let w = &non_derivative.index_axis(Axis(axis_w), idx_w).to_owned(); + let w_square = &w.mapv(|e| e.powi(2)).to_owned(); + + let axis_r = result.shape().len() - 1; + let idx_r = result.shape()[axis_r] - 1; + let wd = &result.index_axis(Axis(axis_r), idx_r).to_owned(); + + for i in 0..self.dimension { + { + let update = &(result.index_axis(Axis(axis_r), i).to_owned() / w + - non_derivative.index_axis(Axis(axis_w), i).to_owned() * wd / w_square); + let mut slice = result.index_axis_mut(Axis(axis_r), i); + slice.assign(update); + } + } + + // delete the last column; some faffing about required to maintain + // C-contiguous ordering. Probably a much better way to do this... + let res_shape = &result.shape().iter().copied().collect::>(); + let mut n_res: usize = res_shape[..axis_r].iter().product(); + n_res *= idx_r; + let idx = (0..axis_r).collect::>(); + + // this ends up being F-contiguous, every time + let res_slice = result.select(Axis(axis_r), &idx[..]).to_owned(); + let raveled = res_slice.to_shape(((n_res,), Order::C)).unwrap(); + let fixed = raveled.to_shape((res_slice.shape(), Order::C))?.to_owned(); + + result = fixed; + } + + Ok(result) + } + + /// Return knots vector + /// + /// If `direction` is given, returns the knots in that direction only. + /// Otherwise, specifying direction as a negative value returns the + /// knots of all directions. + /// + /// ### parameters + /// * direction: Direction number (axis) in which to get the knots. + /// * with_multiplicities: If true, return knots with multiplicities \ + /// (i.e. repeated). + pub fn knots( + &self, + direction: isize, + with_multiplicities: Option, + ) -> Result>, Error> { + let with_multiplicities = with_multiplicities.unwrap_or(false); + let out; + + if direction < 0 { + if with_multiplicities { + out = self + .bases + .iter() + .map(|e| e.knots.clone().to_owned()) + .collect::>(); + } else { + out = self + .bases + .iter() + .map(|e| e.knot_spans(false).to_owned()) + .collect::>(); + } + } else { + let p = direction as usize; + if with_multiplicities { + out = (&[self.bases[p].knots.clone().to_owned()]).to_vec(); + } else { + out = (&[self.bases[p].knot_spans(false).to_owned()]).to_vec(); + } + } + + Ok(out) + } + + /// This will manipulate one or both to ensure that they are both rational + /// or nonrational, and that they lie in the same physical space. + pub fn make_splines_compatible(&mut self, otherspline: &mut SplineObject) -> Result<(), Error> { + if self.rational { + otherspline.force_rational()?; + } else if otherspline.rational { + self.force_rational()?; + } + + if self.dimension > otherspline.dimension { + otherspline.set_dimension(self.dimension)?; + } else { + self.set_dimension(otherspline.dimension)?; + } + + Ok(()) + } + + /// Force a rational representation of the object. + pub fn force_rational(&mut self) -> Result<(), Error> { + if !self.rational { + self.controlpoints = self.insert_phys(&self.controlpoints, 1f64)?; + self.rational = true; + } + + Ok(()) + } + + /// Sets the physical dimension of the object. If increased, the new + /// components are set to zero. + /// + /// ### parameters + /// * new_dim: New dimension + pub fn set_dimension(&mut self, new_dim: usize) -> Result<(), Error> { + let mut dim = self.dimension; + + while new_dim > dim { + self.controlpoints = self.insert_phys(&self.controlpoints, 0f64)?; + dim += 1; + } + + while new_dim < dim { + let axis = if self.rational { -2 } else { -1 }; + self.controlpoints = self.delete_phys(&self.controlpoints, axis)?; + dim -= 1; + } + + self.dimension = new_dim; + + Ok(()) + } + + fn insert_phys(&self, arr: &ArrayD, insert_value: f64) -> Result, Error> { + let mut arr_shape = arr.shape().to_vec(); + let n = arr_shape[arr_shape.len() - 1]; + let arr_prod = arr_shape.iter().product(); + + let raveled = arr.to_shape(((arr_prod,), Order::C))?; + let mut new_arr = Array1::::zeros(0); + + for i in (0..raveled.len()).step_by(n) { + let new_row = concatenate( + Axis(0), + &[ + raveled.slice(s![i..i + n]).view(), + (insert_value * Array1::::ones(1)).view(), + ], + )?; + new_arr = concatenate(Axis(0), &[new_arr.view(), new_row.view()])?; + } + + arr_shape[n] += 1; + let out = new_arr.to_shape((&arr_shape[..], Order::C))?.to_owned(); + + Ok(out) + } + + fn delete_phys(&self, arr: &ArrayD, axis: isize) -> Result, Error> { + let mut arr_shape = arr.shape().to_vec(); + let n = arr_shape[arr_shape.len() - 1]; + let step = (n as isize + axis) as usize; + let arr_prod = arr_shape.iter().product(); + + let raveled = arr.to_shape(((arr_prod,), Order::C))?; + let mut new_arr = Array1::::zeros(0); + + for i in (0..raveled.len()).step_by(n) { + let new_row; + if axis < -1 { + let front = raveled.slice(s![i..i + step]).clone().to_owned(); + let tail = raveled.slice(s![i + step + 1..i + n]).clone().to_owned(); + new_row = concatenate(Axis(0), &[front.view(), tail.view()])?; + } else { + new_row = raveled.slice(s![i..i + step]).clone().to_owned(); + } + new_arr = concatenate(Axis(0), &[new_arr.view(), new_row.view()])?; + } + + arr_shape[n - 1] -= 1; + let out = new_arr.to_shape((&arr_shape[..], Order::C))?.to_owned(); + + Ok(out) + } + + /// Return polynomial order (degree + 1). + /// + /// If `direction` is given, returns the order in that direction only. + /// Otherwise, specifying direction as a negative value returns the + /// order of all directions. + /// + /// ### parameters + /// * direction: Direction in which to get the order. + pub fn order(&self, direction: isize) -> Result, Error> { + let out; + if direction < 0 { + out = self.bases.iter().map(|e| e.order).collect::>(); + } else { + let p = direction as usize; + out = (&[self.bases[p].order]).to_vec(); + } + + Ok(out) + } +} + +fn default_control_points(bases: &[BSplineBasis]) -> Result, Error> { + let mut temp = bases + .iter() + .rev() + .map(|b| { + let mut v = b.greville().into_raw_vec(); + v.reverse(); + v + }) + .multi_cartesian_product() + .collect::>(); + temp.reverse(); + + // because the above is just a little bit incorrect... + for elem in temp.iter_mut() { + elem.reverse(); + } + + let mut data = Vec::new(); + let ncols = temp.first().map_or(0, |row| row.len()); + let mut nrows = 0; + + for elem in temp.iter() { + data.extend_from_slice(elem); + nrows += 1; + } + + let out = Array2::from_shape_vec((nrows, ncols), data)?; + + Ok(out) +} + +/// Custom reshaping function to preserve control points of several +/// dimensions that are stored contiguously. +/// +/// The return value has shape (*newshape, ncomps), where ncomps is +/// the number of components per control point, as inferred by the +/// size of `arr` and the desired shape. +fn reshaper( + arr: Array2, + mut newshape: Vec, + ncomps: usize, +) -> Result, Error> { + newshape.reverse(); + newshape.push(ncomps); + + let mut spec: Vec = (0..newshape.len() - 1).collect(); + spec.reverse(); + spec.push(newshape.len() - 1); + + let tmp = arr.to_shape((&newshape[..], Order::C))?; + let tmp = tmp.to_owned().into_dyn(); + + let out = tmp.view().permuted_axes(&spec[..]).to_owned(); + + Ok(out) +} + +fn determine_shape(bases: &[BSplineBasis]) -> Result, Error> { + let out = bases + .iter() + .map(|e| e.num_functions()) + .collect::>(); + + Ok(out) +} diff --git a/daemon/src/payout_curve/utils.rs b/daemon/src/payout_curve/utils.rs new file mode 100644 index 0000000..aadbdad --- /dev/null +++ b/daemon/src/payout_curve/utils.rs @@ -0,0 +1,168 @@ +use crate::payout_curve::Error; +use ndarray::prelude::*; +use ndarray::s; +use std::cmp::Ordering; +use std::f64::consts::PI; + +pub fn bisect_left(arr: &Array1, val: &f64, mut hi: usize) -> usize { + let mut lo: usize = 0; + while lo < hi { + let mid = (lo + hi) / 2; + if arr[mid] < *val { + lo = mid + 1; + } else { + hi = mid; + } + } + + lo +} + +pub fn bisect_right(arr: &Array1, val: &f64, mut hi: usize) -> usize { + let mut lo: usize = 0; + while lo < hi { + let mid = (lo + hi) / 2; + if *val < arr[mid] { + hi = mid; + } else { + lo = mid + 1; + } + } + + lo +} + +pub fn cmp_f64(a: &f64, b: &f64) -> Ordering { + if a < b { + return Ordering::Less; + } else if a > b { + return Ordering::Greater; + } + Ordering::Equal +} + +/// Gauss-Legendre_Quadrature +/// +/// Could not find a rust implementation of this, so have created one from +/// a C implementation found +/// [here](https://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature#C). +/// +/// The code is well short of optimal, but it gets things moving. Better +/// versions are provided by, for example, +/// [numpy.polynomial.legendre.leggauss](https://github.com/numpy/numpy/blob/v1.21.0/numpy/polynomial/legendre.py#L1519-L1584) +/// but the implementaitons are more involved so we have opted for quick and +/// dirty for the time being. +#[derive(Debug, Clone)] +pub struct GaussLegendreQuadrature { + pub sample_points: Array1, + pub weights: Array1, +} + +impl GaussLegendreQuadrature { + pub fn new(order: usize) -> Result { + if order < 1 { + return Result::Err(Error::DegreeMustBePositive); + } + + let data = legendre_wrapper(&order); + + Ok(GaussLegendreQuadrature { + sample_points: data.0, + weights: data.1, + }) + } +} + +fn legendre_wrapper(order: &usize) -> (Array1, Array1) { + let arr = legendre_coefficients(order); + legendre_roots(&arr) +} + +fn legendre_coefficients(order: &usize) -> Array2 { + let mut lcoef_arr = Array2::::zeros((*order + 1, *order + 1)); + lcoef_arr[[0, 0]] = 1.; + lcoef_arr[[1, 1]] = 1.; + + for n in 2..*order + 1 { + let n_64 = n as f64; + lcoef_arr[[n, 0]] = -(n_64 - 1.) * lcoef_arr[[n - 2, 0]] / n_64; + + for i in 1..n + 1 { + lcoef_arr[[n, i]] = ((2. * n_64 - 1.) * lcoef_arr[[n - 1, i - 1]] + - (n_64 - 1.) * lcoef_arr[[n - 2, i]]) + / n_64; + } + } + + lcoef_arr +} + +fn legendre_eval(coeff_arr: &Array2, n: &usize, x: &f64) -> f64 { + let mut s = coeff_arr[[*n, *n]]; + for i in (1..*n + 1).rev() { + s = s * (*x) + coeff_arr[[*n, i - 1]]; + } + + s +} + +fn legendre_diff(coeff_arr: &Array2, n: &usize, x: &f64) -> f64 { + let n_64 = *n as f64; + n_64 * (x * legendre_eval(coeff_arr, n, x) - legendre_eval(coeff_arr, &(n - 1), x)) + / (x * x - 1.) +} + +fn legendre_roots(coeff_arr: &Array2) -> (Array1, Array1) { + let n = coeff_arr.shape()[0] - 1; + let n_64 = n as f64; + + let mut sample_points_arr = Array1::::zeros(n + 1); + let mut weights_arr = Array1::::zeros(n + 1); + + for i in 1..n + 1 { + let i_64 = i as f64; + let mut x = (PI * (i_64 - 0.25) / (n_64 + 0.5)).cos(); + let mut x1 = x; + x -= legendre_eval(coeff_arr, &n, &x) / legendre_diff(coeff_arr, &n, &x); + + while fdim(&x, &x1) > 2e-16 { + x1 = x; + x -= legendre_eval(coeff_arr, &n, &x) / legendre_diff(coeff_arr, &n, &x); + } + + sample_points_arr[i - 1] = x; + x1 = legendre_diff(coeff_arr, &n, &x); + weights_arr[i - 1] = 2. / ((1. - x * x) * x1 * x1); + } + + // truncate the dummy value off the end + reverse sample points + + // use symmetry to stable things up a bit. + let mut samples = sample_points_arr.slice(s![..n; -1]).to_owned(); + samples = symmetric_samples(&samples); + + let mut weights = weights_arr.slice(s![..n]).to_owned(); + weights = symmetric_weights(&weights); + + (samples, weights) +} + +fn symmetric_samples(arr: &Array1) -> Array1 { + let arr_rev = arr.slice(s![..; -1]).to_owned(); + (arr - &arr_rev) / 2. +} + +fn symmetric_weights(arr: &Array1) -> Array1 { + let s = &arr.sum_axis(Axis(0)); + let arr_rev = arr.slice(s![..; -1]).to_owned(); + (arr + &arr_rev) / s +} + +fn fdim(a: &f64, b: &f64) -> f64 { + let res; + if a - b > 0f64 { + res = a - b; + } else { + res = 0f64; + } + res +} diff --git a/daemon/src/setup_contract.rs b/daemon/src/setup_contract.rs index e3e43d7..c8274ea 100644 --- a/daemon/src/setup_contract.rs +++ b/daemon/src/setup_contract.rs @@ -75,12 +75,7 @@ pub async fn new( let payouts = HashMap::from_iter([( announcement.clone(), - payout_curve::calculate( - cfd.order.price, - cfd.quantity_usd, - params.maker().lock_amount, - (params.taker().lock_amount, cfd.order.leverage), - )?, + payout_curve::calculate(cfd.order.price, cfd.quantity_usd, cfd.order.leverage)?, )]); let own_cfd_txs = create_cfd_transactions( @@ -280,12 +275,7 @@ pub async fn roll_over( id: announcement.id.0, nonce_pks: announcement.nonce_pks.clone(), }, - payout_curve::calculate( - cfd.order.price, - cfd.quantity_usd, - maker_lock_amount, - (taker_lock_amount, cfd.order.leverage), - )?, + payout_curve::calculate(cfd.order.price, cfd.quantity_usd, cfd.order.leverage)?, )]); // unsign lock tx because PartiallySignedTransaction needs an unsigned tx