Julia sets: True shape and escape time

  • 133 Replies
  • 4319 Views

0 Members and 1 Guest are viewing this topic.

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #120 on: January 24, 2020, 09:55:13 AM »
While I was rearranging my IA expressions, I came across this article. The authors rearrange their expressions in an automated manner with respect to rounding error and evaluation costs. So not directly applicable to my goal of reducing overestimation but quite an interesting optimization problem there. Maybe I can use some of their technqiues.

Code: [Select]
SOAP: Structural Optimization of Arithmetic Expressions for High-Level Synthesis
Xitong Gao, Samuel Bayliss, George A. Constantinides

http://cas.ee.ic.ac.uk/people/gac1/pubs/XitongFPT13.pdf

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #121 on: January 26, 2020, 07:01:44 PM »
"A long overlooked optimization opportunity"

A)
While I was answering pauldelbrot's last question I remembered that IA libraries often use sign checks to speed up interval multiplication. E.g. if for [a..b]*[c..d], both intervals lie on the positive axis, the resulting interval will be [a*c..b*d], so one can save 2 multiplications.

Implementing this gave a ~10% speed increase.

B)
Rearranging IA expressions I realized something I have embarrassingly overlooked completely for the last 6 months.
(Reminder: My TSA routine uses the complete memory for cell color storage, 2 bits/cell, so bounding boxes are repeatedly computed as needed. Storing a precomputed bbx would need at least 10 bytes, thereby reducing the maximally reachable refinement level by six, so that was not an option)

For e.g. z^3+A*z+c the expanded component function for the real part is: d + f*x - g*y + x^3 - 3*x*y^2
(A=f+g*i, c=d+e*i)

My routine goes over the image with an outer y-screen coordinate loop and an inner x one.

But (and that's the overlooked part) there are subexpressions which are independent of x, so need not be computed at every bbx call, e.g. d-g*y or y^2. So I could move this computation out of the x loop.

And one could even move those and expressions like f*x+x^3 out of any loop and precompute those once before the actual image calculation starts (as obviously a pixel always represents the same complex interval).

The only expressions to be evaluated at a bbx call are then expressions like x*y^2 or the sum of f*x-g*y which are dependent on both x and y at call time.

That way quite a lot of flops could be avoided, making those expressions having linear (in screen width) call number complexity as opposed to at least quadratic as before.

Using those precomputed values (called "helper") led to a further decrease in running time for a cubic polynomial of 35% for f107. Reduction in time increases relatively when using higher precision datatypes (f161), but are negligible for C++ double (like with almost any attempt optimizing for that datatype).

The image below shows z^3+A*z+c with
c=(2097152+26214400*i) * 2^-25
A=(8388608-3355444*i) * 2^-25
depicting the newly identified interior cells with a rearranged IA formula (blue), exterior (red), and, although not unexpected, not only increasing the black area but also finding new regions (orange arrow, zoomed).

However there's a coding price to pay: I have to manually code the bbx routine and the precomputation parts as opposed to the easy-to-use automatic output of my IA parser I used before. But for about 50% speed up in total, that's worth the once-to-do heavy work.

Now on to the other formulas and especially the 3D ones (new release of juliatsa-core to come when I'm finished).

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #122 on: January 28, 2020, 04:22:06 PM »
"One level earlier"

The first example, where rearranging the IA expression results in detection of a cycle one refinement level earlier than with the classic expanded form. Currently I use the heuristic rule: Removing complete factors from an expression and moving them outside brackets (subdistributivity law).

alpha increases from 0.509 to 0.645 at refinement level 10.

z^5+A*z+c
( A := f+g*i, c = d+e*i )
c=(-14680064-2097152*i) * 2^-25
A=(-34209751+0*i) * 2^-25

E.g. the real part of the bbx was rearranged as:
Code: [Select]
expanded      d - g*y +   f*x + 5*x*y^4 +   x^5 - 2*5*x^3*y^2
rearranged    d - g*y +   x*(f+5*y^4) +     x^3*(x^2-2*5*y^2)


Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #123 on: January 29, 2020, 09:52:40 AM »
"The wrapping effect"

A nice example of two features of IA.

The left image shows the clearest overestimation I've encountered so far: ~1600 pixels in the bbx - and all (details in post #116) exert influence on the cell A whose bbx is depicted. Chances are high, this cell A will never be judged as definite interior or exterior, as one "disturber pixel" in the bbx suffices to impede that.

It shows the wrapping effect of IA calculations mentioned in the literature: Starting with a square parallel to the coordinate axis, one often ends with an angled square of true iterates in the bbx, so even if the bbx were tight here, it would still incorporate a number of unnecessary pixels.

The right is the contrary - almost optimal, only 4 pixels in the bbx, and all contain an iterate.

Used here is the expanded form of z^6+A*z+c with an alpha of 0.467. The best rearranged version has a 0.5 (so only 10% more coverage). I hope I can find something better.

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #124 on: February 02, 2020, 10:52:39 AM »
"Optimization pipeline: The 0-step-approach"
(just text, result: non-improvement)

The squares the complex plane is tiled into for the TSA are closed, so adjacent squares share an edge or at least a corner. If a gray cell lies next to (also diagonally) a white square, at least the shared corner is in the exterior, so the gray cell definitely has a path to white and can be set to "potentially-white" before ever calculating a bbx. And since the TSA goes over (part of) the image repeatedly, setting a cell at the earliest stage will prevent it from being checked in as many subsequent loops as possible.

I chose 4 sets for a test: the basilica as one of the easiest sets to compute, a parabolic one (slow dynamics), the rabbit as an intermediate and a disconnected one to represent the range of possible sets.

I computed all sets to level 15 in the standard way. The refinement step to level 16 was carried out once with ("preset") and once without ("plain") the 0-step-approach. Total computation time was recorded, as was the total number of bbx's calculated.

The idea of presetting sounded to me like a possible improvement, but turned out to be - none, time-wise even the contrary:

Code: [Select]
                         time to compute (sec)                 bbx's computed
                         plain preset rel%             plain            preset      rel%
basilica c=-1+0*i         31     39    +25%         86 773 694        84 668 326    -2,4%
rabbit c=-0.25+0.75*i     47     55    +14%        303 547 648       296 845 861    -2,2%
parabolic c=-0.75+0*i    237    255    +7%       3 080 708 208     3 078 280 433    -0,1%
exterior c=0.25+0.75*i    29     37    +27%         68 340 679        66 539 010    -2,6%

rel% := (preset-plain)/plain in percent

All sets took considerably longer to finish, even though the number of bbx's computation went down as expected.

I figure since the presetting is proportional to the circumference of the set, the number of gray cells however to the area, in higher refinement levels the relative abundance of presettable squares decreases. And additionally there is the cost of going once over the image - even if no bbx is computed, there's still a lot of memory transfer going on, and since I pass over the image row by row, there will be continuous cache misses.

Nevertheless, I find this presetting conceptually quite appealing: Point-sampling computes the whole orbit ("all-step"), the TSA by Figueiredo computes a bbx ("1-step", and backpropagation) and now there's a "0-step" to judge whole squares in the complex plane (however not to a definite color though). The logical next step would then be a quantum computer that knows the answer before I decided what to compute :)

I may still follow this idea and interweave the presetting with loading the raw data, as I have every row in memory there at one point anyways, so no additional memory transfer cost. Maybe that would lead to a (slim) speed improvement.

(Note: I use the terms square/cell/tile interchangeably, sometimes also pixel, as well as the terms raw data/image/2-square).

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #125 on: February 03, 2020, 09:20:17 AM »
"Convergence pipeline: Subdivision - Step 1"

Virtually gridding a square whose bbx is to be computed, usually leads to a smaller overall rectangle, encompassing the iterates ("subdivision", as it is called in the literature) and hopefully to faster convergence.

However gridding during the main loop is not applicable due to very high time costs - gridding a cell into 8x8 means about 64 times the number of bbx to compute.

So I decided to precompute a "shrinkage value" per pixel. It is a small integer the bbx screenrect can be shrunk symmetrically towards the center while still encompassing all iterates. This value will be computed at the desired refinement level once before entering the main loop and stored. Within the main loop, it can be reused every time a specific bbx screenrect is computed to tighten it. I chose to only store an 8-bit integer for a gray pixel to keep memory usage to a minimum, therefore the symmetrical part.

The first image shows for a quartic polynomial 3 bbx screenrects with the shrinkage value (smallest length of a blue line) of such a 8x8 virtual grid of the pixel. For the rightmost example that meant a reduction of the number of bbx pixels from 49 to a mere 9. I hope the reduction will even be more substantial in higher order polynomials.

Initially I intended to compute the shrinkage value at the target refinement level for as many gray pixels until the assigned memory was exhausted. But visualizing it at two test levels surprisingly showed, that the regions of positive shrinkage values are quite constant (2nd image): gray cells have a shrinkage of 0, green ones allow a value of 1 and red ones (upper right) of two.

Next step: This phenomenon can be used as a filter: identifying those regions in a low refinement level and using that information to guide computation in the target level to only a subset of not-yet colored complex plane squares.


Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #126 on: February 04, 2020, 09:21:26 AM »
Example of the virtual gridding method applied to:

z^7+A*z+c
c=( [6802992..6804016] + i*[6805040..6806064] ) * 2^-25
(the seed is a small complex interval only for test purposes)
A=(6805056+6806080*i) * 2^-25

The set was computed to level 15 without using a shrinkage value. The regions where positive values lie were computed at level 10 (colored small image below, colors are in heat-map fashion, highest shrinkage value was 15), and shriinkage values for gray pixels in those regions at level 15 were computed at level 16 in a 8x8 virtual grid before entering the main computation loop.

The result (large image) showed about 1.7% more interior cells (blue) and 3% more exterior cells (red).
(The image was downscaled 16-fold in a trustworthy manner, so a blue or red pixel here represents a 16x16 grid of blue or red pixels in the original size, hence smaller blue and red regions had to be colored back to gray)

Now I have to figure out when to apply the gridding method - at every refinement level or only once in a while?

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #127 on: February 09, 2020, 12:57:27 PM »
"Optimization pipeline: Forever double - 1st step"

Computing an image with the TSA using a datatype that does not provide sufficient precision (but is much faster than a high precision one) might result in an image that does not have mathematical guarantee.

"Might" because:
  • My theoretical analysis is very conservative. Every term adds to the final sum, subtractions are converted to additions, and only absolute values are used. So the bits required might actually be fewer than the calculated one.
  • The bbx's complex plane coordinates are not directly used to judge a cell's fate, but rather the set of pixels it intersects with (the "screenRect"). So even in case of a rounding error, one might still end up with the same set.

If I'm only interested in finding interior - which is often the case when trying to find cycles, there is a quick a posteriori check to validate an image for its mathematical guarantee: Every black cell's bbx is only allowed to intersect with black pixels - and this can be performed by a one-pass over the image's black pixels computing the screenRect with a sufficient datatype. If any one of them intersects with gray or white, the image has to be discarded.

So I'm currently detecting interior in some images using C++'s double datatype, no matter the actual bit requirements. And afterwards I perform the validity check. If it is passed, the image is guaranteed to hold the Figueiredo algorithm's claims: black is definite interior, white is definite exterior and the union of gray contains the Julia set.

However, I can only detect black with that method. So I'm computing the image with the correct datatype to some extent to approximate the Julia set exterior, then judge gray pixels only whether they have a path to white or not (even in case all paths lead to white, this information might be false due to rounding, so I will reduce it to "has a-path-to-white" which means the pixel stays gray forever at that refinement level - and being gray does never affect the mathematical guarantee.

Here's the first example at level 18. Up to level 17, double was sufficient, then my analysis required long double, but the image was judged further (for interior) with just double and afterwards the validity check was performed and passed.

z^3+A*z+c at level 18
c=(2052096-983040*i) * 2^-25
A=(-17310248+29024956*i) * 2^-25

I like the Figueiredo algorithm and what can be "squeezed out" of it in terms of optimiztation.


Offline Adam Majewski

  • *
  • Fractal Furball
  • ***
  • Posts: 226
« Reply #128 on: February 09, 2020, 03:52:39 PM »
Here's the first example at level 18. Up to level 17, double was sufficient, then my analysis required long double, but the image was judged further (for interior) with just double and afterwards the validity check was performed and passed.

z^3+A*z+c at level 18
c=(2052096-983040*i) * 2^-25
A=(-17310248+29024956*i) * 2^-25


Is it possible to remove one parameter ? ( then parameter space will be only 2d complex plane )

similar image :


https://commons.wikimedia.org/wiki/File:Julia_set_of_the_quadratic_polynomial_f(z)_%3D_z%5E2_-_1.12_%2B_0.222i.png

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #129 on: February 10, 2020, 09:15:00 AM »
@Adam: You're dancing rabbit set shows in deed a similarity to my last cubic set. Here's the level 20 version (I only checked for interior at level 20 and used the TSApredictor results, so black may have emerged earlier. Cycle detection was not performed (the reference point method would have worked, but taken probably several days - I'm currently implementing a faster version, so I'll come back to that).

z^2+(-37580964+7449084*i) * 2^25 approx. -1.120000005+0.222000003i

(My decimal-number routine from the parabolic thread cannot handle such depths for the actual parameters -1.12+0.222i, so I took the values from above (and hoped to still be inside the Mset, which, after finding black cells, was confirmed).

I initially chose the cubic one (brute-force found in my quest for cycles inscribed in the hull of other cycles) as it shows an interesting triangle-anchored period-39 cycle, which unfortunately was not detected at level 20 (the current maximum for most sets), so further optimizations are necessary there to "peek" a bit into higher refinement levels.

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #130 on: February 13, 2020, 09:08:15 PM »
"Optimization pipeline: Forever double - step 2"

As I learnt that today's CPUs's optimizations still result in deterministic outcome (forum link), I could now devise a new routine that works with double at all times (if a test is passed) and can be used to detect cycles.

As I mentioned in post #127, the bbx's coordinates are not the primary outcome of the algorithm, but rather the screenRect (intersecting pixels).

Test procedure
I compute all gray cells' bbx once with a sufficient datatype and once with plain double, then convert both to the screenRect. If those two are identical for all cells of interest, double can be used and the mathematical guarantee still holds.

I now perform for a given refinement level, formula and parameters this test before entering the main loop. Then the whole image is computed with double, and no validity check is necessary.

The cost for that optimization is minimal: One pass over the image with the slow datatype and one with the fast. This gets quickly amortized in the repeated loops with double instead of e.g. __float128 later on.

Here's an example Julia set (the same as in post #125):

z^4+A*z+c
c=(-13590528+12922880*i) * 2^-25
A=(2487868-17299244*i) * 2^-25

Bit requirements at level 19 were 76=8 (integer part) + 68 (fractional part), so double normally would not be an option.

At level 19, the immediate basins of two cycles are colored differently (cyan, period 6; purple, period 10), black is the union of the non-immediate basins of all cycles (see details). At this zoom level (32-fold downscaled), only 3 immediate basins of the turquois cycle are visible. The attached zip file contains a 10k version, where all 6 (red arrows in north part of image). Those now hold the record for smallest immediate basins found.

Technical details
  • Equality of screenrects for two number types is only checked if both bbxes completely lie in the R-square of interest. For bbx's completely outside, it is only demanded that both lie outside, but the actual position is not compared. If both lie partially outside, partially inside, only the inside part has to be identical (as both have a path to the exterior). Any other case results in failing the test.
  • I used a new version of the reference point-based periodicity detection method where I had to trade in the possibility to color non-immediate basins to gain a substantial speed advantage, hence the black regions.

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #131 on: February 18, 2020, 09:42:29 AM »
Cubic Julia set with a detected periodic-54 cycle. Double number type was validated to be mathematically safe at levels 19 and 20.

Interestingly the set has a second period-54 cycle, whose immediate basins are in between the red detected ones, however not emerging till level 20.

c=(-1228800-307200*i) * 2^-25
A=(-35378933+7505309*i) * 2^-25

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #132 on: February 20, 2020, 09:13:16 PM »
"Peek at level 21"

In line with all the optimizations of the TSA, I implemented a new version of the TSApredictor to take a peek into higher refinement levels.

The predictor computes (part of) an image of size 2^REFINEMENTLEVEL x 2^REFINEMENTLEVEL pixels: Around each poiint-sampled identified periodic point I set a pixel square of size 2048 x 2048 (called a "region"). The outside is defined as the outside of the union of all (disjoint) regions.

Each pixel in a region is analyzed whether it has a path to the outside and this information is backpropagated (as in the original article). If no further change occurs, the still gray pixels are bounded, land only in other regions and are contained in the interior of the Julia set.

The cubic set z^3+A*z+c

c=(9142272+1933312*i) * 2^-25
A=(-1307924-30649211*i) * 2^-25

has (point-sampled) a period-4 and period-28 cycle. The former being readily emergent at level 12, the latter not at level 17, and the predictor shows at level 21 interior for all 32 periodic points' regions (calculation time only 6 min with float128 using precomputed and stored screenRects).

The image below shows TSA level 17, the period-4 cycle (green) and the periodic points (point-sampled identified, red squares, afterwards copied into image).

The right part of the image shows the 32 periodic points' regions (in no particular order). Red are pixels having a path to the outside of all regions, black are pixels that land only in black of another region and are therefore bounded.

Next step is to introduce a cycle detection method in the predictor.

Offline marcm200

  • *
  • Fractal Frankfurter
  • *
  • Posts: 565
« Reply #133 on: Yesterday at 10:12:01 AM »
"All in"

Quartic Julia set z^4+A*z+c with

c=(-28418704-3814114*i) * 2^-25
A=(-15226452+4431922*i) * 2^-25

showing 3 cycles of length 2 (immediate basins in green), 4 (purple), 21 (red, the wheel-like structure between the two green large Fatou components) at level 20 in the 2-square. The image is downscaled 128-fold, so it might have happened that no longer all 21 small red basins are visible (I have not counted them).

I had to go all-in with that image: Levels 18-20 were validated for the usage of double instead of a higher precision number type, almost all of the recent optimizations were in place.

Technical details
  • The new TSApredictor was used to check whether the 21-cycle could be detected in level 20 (using a pixel-rectangle of size 5120x5120 around every periodic point was necessary as opposed to the standard 2k width).
  • As memory was not a problem until level 20, in earlier levels I precomputed whole screenRects and stored those for reusage.
  • Cycle detection was performed with sufficient fpa datatype (fixed point, 32 bit integer part, 96 bits fractional part) as there was already interier present in the level 20 start and number type verification only checks gray cells.


clip
3d Julia sets: True shape

Started by marcm200 on Image Threads

31 Replies
896 Views
Last post January 20, 2020, 11:42:23 AM
by marcm200
clip
True-shape based C++ oracle for the int-/exterior of Julia sets

Started by marcm200 on Programming

6 Replies
280 Views
Last post August 23, 2019, 10:47:10 AM
by marcm200
clip
3D M-set for tricomplex numbers: Escape time and shape change?

Started by marcm200 on Fractal Mathematics And New Theories

2 Replies
230 Views
Last post May 04, 2019, 09:58:37 PM
by marcm200
xx
True shape/IA 2D/3D Julia set source code

Started by marcm200 on Programming

12 Replies
534 Views
Last post January 21, 2020, 01:35:39 PM
by marcm200
clip
Escape Time Blur

Started by AlexH on Share a fractal

0 Replies
389 Views
Last post August 15, 2018, 05:42:06 AM
by AlexH