automatic skew adjustment for Burning Ship etc

  • 18 Replies
  • 630 Views

0 Members and 1 Guest are viewing this topic.

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1063
    • mathr.co.uk
« on: August 26, 2018, 11:07:57 PM »
I think I figured out how to automatically unskew the Burning Ship:

Starting from the size estimate for Mandelbrot set mini-sets:
https://mathr.co.uk/mandelbrot/book-draft/#size-estimate
Make L and B be 2x2 matrices instead of complex numbers. Note that (for the Burning Ship power 2 at least) L will not have any non-conformal stretching, so at the end you just need to consider B instead of (L2 B).  Other fractals might need diagonalisation to handle \( L^\frac{d}{d-1}B \) Then the rotation and skew adjustment matrix (to be applied to delta-C coordinates) is:

\[ \frac{B}{\sqrt{|\det B|}}^{-1} \]

A future version of KF might include this as part of the Newton-Raphson zooming (as it is mini-set dependent).  It would need a change to the .kfr file format to store a 2x2 rotation and stretch matrix (with determinant 1), instead of the current rotation and skew ratio.
« Last Edit: August 27, 2018, 04:17:09 PM by claude, Reason: markup »

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1063
    • mathr.co.uk
« Reply #1 on: August 27, 2018, 02:46:22 PM »
Maybe I was wrong about L not being non-conformal, these are rendered using L2B to unskew.  Higher powers might be problematic with multivalued rational matrix powers?

image 1
Code: [Select]
re"-1.7376030903855331822680777237384033625"
im"-1.1874559032877022930121132805164953765e-2"
size"1.8313472770622653e-21"
("transform",(3.623779526375508,-1.0143835039976357,-11.626049226570013,3.5303672475055015))
stretch =  12.719
stretchAngle =  2.8487
rotation = -0.97761

image 2
Code: [Select]
re"0.8741850936284249581163187788301509472"
im"-1.514132440804810318582352545762872942"
size"1.5340280079199914e-20"
("transform",(3.8650265392253633,1.796297826808982,1.6408887811634636,1.021344850178012))
stretch =  4.6749
stretchAngle =  0.45570
rotation = -0.031794

image 3
Code: [Select]
re"-0.980205145359994589554021418837164241"
im"-0.87926047680022105334908200540068974"
size"4.1432554188574579e-19"
("transform",(1.9402951718628327,0.26773345041536056,-1.1286194048356972,0.35965189351455074))
stretch =  2.2453
stretchA =  0.023460
rotation = -0.54564
« Last Edit: August 27, 2018, 02:57:43 PM by claude, Reason: markup »

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1063
    • mathr.co.uk
« Reply #2 on: August 27, 2018, 04:16:05 PM »
Cubic Burning Ship, auto-unskewed using just B, seems to work in these tests.  Using L^(3/2)B might be more proper (?), but in earlier tests I got into all sorts of bother with complex eigenvalues when trying to diagonalize L to do the non-natural power.  I was too lazy to paste the transform matrices into my Octave polar decomposition code to get the rotation and stretch parameters this time, the code is at https://mathr.co.uk/blog/2018-03-06_polar_decomposition_in_2d.html

Image 1
Code: [Select]
re"-1.331545877755506536679413667"
im"-1.866738708385317784521630591e-2"
size"5.6873069705376381e-11"
("transform",(1.9269337597628975,0.9630731164146824,0.27956877402716035,0.6586864566692452))

Image 2
Code: [Select]
re"0.69915388036874369446264475638759988"
im"-1.0112702532669882764602180207314239"
size"5.7138649328121838e-18"
("transform",(2.2936795502913774,-2.8047825330222804,-1.647246495315255,2.4502848259376235))

Image 3
Code: [Select]
re"-0.655414004008662016315142684293653933"
im"-0.66290269845631425263492772004157505"
size"9.6313199324963481e-20"
("transform",(0.9311451388482219,0.37639007954903514,-0.31203652645385366,0.9478141593237559))

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1063
    • mathr.co.uk
« Reply #3 on: August 27, 2018, 04:28:07 PM »
Mandelbar (power 2) unskewed with B:

Image 1
Code: [Select]
re"0.409949752902015175924"
im"-1.13229287695124820074"
size"7.6341550596726096e-5"
("transform",(0.8222675871469319,-0.6313079716166592,-9.501650393705972e-2,1.289099428141756))

Image2
Code: [Select]
re"0.76163946268372037852553351"
im"-1.1760393075628667190689237"
size"2.6795187395230754e-9"
("transform",(0.3385326994408471,0.42189527498794765,0.6060653725043186,3.7092313949801894))

Image 3
Code: [Select]
re"0.761911347420476081235504"
im"-1.176290433766432477647868"
size"7.9635497079833018e-9"
("transform",(4.469684871193343,-2.55715806999073,-1.7282119877629967,1.21245935392219))

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1632
« Reply #4 on: August 28, 2018, 04:37:31 AM »
Looks perfect. Not sure if it can be done in a rigorously correct way. I jotted down some formulas below.

Let the period p orbit at supercritical point \( c_0 \) be \( z_0, \ldots, z_p \) with \( z_0=z_p=0 \).
(z is a 2-vector.) Iteration is \( z \leftarrow f(z) + c \) and Jacobian is \( J_k=\frac{\partial f}{\partial z}(z_k) \).

I get for the linearized p-times iteration
\( z \leftarrow L(f(z)+B(c-c0)) \) where:
\( q_m = \prod_{k=1}^m J_k^{-1}\\
L = (q_{p-1})^{-1}\\
B = 1 + \sum_{k=1}^{p-1} q_k.
 \)

Problem is the L factor can not be generally absorbed in \( f(z) \) with a change of variables \( z \rightarrow A z \) for some matrix A, unless f(z) has a special nice form. Maybe some heuristic like considering the asymptotic behavior (large z) will work generally. It seems from the results that the method is quite robust against "not doing it rigorously" which is good news.

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1063
    • mathr.co.uk
« Reply #5 on: August 28, 2018, 10:14:24 AM »
I think the asymptotic behaviour for small Z is what matters for renormalization (not large Z).  In 'et' I compute this by simply trying a small Z, and assuming the asymptotic behaviour is a small positive integer power.  Not great, but seems to work for most of the formulas I've tried.

I added the auto-skew as an option to KF's NR-zoom dialog, coming in the next release (probably later today).

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1632
« Reply #6 on: September 05, 2018, 04:24:14 AM »
The auto-skew works great if you can find a mini, but there are skewed regions where a mini can't be found.
I haven't thought this trough, so ignore following if nonsense.

Consider the attachment of a burning ship location and corresponding Julia set. Parameter space needs skewing, but Julia does not.
As the z-iterates around a tiny region around a point c0 are related to the z iterates around a tiny region in the center of the corresponding Julia set
(by something like \( \frac{\partial z}{\partial z_0}/\frac{\partial z}{\partial c} \) I think) would it be possible to use that ratio (matrices of course) to unskew if no mini is available?

Third image is manually unskewed BS parameter space.

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1632
« Reply #7 on: September 05, 2018, 07:12:39 PM »
I tried it and it works very well though I just tried on Mandelbar.
If c0 is your escaping reference location, compute orbit and matrices dfdc and dfdz (z = (x0,y0)).
Skew matrix is then S = inv(dfdc)*dfdz, S -> S/sqrt(abs(det(S))) evaluated at escape iteration.
Then c <-- S*(c-c0)+c0.
Example below.
PS Seems to work also for non-escaping points.
« Last Edit: September 05, 2018, 08:01:10 PM by gerrit »

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1063
    • mathr.co.uk
« Reply #8 on: September 05, 2018, 08:18:45 PM »
Great research!  I'll see if I can add it for all formulas to the next release of KF.

Here is how I convert the matrix S (with det(S)=1) into KF's rotation and skew ratio - unfortunately KF doesn't support rotation after skew (only before) so if you need skew then the rotation angle is fixed by the skew direction.  And skew is represented strangely in KF - it only stretches horizontally, and 360 is unskewed (stretch factor 1.0).  The key algorithm is "polar decomposition", https://mathr.co.uk/blog/2018-03-06_polar_decomposition_in_2d.html

https://code.mathr.co.uk/kalles-fraktaler-2/blob/kf-2.13.8:/common/matrix.cpp#l37 polar_decomposition
https://code.mathr.co.uk/kalles-fraktaler-2/blob/kf-2.13.8:/fraktal_sft/fraktal_sft.h#l478 SetTransformPolar
https://code.mathr.co.uk/kalles-fraktaler-2/blob/kf-2.13.8:/fraktal_sft/fraktal_sft.cpp#l2815 SetRatio

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1632
« Reply #9 on: September 06, 2018, 01:16:08 AM »
Works fine for burning ship too. Here's a location on the Eastern shore.

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1063
    • mathr.co.uk
« Reply #10 on: September 06, 2018, 04:46:11 AM »
tested gerrit's method in 'et' and it works* for:
- Burning Ship power 2 : (|x|+i|y|)^2+c
- Burning Ship power 3 : (|x|+i|y|)^3+c
- Mandelbar power 2 : (x-iy)^2+c
- Mandelbar power 3 : (x-iy)^3+c
- Perpendicular Buffalo power 2 : w := (x - i |y|)^2 ; z := (|u| + i v) + c
- Perpendicular Buffalo power 4 : w := (x - i |y|)^4 ; z := (|u| + i v) + c  (not sure if this formula is in KF, maybe under a different name?)

tested and fails (makes it worse instead of better) for
- SimonBrot 4th :  z := z^2 (|x| + i |y|)^2 + c
- SimonBrot2 4th : w := z^2 ; z := w (|u| + i |v|) + c
- 4th Celtic False Quasi Heart :  z := (|x^4 + y^4 - 6 x^2 y2| + i 4 x y |x^2 - y^2|) + c

sometimes I get huge skew matrix elements, or the rotation is not accurate enough, or issues with reflection? (det(S) comes out as -1 instead of +1, not sure if anything different is needed to handle this case properly, anyway some failure cases don't have negative det(S))
one example: det(S) = (-6.8951363957169e7) * (-6.3092208995121166e7) - 2.824506850185153e8 * 1.5401958982677631e7 = -1.0

maybe something is wrong with my calculation of S? perhaps one or both matrix is transposed?  I always get confused which way Jacobian matrices should go... eyes appreciated at my code here:
https://code.mathr.co.uk/et/blob/refs/heads/master:/src/Fractal/EscapeTime/Algorithms/R2/Skew.hs#l135


* if you try and use this skew algorithm in the exact center of a mini, it can break - I think this is because z hits 0 0 and you end up with a 0 0 0 0 matrix with all the information lost... maybe it would work to bail out before this happens?

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1632
« Reply #11 on: September 06, 2018, 06:22:20 AM »
Tried Simonbrot4. Numerically I find dfdz to be pathologically ill-conditioned (also dfdc). If do a few iterations it's singular. (Rounding errors must make it non-singular.)

Jacobian fp of the map seems to be the problem:
Code: [Select]
    fp(1,1) = 2*(x^2-y^2)*2*x-4*y*abs(x*y)-4*x*y*abs(y)*sign(x);
    fp(1,2) = 2*(x^2-y^2)*2*y-4*x*abs(x*y)-4*x*y*abs(x)*sign(y);
    fp(2,1) = 4*x*(x*y+abs(x*y)) + 2*(x^2-y^2)*(y+abs(y)*sign(x));
    fp(2,2) = -4*y*(x*y+abs(x*y)) + 2*(x^2-y^2)*(x+abs(x)*sign(y));

You can see if x and y have opposite sign the lower row is zero, hence fp is singular.

Offline claude

  • *
  • 3e
  • *****
  • Posts: 1063
    • mathr.co.uk
« Reply #12 on: September 06, 2018, 07:52:23 AM »
Ouch.  Is there maybe some way to rescue it?  Like if A and B are singular, so C = A-1 B is a meaningless calculation, but maybe a semantically relevant C could be calculated in a different way (similar to analytic continuation)?  Perhaps something like L'Hpital's Rule for indeterminant limit ratios of real numbers, extended to matrices? https://en.wikipedia.org/wiki/L%27H%C3%B4pital%27s_rule

EDIT found https://arxiv.org/abs/1209.0363v1"A L'hospital's rule for multivariable functions" by Gary R. Lawlor 2012
Quote
We prove a version of L'hospital's rule for multivariable functions, which holds for non-isolated singular points. We also give an algorithm for resolving many indeterminate limits with isolated singular points.
However, it's analysis/calculus stuff, not directly implementable constructive algorithms.  I don't know enough to even know if it's relevant.
« Last Edit: September 06, 2018, 09:07:30 AM by claude, Reason: paper »

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1632
« Reply #13 on: September 06, 2018, 05:08:44 PM »
Problem is in the II and IV quadrant the function collapses to a real function \( \pm |z|^4 \) so it's really pathological, so the whole idea of local Julia/parameter space patching collapses.

Offline gerrit

  • *
  • 3f
  • ******
  • Posts: 1632
« Reply #14 on: September 08, 2018, 12:16:56 AM »
Been messing around with this some more.

Using S = inv(df/dc)*df/dz (normalized to have |det| 1) makes the parameter space image locally similar to the Julia set, based on observation Julia set is not stretched.

Why not just use S = inv(df/dc), which makes the par. space locally isotropic, independent of Julia set?

It seems to make no visual difference, but the problem of "df/dz=0" when your image center is at nucleus center is gone and it's of course simpler.

Note added: I did find some cases where not including df/dz fails (\( z \leftarrow (x^2-y^2)+i(x^2+2xy) \)). See attached "test".
Also found an example of the opposite  (\( z \leftarrow (x^2-y^2)+i(x^2-2xy) \)) but the skew is rather trivial there. Also attached "test2" (except fail image with Julia which goes over #attach limit, it's just totally black).
« Last Edit: September 08, 2018, 02:10:30 AM by gerrit »


clip
el day de lay (burning ship)

Started by claude on Fractal movie gallery

0 Replies
134 Views
Last post April 30, 2018, 09:12:50 PM
by claude
xx
Burning Ship

Started by buddhi on Fractal Image Gallery

5 Replies
267 Views
Last post January 24, 2019, 03:14:44 PM
by gerson
clip
Burning Ship

Started by Bill Snowzell on Fractal Image Gallery

2 Replies
200 Views
Last post February 22, 2018, 03:38:51 PM
by Bill Snowzell
clip
Burning Ship

Started by Bill Snowzell on Fractal Image Gallery

5 Replies
323 Views
Last post October 02, 2017, 06:23:14 PM
by Sabine62
xx
Burning Ship

Started by Bill Snowzell on Fractal Image Gallery

3 Replies
262 Views
Last post September 13, 2017, 11:52:54 PM
by Fraktalist