I added a "disable texture resize" option to KF, which means you can access actual texture image pixels in OpenGL, which makes the "legendary colour palette" technique much easier (no need to convert the image to a linear palette in text format first).

colouring algorithm code:

```

vec3 colour()

{

if (getInterior()) return KFP_InteriorColor;

int w = textureSize(KFP_Texture, 0).x;

int h = textureSize(KFP_Texture, 0).y;

int N = to_int(floor(div(getN(), KFP_IterDiv)));

N += int(KFP_ColorOffset) * w * h / 1024;

int i = N % w;

int j = (N / w) % h;

return texelFetch(KFP_Texture, ivec2(i, j), 0).rgb;

}

```

(if you try this in latest released KF you will probably be disappointed, as the texture image is resized to match the output image size)

Rodney (selection) https://diode.zone/videos/watch/350b50b9-09f4-483d-929a-ff00aedad309

@celesteh https://kfjc.org/listen/archives has archives of both programs for the next week or so (search in the page for Cousin Mary)

fractal perturbation maths summary

I finally ixed the bug that made OpenCL crash when run in a background thread. Typically, it had nothing to do with OpenCL: when the user interface is wiggled during rendering it sets a "stop" flag to true which the various parts of the renderer inspect often to know their work is no longer needed, but the OpenCL part didn't notice, so after the reference and series approximation calculations had given up the internal state wasn't quite pristine.

The symptoms were crashing on an assertion failure (series approximation iteration count was not less than total number of iterations), but the fix was a few words, essentially "if (stop) return;" at the start of the OpenCL rendering function.

With OpenCL in a background thread, it actually becomes useable for exploration - previously the user interface blocked until rendering was complete (after all references).

fractal perturbation maths summary

benchmarks for one Burning Ship showed:

rescaled double was fastest on CPU; location was not "deep needle" relative to double's range

extended float was fastest on GPU; location was "deep needle" relative to float's range so rescaled float performed poorly. rescaled double was a close second

rescaled double was fastest on GPU for hybrids

hybrids on CPU are super slow (they are "interpreted", there are many branches on the (constant) formula parameters in the inner loop; when generating OpenCL source code at runtime all these branches can be removed as it is "compiled").

fractal perturbation maths summary

concretely, those benchmarks showed:

6.3s wall, 45.2s CPU for rescaled double with 3 series terms

10.0s wall, 21.1s CPU for rescaled double with 10 series terms

10.5s wall, 6.4s CPU for rescaled float hybrid on GPU with OpenCL

fractal perturbation maths summary

Earlier I implemented rescaled iterations for hybrids with and without derivatives for OpenCL.

Meanwhile I discovered series approximation breaks if the reference orbit is stored in single precision; but later rounding the coefficients calculated with a double precision orbit to single precision seems to work ok. Yesterday I hacked the reference storage abstraction to store double in addition to single if single is requested.

Today I added detailed timing information collection. It reports wall clock time and CPU time for reference calculations, series approximation calculations, and perturbation calculations, as well as a total (which is more than the sum of parts because it include loading/colouring/saving).

Benchmarks on one location were revealing:

- increasing the number of series approximation terms is good for efficiency (lower CPU time) but bad for racing (higher wall clock time)

- series approximation takes a surprisingly large amount of time

- rescaled single precision float hybrids on GPU using OpenCL are competitive speed-wise with built in Mandelbrot on CPU (without OpenCL), even though neither threaded reference calculations nor series approximation can be used.

- hybrids without OpenCL are too slow and inefficient to be used at all

fractal perturbation maths summary

A preliminary benchmark at one location (zoom 10^533) with derivatives:

scaled double: 3m50s

long double: 11m51s

So scaled double is a big win.

As it's just a Mandelbrot set (just configured with the hybrid formula editor) it can also be rendered with a hardcoded algorithm:

scaled double: 18.5s

long double: 54.1s

So the hardcoded algorithm is very much faster (it can use complex differentiability instead of 2x2 Jacobian, plus no branches or loops in the iteration step.

Complex differentiability also means that series approximation works, which speeds it up a bunch even more:

scaled double: 10.4s

long double: 19.4s

So altogether the hybrid is 24x slower than the built in thing for the same image. I think OpenCL might narrow the gap, so I'll implement that next.

fractal perturbation maths summary

For distance estimation of hybrid formulas I use dual numbers for automatic differentiation.

One small adjustment was needed for it to work with rescaled iterations: instead of initializing the dual parts (before iteration) with 1 and scaling by the pixel spacing at the end, initialize the dual parts with the pixel spacing and don't scale at the end. This avoids overflow of the derivative, and the same rescaling factor can be used for regular and dual parts.

I did try something else with rescaling the regular and dual parts with different factors, but didn't get it working.

C++ programming (-)

I had T operator*(int, T) and T operator *(T, int) and T operator*(double, T) and NO T operator*(T, double), so multiplying would work fine in one order and the other order would silently truncate my double to int without any error or even warning. Cue hilarious bugs.

Yes, g++ -Wconversion exists, but the vast flood of warnings on this large legacy code base would be just as useless.

Code is fixed now, but that's some hours I won't get back.

fractal perturbation maths summary

Rescaling hybrid iterations seems like a big challenge, but might not be so hard after all (I've not implemented it yet but this is the algorithm/maths that I've worked out so far):

if either or both the real and imaginary parts of the reference orbit Z are small, one needs to do a full range iteration with floatexp and recalculate the scale factor afterwards, as with refgular formulas like Burning Ship.

otherwise, thread s through from the top level down to the operators.

initialize with

W := Z + z*s

if absolute x

re(z) := diffabs(re(Z/s), re(z))

similarly for aboslute y with im

When combining operators (this subterm only occurs with multiplication) replace

f(op1, Z + z) with op1(Z + z*s)

and that's all the changes that need to be made as far as I can tell!

I thought it'd be much harder.

Let's see what goes wrong when I try to implement it tomorrow...

fractal perturbation maths summary

Hybrid fractals in KF are built from stanzas, each has some lines, each line has two operators, and each operator has controls for absolute x, absolute y, negate x, negate y, integer power, complex multiplier. The two operators in a line can be combined by addition, subtraction or multiplication, and currently the number of lines in a stanza can be either 1 or 2 and there can be 1, 2, 3 or 4 stanzas. The output of each line is fed into the next, and at the end of each stanza the +c part of the formula happens. There are controls to choose how many times to repeat each stanza, and which stanza to continue from after reaching the end.

Implementing perturbation for this is quite methodical. Start from an operator, with input Z and z being the reference and delta.

Set mutable variables

z (as input)

W := Z + z

B := Z

If absolute x, then

re(z) := diffabs(re(Z), re(z))

re(W) := abs(W)

re(B) := abs(B)

similarly for y and im

If negate x, then

re(z) := -re(z)

etc for W and B

similarly for y and im

Now compute

S = sum_{i=0}^{power-1} W^i B^{power-1 - i}

and return a*z*S

Combining operators into lines may be done by https://mathr.co.uk/blog/2018-03-12_perturbation_algebra.html

Combining lines into stanzas can be done by iterating unperturbed Z alongsize perturbed Z, z; only the +C needs high precision.

fractal perturbation maths summary

ugh, left image is incorrect, right image is correct, at least the image descriptions are ok

fractal perturbation maths summary

Here's an example of the kind of glitches that can happen (in this case, missing details in the right hand image with KF's default threshold), with a lower tolerance for glitches it takes much longer but the image is hopefully correct.

Test location in the Mandelbrot set taken from here https://www.deviantart.com/olbaid-st/art/Deep-Mandelbrot-Set-023-10-1086-695492737

Show thread

fractal perturbation maths summary

About glitches: perturbation assumes exact maths, but some images have glitches when naively using perturbation in low precision.

A solution was discovered/invented by Pauldelbrot in http://www.fractalforums.com/announcements-and-news/pertubation-theory-glitches-improvement/msg73027/#msg73027

One has perturbed iteration as in [3] (recap: z -> 2 Z z + z^2 + c)

Then one perturbs this with z -> z + e, c -> c + f:

[13]: e -> (2 (Z + z) + e) e + f

We are interested what happens to the ratio e/z under iteration, so rewrite [3] as

[14]: z -> (2 Z + z) z + c

Pattern matching, the interesting part (assuming c and f are small) of e/z is 2(Z + z) / 2 Z.

When e/z is small, the nearby pixels "stick together" and there is not enough precision in the number type to distinguish them, which makes a glitch.

So a glitch can be detected when

[15]: |Z + z|^2 < G |Z|^2

Where G is a threshold (somewhere between 1e-2 and 1e-8, depending how strict you want to be).

This does not add much cost, as |Z+z|^2 already needs to be computed for escape test, and G|Z^2| can be computed once for each iteration of the reference orbit and stored.

The problem now is: How to choose G? Too big and it takes forever as glitches are detected all over, too small and some glitches can be missed leading to bad images.

The glitched pixels can be recalculated with a more appropriate reference point.

@jayrope compare my https://code.mathr.co.uk/rodney/blob/6faf3d9ebdb043c79ce38d90179fe8517bf6518a:/sonify.cl#l28 vcf function with the inner loop of Pd's https://github.com/pure-data/pure-data/blob/eeef4ba9130d3182146927c37fa57d61bbff0f0b/src/d_osc.c#L386-L447 (Pd has a linear interpolated table lookup for cos() which is faster but less accurate; I just use the C math library which makes it simpler but slower - lines 415,419-433 of Pd's code are implementing something like `cos(cf)` and `sin(cf)`)

the process is extracting the core mathematics/algorithms, and hardcoding it in C without any of the wrapping API layers that make it useable from a patcher system

https://code.mathr.co.uk/rodney/tree/6faf3d9ebdb043c79ce38d90179fe8517bf6518a check sonify.cl for the DSP, sonify.c for the control code, and sonify.sh for how to use it

Just saw this in my inbox from KFJC radio station:

---8<---

KFJC Month of Mayhem 2021

Music From Mills, Part 1

Saturday Mayhem 1 6PM-8PM

Music from Mills Part 2

Saturday Mayhem 8 at 6PM-8PM

Hosted by Cousin Mary

Recent news of the scheduled closure of Mills College has hit the Bay Area experimental music community especially hard. Join Cousin Mary in an exploration of the music associated with Mills College graduates and faculty.

---8<---

making art with maths and algorithms

Joined May 2018