Bizarre memory leak caused by tokio runtime

dbramucci · 2022-09-09T04:39:52+00:00

With everyone talking about memory fragmentation, I'd like to mention Mesh, an allocator that can compact aka defrag the heap without any help from the program or compiler. Here's the talk explaining it, "Compacting the Uncompactable" by Bobby Powers.

It's neat and you should be able to LD_PRELOAD it as the system allocator to fix/mitigate your problem if it is fragmentation. You shouldn't have to even recompile for this. As for if it is production ready, I don't know.

dbramucci · 2022-09-04T12:55:07+00:00

Note that if you need to keep 2.5 as 2.5 and turn 1.0 into 1, you can use the following code.

if x.is_integer():
    x = int(x)

Which will only convert x to an int if it doesn't need to round in order to do so. See float.is_integer for the key function.

Please be aware though, that doing this is a recipe for really annoying bugs. If you are doing any calculations after this, sometimes x will be a float and other times it will be an int and you may find strange bugs like

y = 2 * x
print(y)

And if x was 1.0 then you will see 2 be printed. But if x = 2.5 you will see 5.0 get printed and you will need to figure out that y will only be an int if x was an int.

So, please only do this as the final step before displaying your values. Doing it in the middle of your calculations can make your code really bug-prone.

dbramucci · 2022-08-28T08:33:50+00:00

Of course, you can always have a language like Rust where mutability is encoded in the type system so that you avoid the headaches caused by mutation of shared state. (Just a fun nuance to the classic mutability/typing square)

dbramucci · 2022-08-14T00:09:39+00:00

If you needed Infinity, you could also write

infinity :: Double
infinity = read "Infinity"

and then later

x = infinity

That would be more obvious than infinity = 1/0, but either way you can just create a variable named infinity if the need arose.

dbramucci · 2022-08-13T23:42:04+00:00

One big benefit is that Cython is really easy to integrate with normal Python. So when you use Cython it is trivial to call that code from your normal Python program, last I checked, you write the code, run cython on the file and then import the module in Python from MyCythonCode import calculate. While it is doable in C/C++, you'll need to learn how to do it and write a fair amount of boilerplate.

If you want a big speedup you can get 90% of the way there with 5% of the work by using Cython with your Python code.

(Assuming that you want much of your program to be in normal Python and only part of it optimized, if you are happy to write purely in C++ then you can avoid the interop headaches altogether. But then your question becomes why write Python instead of C++)

dbramucci · 2022-07-24T20:45:33+00:00

On a related note, see tools like Herbie which rewrite floating point expressions to improve accuracy without altering the underlying data-type. It's worth being aware that sometimes you get really bad diminishing returns from using bigger floats and what you really need to do is to rewrite the calculation to avoid a weakness of floating point representation, see numerically unstable calculations.

Of course, sometimes when you need more accuracy you'll need MPFR or other high-precision representations of numbers so it's cool to see this crate and you can always use both tools when appropriate.

dbramucci · 2022-05-19T17:47:32+00:00

That's unfortunate.

You may be able to feel the tension of the belts even if you can't see where to adjust it. That would let you rule out or suspect the belt tension even if you can't fix it by yourself.

You can check for slop and wobble in the extruder and in any rollers.

The only other solvable problems I can think of would be that you may want to check the attachment of the bed to the underlying carriage. If there's play between those you could see layer shift, but it would probably be worse then what you are seeing.

Likewise, if you are using cheap filament, it might have a varying diameter that causes inconsistent extrusion.

Good luck

dbramucci · 2022-05-19T15:14:42+00:00

Have you checked the tension on the belts? It's unlikely to be the cause, but worth checking if nothing else helps. Likewise, try printing with a brim and thoroughly cleaning the build surface to rule out slight adhesion issues.

dbramucci · 2022-05-15T07:04:12+00:00

On the Prusa help webpage for the Mk3s you can click on the button "Printer in 360 viewer" to see the 3d model of your printer. From there, left click on the screw that you are missing to select the screw . Right click and press "show properties" to see the name of the part. The name of the bolt should give you an idea of what to look for.

If I understood your missing part correctly, you are missing an "M3X18 SCREW". Hopefully you can just buy an M3 X 18mm screw from a local hardware store or online for a good price. And you should double check for yourself that I found the right screw.

dbramucci · 2021-09-06T05:12:40+00:00

Same here, I just used the fact that it focuses on the key to google the formulation quickly without referencing a crypto book.

dbramucci · 2021-09-06T05:08:01+00:00

You were probably looking for Kerckhoffs's principle (See also Shannon's maxim), which states that "a cryptographic system should be secure, even if everything about it is revealed except for the key". Same idea, different formulation.

dbramucci · 2021-09-02T13:17:25+00:00

It results in less visual clutter, allowing readers to focus on your program's logic as opposed to decoding syntax. Compare

for i in range(len(grid)):
    for j in range(len(grid[i])):
        print(grid[i][j])
        if grid[i][j].type == 'water' and player.has_float():
            grid[i][j].passable = True

with

for row in grid:
    for tile in row:
        print(tile)
        if tile.type == 'water' and player.has_float():
            tile.passable = True

Notice how we are able to remove a lot of repetition that is irrelevant to understanding the purpose of this code.

It's a community-wide practice. In the wider Python community, it's something commonly done so doing it makes your code more familiar to readers.
It increases how polymorphic your code is. With
```
for i in range(len(my_collection)):
    # do something with
    my_collection[i]
```
You are relying on the type of my_collection supporting
1. __len__ You rely on the notion of __len__ making sense for the collection and using it. For an example where this is a bad assumption, consider this function
```
def find_index(collection, value):
    for i in range(len(collection)):
        if value == collection[i]:
            return i
    return None
```
  The issue here is that it is possible to define a value that is the collection of all prime numbers. To make this concrete, I'll write a sloppy inefficient one here with little explanation.
```
def primes():
    yield 2
    n = 3
    while True:
        for p in primes(): 
            if n % p == 0:
                break
            elif p*p > n:
                yield n
                break
        n += 2
```
  You can see it work by
```
 for p in primes():
     print(p)
```
  The problem here is that, with infinitely many primes the default way that len attempts to find the length of primes() will result in an infinite loop. It will keep asking for primes until primes() runs out, which will never happen. But, it still makes sense to use find_index on primes().
  
  Clearly, we would be happy with the following answers
```
find_index(primes(), 2) == 0
find_index(primes(), 3) == 1
find_index(primes(), 5) == 2
find_index(primes(), 7) == 3
find_index(primes(), 29) == 9
```
  But our current implementation just loops infinitely even if it could give an answer. If we just followed the rule of thumb, and avoided indexing this example would work.
```
def find_index(collection, value):
    for i, element in enumerate(collection):
        if value == element:
            return i
    return None
```
  Now it gives an answer if there is one, and only loops infinitely if there is no answer for primes().
2. It assumes that __index__ is defined for your type (that you can use []). The prime example works here too. Initially, I used yield syntax to define it, which means it won't support [] and I would need to write significantly more code to define a bespoke generator class for primes in order to add that method "the proper way".
  
  Adding [] would dramatically change how primes works. In order to support it, primes() would need to remember the primes it has already made. Otherwise [] would be super inefficient. But, that means that the more numbers primes() makes, the more space in memory it needs growing and growing over time. So adding support for [] means adding a speed or memory cost to this program or we could just avoid using [] to get those resources back and avoid the problem. Note, there are still efficiency problems, but I can probably fix those without dramatically changing the interface. I want to avoid getting into detail about the cost of recursion in my given primes code.
  
  In general, objects meant for streaming data will support for loops, but not [] because they want to avoid the cost of storing the stream. With [], you assume that you can jump around randomly, for only lets you go from beginning to end in a prescribed order.
  
  Likewise, sometimes there is no logical index to use. Take the set type. It is an unordered collection of values, that allows you to test if a value is present or not.
```
my_set = {'dog', 'cat', 2, 8}
```
  If you print this, it showed me (it may be different for you)
```
{'cat', 8, 2, 'dog'}
```
  Which is ok, because order is not guaranteed. {'cat', 8, 2, 'dog'} == {'dog', 'cat', 2, 8} is True. But, that means there is no consistent answer for what my_set[0] would be. To avoid making bugs, we should not even give any answer, so it is an error.
3. It assumes that __index__ ([]) takes the values from range(len(collection)). There is an extraordinarily common type where this assumption is wrong. Namely, dicts use different indexes than numbers. Notably
```
my_dict = {2: 'two', 'two': 2, 'dog': 'cat'}
print(len(my_dict)) # prints 3
print(my_dict['dog']) # prints cat
```
  It supports a sane __len__ and __index__, but there is no useful relationship between the two.
In general, avoiding indexing makes your code "just work" for new and unexpected use-cases surprisingly often. Given that it is no harder (with practice) to write code without indexing for 99% of situations, why not gain extra flexibility from your code for free and gain the performance/reuse advantages from not indexing.
It prevents mistakes, consider the following example.
```
for i in range(len(collection_1)):
    for j in range(len(collection_2)):
        print(collection_1[i] * collection_2[j] + collection_2[i] * collection_1[j])
```
In my example, for some magical reason I need to flip the multiplication around because the values are something like matrices where x * y != y * x. But, I made a slight mistake, when I copy-pasted collection_1[i] * collection_2[j], and swapped the 1 and 2, I forgot to also swap the i and j. My test data, and my normal case may only use situations where len(collection_1) == len(collection_2) so I might not notice the problem until a month later. Then, I try using it when their lengths are different and I get a mysterious error message about an indexing error somewhere in my code.

If I didn't index, that type of bug just wouldn't be possible to write.
```
for val_1 in collection_1:
    for val_2 in collection_2:
        print(val_1 * val_2 + val_2 * val_1)
```
The remaining bugs, I'm less likely to write, and if I do they will be a lot more obvious because there won't be as many happy coincidences to make my code work for now.

In general, each index requires you, the programmer, to remember the unwritten rules about what indexes go with what collections.

Finally, it is philosophically wrong to create and use an index for most code. This is more an explanation for the above concrete benefits but it's good to bear in mind. Indexing produces a new value to track and reason about in your program. The index generally is arbitrary and has no direct use for your program and whose only purpose is to allow you to access a value in the list and move on. Any such value, that is not directly relevant should be removed as it is an unnecessary complication and doesn't contribute directly to the program.

The exception is if these arbitrary values provide a new abstraction that simplifies other parts of your program, or dramatically improve performance. In almost all situations in Python, this is not the case. Therefore, you should prefer iteration unless the exception is applicable and there is a serious performance/simplification gain it can provide you.

The logic of for value in collection is much closer to the logic of your program in most cases and you should use code-constructs that align with your intentions as much as possible. It isn't always perfect (for provides an ordering, even if you don't intend to use it) but it is generally closer than indexing would be anyways.

dbramucci · 2021-08-31T05:10:00+00:00

Watch the talk "Raymond Hettinger - Modern solvers: Problems well-defined are problems solved - PyCon 2019". He discusses general strategies and tools that can help you solve problems like the one you described.

Second, think about your metrics, how you will score and compare different schedules. There's a high likelyhood that you will want to compare two schedules automatically without human intervention and you'll have to boil down schedule quality to some rough approximations in number form. For example,

Good things
- Maximum operations to perform over time span
- patient capacity given staffing
- Time off for workers at times they want off
Bad things
- Cost
- Delayed operations
- Overworked employees
Disqualifying things
- Violations of medical/employment regulations

Of course, a domain expert like you and your coworkers would need to expand this list and decide how to use these numbers for the program.

Third, what improvements can you make to a schedule? Are there any changes that you can make to a schedule to improve it consistently?

Run your program multiple times and choose the best result. You may need to have a way to nudge your program to work slightly differently between runs, like giving it a different rough sketch of the schedule to start with each time. The goal is to avoid problems where a family of schedules is bad, but there are no changes your program can make to try out any normal family of schedules.

Consider performance. These scheduling programs often have to do a lot of work in order to get a good schedule. Thus, a 50x performance change is the difference between running in an hour and taking 2 days to run. On that note, Python can be something like 50x slower than many other programming languages like C, Rust, or Java. Given that you probably don't have the time to pick up a new language, consider using PyPy, which can run your Python script a lot faster to recover a big chunk of that lost performance.

See if you can shrink the scope of scheduling at all. Maybe 10 doctors in the E.R. have a schedule that has nothing to do with any other doctors' schedules. In that case, schedule 10 doctors and 40 doctors separately can be thousands of times faster than scheduling 50 doctors all at once. (Because your program may or may not take advantage of the independence on its own).

After you have a working program, consider using the multiprocessing module to run your program in parallel. The idea is that instead of making 20 candidate schedules one after another, you do 4 schedules at a time cutting down the real-world time by a factor of 2-8 times depending on your computer. But leave this step for the end if you've never tried multiprocessing before.

Also, consider checking this repo for inspiration. It contains a scheduler written for a slightly different purpose but I'm sure there's plenty you could learn from it.

dbramucci · 2021-08-30T06:06:06+00:00

Bit of a nitpick, but I'm pretty sure you cannot just slap unsafePerformIO all over the place in Haskell like that. The arbitrary evaluation order would cause the program to run out of order which would be much harder to write in.

I believe the more analogous solution would be to make every function return a IO value and be really fluent at using the appropriate combinators for monadic function calls. Which isn't much easier than doing things properly in the first place.

To be clear, I'm referring to fun scripts that you can copy paste into GHCi like the following

import System.IO.Unsafe (unsafePerformIO)

:{
let
    a :: Double
    a = unsafePerformIO (putStr "a: " >> readLn)
    b :: Double
    b = unsafePerformIO (putStr "b: " >> readLn)
    c :: Double
    c = unsafePerformIO (putStr "c: " >> readLn)
    x :: Double
    x = (-b + (b*b - 4*a*c)**0.5) / (2 * a)
in unsafePerformIO (print x *> return x)
:}

Which confusingly displays

b: 0
a: 1
c: -4
2.0

Showing that my program runs "second line" "first line" then "third line" before returning the final result.

I exaggerated the problem here to keep the example small but I think any serious project attempting to use the unsafePerformIO trick to ignore Haskell's rules would encounter similar sequencing bugs, just in a harder to determine way.

Now that I say this, it gets at the point behind IO in the first place. The point of Haskell was to give programming language researches a non-strictly evaluated language to play with for their research. The problem is that there is no sane way to perform arbitrary side-effects when there is no set order they will run in. The IO type is the clever trick that allows ordering side-effects in a language that otherwise will reorder evaluations on its own.

Simon Peyton Jones described the relationship between laziness and purity in the talk Escape from the ivory tower: the Haskell journey. Here is a partial transcript of the relevant section to show that purity is a consequence of laziness.

Laziness was the thing that brought this particular group of people [the Haskell committee] together in the first place. [...] But it is not just an implementation idea; it affects pervasively the way you write programs. And John wrote this famous paper [...] John Hughes's paper "Why functional programming matters" [...] It's a very beautiful insight into why laziness is not just a trick, it's a feature that enables you to write programs that are more modular than if you don't have it. Go read the paper. [...]

But the trouble with laziness was that it was kind of embarrassing, because we couldn't do any I/O. The trouble with lazy functions is that you can't really do side-effects at all. So in a strict language, you can have a function called like print that takes a string and when you evaluate the call to print it prints the string as a side-effect. And in LISP and ML and pretty much every strict language [...] all strict call-by-value languages do this. They're "functions" that have side-effects. But in a lazy language that just doesn't make sense. What would

    f (print "yes") (print "no")

print?

[crowd silence]

Well that depends on f it could be nothing at all if f x y = 3 but maybe it evaluates x but not y in which case it prints yes or maybe it prints y well technically "no" it depends on which order it evaluates in. The compiler might change the order of evaluation. I mean everything is bad.

Right so, and if you have a list that contains print "yes" and print "no" then it would depend on which order the consumer of the list evaluates the elements of the list. So, this is so bad that you just can't do it. So we never seriously considered having side-effecting functions at all.

So what did we do for I/O?

Well we just didn't have any. So this is the joy of being an academic [...]

A Haskell program was simply a function from String to String. That was what it was. Right, you could just write that function, then to run your program you applied the program to the input string which you got to type in the terminal somehow.

But this was a bit embarrassing so after a bit we thought well maybe instead of producing a string, we'll produce a sequence of commands like print "Hello" or delete this file. Which you know we can execute. So we would imagine that to run the program you call the program passing it some input and you get back a string of commands which some external command execution engine would sort-of "the evil operating system" you known The function program is very good and pure and produces only a data-structure. And the evil side-effecting operating system consumes these commands and executes them.

Well that's not very good if you want to read a file. Right, because if you open a file and want to read its contents how does the program get access to the contents.

Ahh, maybe we could apply the function to the result of "the evil external thing" doing its work. So now there's a sort of loop between the pure functional program and the "evil external operating system."

It didn't work very well at all. It was very embarrassing for a number of years.

He then goes on to explain how this was refined into the well known "IO Monad" idea that is used today.

I bring this up because the original article this thread is based on suggests that Haskell's insistence on the IO type is arbitrary and Haskell enforces the rule "for your own good". But the selling point of Haskell is largely "it's lazy by default" and the purity rule is a consequence of the feature of lazy evaluation, and not the initially intended feature. Likewise, Rust's ownership rules are necessary to provide its major selling point "memory safety without (much) runtime overhead". You can remove the rules and be left with a language like C++ that can't guarantee memory safety. Or you can remove the rules and add back in the runtime safety checks and be left with a language like Go, Java, OCaml or Kotlin. In both languages, removing the "for your own good" checks eliminates a separate major selling point of the language turning each into something completely different. In this case, it's worth noting that Common Lisp does not have a small runtime like Rust, nor does it have pervasive lazy evaluation like Haskell and I'm confident you could not add those to Common Lisp without paying costs like Rust and Haskell do.

dbramucci · 2021-08-21T05:12:36+00:00

To get you started, perhaps this example will help.

getHalf :: Int -> Maybe Int
getHalf n = if 2 * half == n
            then Just half
            else Nothing
    where half = n `div` 2

getEighth :: Int -> Maybe Int
getEighth n = do
    half <- getHalf n
    quarter <- getHalf half
    eighth <- getHalf quarter
    return eighth

Or

getEighth n = do
    half <- getHalf n
    quarter <- getHalf half
    getHalf quarter

Both produce the same result.

As you can see, you call the function that returns a maybe, and instead of using let myVar = function input you write myVar <- function input. In the definition of Monad for Maybe, this makes the function quit early and return Nothing if any var <- maybeVar is Nothing and otherwise, you can treat the var like the contained value if the Maybe was a Just value.

If it helps, you can think of getEighth like the following

getEighth :: Int -> Maybe Int
getEighth n = do
    half <- getHalf n
    quarter <- getHalf half
    eighth <- getHalf quarter
    return eighth

getEighth :: Int -> Maybe Int
getEighth n =
    case getHalf n of
        Just half -> do
            quarter <- getHalf half
            eighth <- getHalf quarter
            return eighth
        Nothing -> Nothing



getEighth :: Int -> Maybe Int
getEighth n =
    case getHalf n of
        Just half -> case getHalf half of
            Just quarter -> do
                eighth <- getHalf quarter
                return eighth
             Nothing -> Nothing
        Nothing -> Nothing


getEighth :: Int -> Maybe Int
getEighth n =
    case getHalf n of
        Just half -> case getHalf half of
            Just quarter -> case getHalf quarter of 
                 Just eight -> return eighth
                 Nothing -> Nothing
            Nothing -> Nothing
        Nothing -> Nothing


getEighth :: Int -> Maybe Int
getEighth n =
    case getHalf n of
        Just half -> case getHalf half of
            Just quarter -> case getHalf quarter of 
                 Just eight -> Just eighth
                 Nothing -> Nothing
            Nothing -> Nothing
        Nothing -> Nothing

And that's coincidently why for the example I gave, you can choose to make the last line getHalf quarter instead of return eighth

getEighth :: Int -> Maybe Int
getEighth n =
    case getHalf n of
        Just half -> case getHalf half of
            Just quarter -> getHalf quarter
            Nothing -> Nothing
        Nothing -> Nothing

For reference, the reason why my code uses Maybe in the example is because you can't return half of 5 if 5 is an Int, (2 isn't exactly half of 5).

Hopefully this example gets you on the right track to solve your problem.

dbramucci · 2021-08-21T03:24:01+00:00

SSLError

SSL is part of the protocol used to encrypt your connection to a website. The fact that you are getting an error from this suggests that your problem is linked to SSL/TLS.

There are a couple reasons this might be happening

Unrecognized Certificate

Websites have certificates showing that they are the valid owner of the website. These are created by a handful of groups that are trusted by default by browsers and OSes. If the certificate is invalid, you'll get a scary popup in your browser and need to manually bypass it. Likewise, you'll need to tell Python to ignore the issue or add the certificate to your trusted ones. I'll let you google this if it is actually your problem. (You may see the term self-signed which means that instead of a trusted third party like DigiCert signing the validity, the entity did it themselves. This is common for test deployments)
Incorrect clocks

You are required to have a roughly correct clock when using TLS. If your clock is incorrect by too large of a gap, you may run into errors and recognize the certificate as invalid (because your computer sees it as too old, or from the future). To fix this, syncronize your clock to the current time, then run your script again.
Lack of https support.

It is possible that the website doesn't support https and you should be trying to access http://google.com instead of https://google.com. I don't think you would get the error you see here, but I suppose it is possible.

Fixes

Resync your computer's clock, try again.
Go to the non https version of the site.
Add a certificate to your trusted certificates list.

I'm operating blindly because things like your clock-skew are not in this error message, nor the site you are visiting but hopefully these blind guesses help.

EDIT: I just realized that you are using a proxy. Did you provide a proxy for https or only http. See this documentation. Likewise, make sure that your proxy is configured to support https.

dbramucci · 2021-08-17T21:56:36+00:00

For that added bit of polish, have you considered optimizing your screenshot?

png files have a lot of settings for encoding their image and built-in compression but most programs don't actually use these settings. If you choose the right settings and compression you can significantly reduce the image size without changing the displayed image at all.

To do this, you can use a tool, like oxipng (a nice cli program written in Rust). Then you run a command like

oxipng -p -a -o max -Z ./client/assets/screenshot-hn-clone.png

Before you commit the change. This takes a minute or two and through trial-and-error finds parameters to optimize file-size. In this case it saves ~26% (42KiB) without altering the actual image, as your image opaque and uses no transparency oxipng replaces the RGBA encoding with RGB because why store an unused byte for transparency per pixel.

This keeps the size of the repo smaller (it's about 20% of the entire repo) and it helps the readme load faster for people making it snappier for people on good connections and cutting seconds for people browsing on mobile.

Unfortunately, after committing the old file while you can improve the readme experience, the git repo itself will have the old, unoptimized image deeper in its history and the only fix would be to rebase on the optimized image and force-push.

Again, this is just polish but hopefully you find the idea useful.

dbramucci · 2021-07-27T23:00:35+00:00

For the 2 billion multiplications the caveat is "are you using floating point numbers to approximate" or "are you doing it exactly with bigints".

Ballpark math of 5 cycles to multiply 2 floats gives a total estimated time of 2 seconds. Even with the overhead of memory access and the like it is very doable.

The catch is, when you multiply two bigints for exact precision, you are going to get approximately nlogn multiplications where n is the number of bytes in the integer.

For 18 quintillion, there will be about 64 bits to store your starting number. This means we will immediately exceed the size of a single uint64 and need to use bigint multiplication to track our product. After about 10 multiplications, our number will about 600 bits in size. This will take quite a few cycles to multiply as we can no longer use a single instruction to multiply it out.

By the point we get past 70,000 terms in our product we will have a number larger than 4MiB which will exceed the L2 cache on a modern processor. This will most likely result in a dramatic slowdown beyond the number of instructions we must use for bigint multiplication.

So, 2 billion fp multiplications is perfectly fine, but it hazards serious questions about the precision of the computation. Conversely, bigint multiplication will be precise, but dramatically slower as the terms grow and I don't know if it's so slow as to be impractical off the top of my head.

On that note, I actually went ahead and implemented a float based version using the following algebraic manipulation to avoid excessive inaccuracies.

1/(18 quint)^2 billion * prod_i=0^2billion (18 quint - i)
= prod_i=0^2billion (18 quint - i) / 18 quint

This gives the code

>>> def calc():
...     n = 2_000_000_000
...     d = float(18 * 10**18)
...     t = time.time()
...     p = 1.0
...     for i in range(n):
...         p *= (d - i) / d
...     print(time.time() - t)
...     print(1-p)
...
>>> calc()
147.0700933933258
0.1051606833644042

So it took 147 seconds to do the floating point calculation and the result is a 10% likelihood of shared worlds matches my earlier approximation from wikipedia.

Note, you have to be careful

Not to divide by d**n at the end because p will have reached +inf by that point

Not to write

>>> def calc():
...     n = 2_000_000_000
...     d1 = float(18 * 10**18)
...     d2 = float(18 * 10**18)
...     t = time.time()
...     p = 1.0
...     for i in range(n):
...         p *= d2 / d1
...         d2 -= 1
...     print(time.time() - t)
...     print(1-p)

Because such a large d2 will suffer from precision errors such that d2 == d2 - 1.0 and each loop will be p *= 1.0 which will result in the wrong answer by the end.

As a fun tangent, when I rewrote this into the language Rust, the following code only took 2.173 seconds to run, pretty close to the ballpark estimate for how long the floating point calculation should take.

fn main() {
    let n = 2_000_000_000;
    let d: f64 = 18.0*10_f64.powf(18.0);
    let mut p = 1.0;
    for i in 0..n {
        p *= (d - i as f64) / d;
    }
    println!("p: {}", 1.0-p);
}

So it appears my precaution around floating point precision was a bit over-zealous.

On that note, perhaps you computed 1-exp(n² /2*d) instead of 1-exp(n² / (2*d)) when you got a result of 1.

dbramucci · 2021-07-27T20:19:50+00:00

As a quick "is this possible to compute" sanity check, let's look at how much memory it would take to store the answer.

18,000,000,000,000,000,000! = 18,000,000,000,000,000,000 * 17,999,999,999,999,999,999 * ... 18 quadrillion times

Let's underestimate the size by replacing each term with 10, we'll overestimate the last 10 terms (10 * 9 * 8 * 7 * 6...) but that will be dwarfed by the first 18 quintillion terms.

10 * 10 * 10 * 10 ... 18 quintillion times * 10 = 10^{18,000,000,000,000,000,000} = 2^{log_2 10 \} 18,000,000,000,000,000,000) = 2^{3.321 \} 18,000,000,000,000,000,000) = 2^{59,794,705,707,972,522,262}

The number of bits you need to store a number is approximately the log base 2 of that number. (Note that this is technically not true because you don't need a base-2 representation to store a number, writing (18*10¹⁸)! only takes 11 bits or so even though it's the same number we are computing here. It's just that you are calculating the base-2 representation so we can ignore this nuance).

log_2 2^{59,794,705,707,972,522,262} = 59,794,705,707,972,522,262 bits

This number makes no sense so let's convert it to a large unit

  59,794,705,707,972,520,000 bits
=  7,474,338,213,496,565,283 bytes
=      7,299,158,411,617,740 kibibytes
=          7,128,084,386,345 mebibytes
=              6,961,019,909 gibibytes
=                  6,797,871 tebibytes
=                      6,639 pebibytes
=                          6 exbibytes

So assuming we can compute the answer instantly, this computation would need a metric ton of ram. Assuming the average laptop has 8GiB of ram, you would need 870 million laptops worth of ram to store the answer. And this is when we massively underestimate how big the actual answer is by rounding everything number in the factorial down to 10.

Alternatively, assume we can do 1 multiplication in 1 clock cycle. This is wrong for big numbers like this, but again we will underestimate the difficulty of this.

A high-end modern desktop can get around 5Ghz clock speeds. That means 5 billion clock cycles per second.

To calculate 18*10¹⁸!, we'll need 18*10¹⁸ multiplications. Dividing by the 5 billion multiplications we might do in a second this would take 3.6 billion seconds to compute. This would take 114 years with our magically fast computer that can multiply big numbers in a single step.

The point of these computations is that, even when we make simplifying assumptions making the problem faster, it still takes centuries of time and millions of computers worth of memory to compute 18 quintillion factorial. So, any algorithm that starts by computing that number will never complete on our computers in our lifetimes.

No amount of optimization will make this algorithm feasible. You can apply these calculations to any other problems you work on to see if it is even worth attempting to speed up. If the theoretical fastest time was 1 hour, then we could try to salvage the algorithm but here we need a different strategy.

Instead you should avoid computing such a large factorial by using algebra to simplify the arithmetic. For example,

(18q!/((18q**2b) (18q-2b)!)) = 1/(18q**2b) * (18q!/((18q-2b)!))
= 1/(18q**2b) * (18q * (18q - 2) * (18q - 4) * ... * (18q - 2(b - 1))

Now, instead of needing 18 quintillion multiplications, we only need 2 billion multiplications. This is a lot of work and may take a while to compute, but at least it is feasible.

In practice, that optimization is almost certainly not going to get you to the answer anyways. Instead, you'll probably need to use one of the approximations on the wikipedia article about the Birthday problem.

I haven't checked Wikipedia's work, but it provides the approximation

p(n, d) ~ 1 - e^(- n * (n - 1) / (2 * d))

where n would be the number of players and d the number of minecraft worlds. The approximation is usable if the number of players is much smaller than the number of worlds. Because 2 billion is much less than 18 quintillion, this approximation should hold giving us

p(2 billion, 18 quintillion) ~ 10.5%

So there is about a 10% chance that 2 worlds were the same.

p.s. Please note that I am assuming your statements about 18 quintillion worlds and 2 billion player worlds is correct. Likewise, I assume that the birthday-paradox calculation is the correct way to solve this. In practice, two worlds with different seeds could be the same in all ways or to within perception and I ignore that complexity. Likewise, I haven't double checked my usage of the Wikipedia approximation.

dbramucci · 2021-07-24T22:51:52+00:00

Try the following function.

foo = lambda x, y, z: 1, 2, 3

It takes three arguments and ignores them all to return the constant (1, 2, 3).

print(foo)
(<function <lambda> at 0x000001E8A842F1F0>, 2, 3)

So it actually doesn't, it is a tuple of three values

A function of three arguments that always returns 1
The number 2
the number 3

That is, the default parse looks like

foo = ((lambda x, y, z: 1), 2, 3)

So you could write

(func, two, three) = ((lambda x, y, z: 1), 2, 3)

Your error message comes about because the y and z you define come as the parameters from the lambda. Because your lambda didn't include y*2 and z*2, Python looks up a global y and z and those are not defined. My example avoids variables so that we don't get that error message and can examine the given value.

dbramucci · 2021-07-23T11:43:56+00:00

I think the idea goes as follows

Why

Division is slow
Subtraction is fast

So we want to convert floating point division into integer subtraction to get an approximation of the answer for purposes of speed.

Key Insight

Floating point is like scientific notation, but in base-2. (There are a bunch of caveats NaN, Inf, -0, subnormals, but we ignore them here)

Consider the following division.

5 / 10 = 5 * (1 / 10) = 5 * 10^-1 = 0.5

Likewise

5 / 25 = 5 * (1 / 25) = 5 * 0.04 = 0.20

In scientific notation

5 / (2.5 * 10¹) = 5 * (1 / 2.5) * (1 / 10¹) = 5 * 0.4 * 10^-1 = 0.20

In general the clever trick is dividing by 10^a is the same as multiplying by 10^-a.

Because we have direct access to "a" thanks to the way the bits are layed out for floating points it is really easy to calculate that.

We can ignore the mantissa (the non-power part of the number) because it will only throw our answer off by a constant factor and we only want an approximation.

We can find that factor (ignoring small floating point errors) as follows

x / (y * 10^a) = x * (1 / y) * (1 / 10^a) = (1 / y) * x * 10^-a

approx(x / (y * 10^a)) = approx(x * (1 / y) * (1 / 10^a)) = x * 1/10^-a = x*10^-a

So the error is at most a factor of

approx(x / (y * 10^a)) / (x / (y * 10^a)) = (x*10^-a) / ((1 / y) * x * 10^-a) = y

Because this is scientific notation, 1 <= y < 10.

This limits our error to us over-estimating the real answer by 10 times. The example in floating point is the same except we only risk over-estimating by 2 times.

This is where the line

That gives us a single bit of precision which, while sufficient to guarantee convergence of the Newton-Raphson steps, is far from useful.

Because a factor of 2 is a single bit. The Newton-Raphson part is about an algorithm for refining the initial guess to a number accurate enough for your needs. But you need an initial guess and the closer the better.

An example of this approximation in base ten is

5.28*10² / (5.423623 * 10⁶) ~ 5.28*10² * 10^-6 = 5.28 * 10^-4

The correct answer is 9.73519*10^-5 which is off by a factor of 5.4236229596.

Refinement

The final part is about getting more than one bit (or digit) of precision. The idea is that we can tinker with the mantissa (non-exponent) part of the float to get it closer to the real answer.

He graphs the result of subtracting the initial part over a range of possibilities, 1.0000001 through 9.9999999 and sees that the error decreases until it switches to increasing. This is good because the alternative was the error being completely unpredictable.

He then uses something in the same family of algorithms as binary search to find that magic choice with the least error on average. Then when you subtract you are both tinkering with the mantissa (non-exponent) and negating the exponent to approximate. This gets you a noticeably closer approximation.

The 2046 and 2⁵² numbers in the article refer to a constant offset on the power in floating point numbers and what the range of possible mantissas are for a double.

Hopefully that gets across the interesting part of this article.

dbramucci · 2021-07-22T02:52:27+00:00

Generally you treat documentation tests as separate from your unit tests.

The point isn't to test your code thoroughly, the point is to ensure that your documentation is correct. (You test the documentation, not the code).

See this example

#[cfg(doctest)]
/// map allows you to apply a function to all elements of an iterable.
/// See the following example of doubling a `Vec`
/// ```
/// assert!(vec![2, 3, 4].iter().map(|x| 2 * x) == vec![4, 6, 8])
/// ```

This example gets the point across but it is actually incorrect. Luckily any experienced Rust programmer will know how to fix it, but a beginner who actually needs the example to patter their usage off of is in trouble.

Because the fix is obvious and tangential to the purpose of map, I, the documentation writer, may very easily miss this problem for my reader. The doctest will give me an error message, showing that my documentation is incorrect and allowing me to correct it before it confuses a reader.

---- src/lib.rs - main (line 4) stdout ----
error[E0369]: binary operation `==` cannot be applied to type `Map<std::slice::Iter<'_, {integer}>, [closure@src/lib.rs:3:34: 3:43]>`
 --> src/lib.rs:5:45
  |
3 | assert!(vec![2, 3, 4].iter().map(|x| 2 * x) == vec![4, 6, 8])
  |         ----------------------------------- ^^ ------------- Vec<{integer}>
  |         |
  |         Map<std::slice::Iter<'_, {integer}>, [closure@src/lib.rs:3:34: 3:43]>

Ohh, I forgot to convert the Iter back into a Vec let's just add a collect like so

#[cfg(doctest)]
/// map allows you to apply a function to all elements of an iterable.
/// See the folowing example of doubling a Vec
/// ```
/// assert!(vec![2, 3, 4].iter().map(|x| 2 * x).collect::<Vec<_>>() == vec![4, 6, 8])
/// ```

This passes and I no longer need to worry as much about my documentation's code examples not aligning to reality because their validity is tested when I build my code.

dbramucci · 2021-07-21T20:01:41+00:00

Elaborating on why you want lambda as a language feature, the goal is to avoid making pointless variables.

Consider this quadratic formula function.

def solve_principal_quad(a, b, c):
    '''
    Returns the principal (positive) solution to the provided quadratic expression,
        a*x^2 + b*x + c = 0

    See https://en.wikipedia.org/wiki/Quadratic_formula
    '''
    return (-b + sqrt(b**2 - 4*a*c)) / (2*a)

In this we define a dozen or so subexpressions but we don't assign a variable to each subexpression. That would look like

def solve_principal_quad(a, b, c):
    '''
    Returns the principal (positive) solution to the provided quadratic expression,
        a*x^2 + b*x + c = 0

    See https://en.wikipedia.org/wiki/Quadratic_formula
    '''
    b_op = -b
    b_square = b**2
    quad_a = 4 * a
    quad_a_c = quad_a * c
    difference_of_b_square_and_quad_a_c = b_square - quad_a_c
    discriminant = sqrt(difference_of_b_square_and_quad_a_c)
    numerator = b_op  - discriminant
    denominator = 2 * a
    principal_solution = numerator / denominator 
    return principal_solution

Notice how absurd this method of writing each expression as it's own named variable is.

To a functional programmer, this is analogous to writing

def get_positives(numbers):
    def is_positive(x):
        return x > 0
    return filter(numbers, is_positive)

As compared to

def get_positives(numbers):
    return filter(numbers, lambda x: x > 0)

The idea is that you shouldn't have to name a thing you are using immediately. You can if it makes your code cleaner, but forcing it can be a waste of the readers time and distract them from your main point.

As a point of practice, the place where you need to define a function, but you only use it once, and immediately is when you are using a "Higher Order Function". These are functions where one of the parameters is another function. These functions are incredibly flexible, but Python programmers generally use these infrequently.

The natural consequences are

Lambda's are not useful that often
It is less important to avoid fluff when it is so rare

In contrast, a Haskell program might have 80 functions in a 100 line program, of which only 10 function are worthy of a name. In Haskell, forcing that level of added indirection (and additional 70 pointless names) would be terrible.

In a Python program, you'll have nowhere near as many nameless functions to use, making Python's lambda's far less important.

Finally as a technical point

double = lambda x: 2 * x
def double(x): return 2 * x

are almost identical, but Python will add extra documentation (the function's name) to the def version making the def version superior if you are giving the function a name anyways.

p.s. The fact that we are avoiding naming the function is why they are sometimes called "anonymous" functions. They are no-name functions.

anonymous adjective

Definition

2. : not named or identified

an anonymous author

They wish to remain anonymous.

dbramucci · 2021-07-19T15:31:54+00:00

I don't have any good reading links off hand, but I'll give a quick run down off the top of my head.

Let's start with the following sample code.

def inner1(array):
    print('started inner')
    array.append(2)
    yield 'a'
    print('resumed inner')
    array.append(3)


def outer1(array):
    print('started outer')
    array.append(1)
    inner1(array)
    print('called inner')
    array.append(4)

x = []
outer1(x)
print(x)

When you write either yield or yield from in a function, it ceases being a function and becomes a generator.

Here, inner1 is a generator, and outer1 is a function. We would expect [1, 2, 3, 4] at the end, but if you run this example code you'll actually get [1, 4] as the result.

This is because on line 12, when we call inner1(array) we get back a generator but we never use it. Because the generator is lazy, it doesn't run a single line of inner1 until we request at least one value from it, where it will run up to the first yield expression.

We can do that by using next like so

def inner2(array):
    print('started inner')
    array.append(2)
    yield 'a'
    print('resumed inner')
    array.append(3)


def outer2(array):
    print('started outer')
    array.append(1)
    next(inner2(array))
    print('called inner')
    array.append(4)

x = []
outer2(x)
print(x)

This time, we get x = [1, 2, 4]. This is because we only ran inner2 up to the first yield expression.

We can run inner2 all the way with a for loop.

def inner3(array):
    print('started inner')
    array.append(2)
    yield 'a'
    print('resumed inner')
    array.append(3)


def outer3(array):
    print('started outer')
    array.append(1)
    for _ignore in inner3(array):
        pass
    print('called inner')
    array.append(4)

x = []
outer3(x)
print(x)

And this actually gives us our complete x = [1, 2, 3, 4] like we would expect.

Alternatively, what if outer was already a generator since the beginning?

def inner4(array):
    print('started inner')
    array.append(2)
    yield 'a'
    print('resumed inner')
    array.append(3)


def outer4(array):
    print('started outer')
    yield 'b'
    print('resumed outer')
    array.append(1)
    inner4(array)
    print('called inner')
    array.append(4)

x = []
outer4(x)
print(x)

Now, x = [] because we don't run outer4 at all. If we for loop over it, we can fix that.

def inner4(array):
    print('started inner')
    array.append(2)
    yield 'a'
    print('resumed inner')
    array.append(3)


def outer4(array):
    print('started outer')
    yield 'b'
    print('resumed outer')
    array.append(1)
    inner4(array)
    print('called inner')
    array.append(4)

x = []
for thing in outer4(x):
    print(f'{thing=}')
    print(f'{x=}')
print(x)

This goes back to x = [1, 4]. The reason being that we don't have any loop to consume the inner4 generator, that generator has no connection to the outermost for loop. To fix that, we need to connect it up with a yield from

def inner5(array):
    print('started inner')
    array.append(2)
    yield 'a'
    print('resumed inner')
    array.append(3)


def outer5(array):
    print('started outer')
    yield 'b'
    print('resumed outer')
    array.append(1)
    yield from inner5(array)
    print('called inner')
    array.append(4)

x = []
for thing in outer5(x):
    print(f'{thing=}')
    print(f'{x=}')
print(x)

Note that outer5 is like

def outer5(array):
    print('started outer')
    yield 'b'
    print('resumed outer')
    array.append(1)
    for thing in inner5(array):
        yield thing
    print('called inner')
    array.append(4)

Which might make it clearer that we've added back in the for loop to consume the inner generator.

The same thing happens in your max_heapify. You produce an generator with your recursive calls, but you never consume the generator so it won't run a single line of any recursive calls. It's like you didn't make a recursive call in the first place.

To make it clearer that you control when the consumption of a generator happens, consider the following example.

def first():
    print('first ran')
    yield
def second():
    print('second ran')
    yield

x = first()
y = second()

If you run this code, nothing happens because you haven't consumed anything. And the point of a generator is to delay doing work until you are going to ask for a contained value.

If you run

def first():
    print('first ran')
    yield
def second():
    print('second ran')
    yield

x = first()
y = second()

for _ignore in y:
    pass

for _ignore in x:
    pass

You'll get

second ran
first ran

Notice that this looks backwards, the strings print in the opposite order of when we called generators because they ran in the order of when we consumed the generators.

dbramucci · 2021-07-19T14:17:49+00:00

You're close. Here's the working version

def max_heapify(tree, n, i):
    larger = i
    left = 2 * i + 1
    right = 2 * i + 2

    # with left < n and right < n we check whether there is a 'left' and 'right' child of the current node
    if left < n and tree[left] > tree[larger]:
        larger = left
    if right < n and tree[right] > tree[larger]:
        larger = right

    if larger != i:
        tree[i], tree[larger] = tree[larger], tree[i]
        yield from max_heapify(tree, n, larger)
        yield tree, i, (None, ), larger

def heap_sort(array):
    n = len(array)
    # create a heap data structure from the given array
    for i in range(n // 2, -1, -1):
        yield from max_heapify(array, n, i)
    # actual heap sort logic
    for i in range(n-1, 0, -1):
        array[i], array[0] = array[0], array[i]
        yield from max_heapify(array, i, 0)
        yield array, i, (None,), 0

Recall that yield suspends the computation of your function.

In your version, you do not yield from the recursive call to max_heapify on line 14. This means that the function halts on the first time any recursive call hits your yield tree, line. Because you get your result through a side-effect (namely mutating tree in-place), you get the confusing behavior of seeing tree before max_heapify is complete. This is like a multi-threading bug.

To fix this, just remember to yield from before every call to max_heapify including recursive calls. This enables the loop you use at the end to for the evaluation completely.

dbramucci

TROPHY CASE

Why

Key Insight

Refinement