Precision tests

Makefile

@@ -0,0 +1,54 @@
+
+CFLAGS=-Wall 
+#INC=-I/usr/local/include/flint
+INC=
+#LIBS=-lm -lflint
+
+all: arb_pythonlib mpfr_pythonlib
+
+# In python, "import arb_fractalmath"
+# arb_fractal_lib.c
+# arb_fractal_lib.h
+# (compiles into)
+# libarbfractalmath.a
+#
+# arb_fractalmath.pyx
+# (makes)
+# arb_fractalmath.c
+# arb_fractalmath.so
+
+# In python, "import mpfr_fractalmath"
+# mpfr_fractal_lib.c
+# mpfr_fractal_lib.h
+# (compiles into)
+# libmpfrfractalmath.a
+#
+# mpfr_fractalmath.pyx
+# (makes)
+# mpfr_fractalmath.c
+# mpfr_fractalmath.so
+
+arb_pythonlib: libarbfractalmath.a arb_fractalmath.pyx
+	python3.9 setup_arb.py build_ext --inplace
+
+libarbfractalmath.a: arb_fractal_lib.o
+	ar -ru $@ $^ $(LIBS)
+	ranlib $@
+#	$(CC) -shared $(LDFLAGS) -o $@ $^ $(LIBS)
+
+mpfr_pythonlib: libmpfrfractalmath.a mpfr_fractalmath.pyx
+	python3.9 setup_mpfr.py build_ext --inplace
+
+libmpfrfractalmath.a: mpfr_fractal_lib.o
+	ar -ru $@ $^ $(LIBS)
+	ranlib $@
+
+%.o: %.c
+	$(CC) $(INC) $(CFLAGS) -c -o $@ $<
+
+clean:
+	python3.9 setup_arb.py clean --all
+	python3.9 setup_mpfr.py clean --all
+	rm -f *.o *.d *~ *.so *.a 
+	rm -rf build
+

README.md

@@ -1 +1,112 @@
 # mandl
+
+
+# Constructing A Dive
+
+## Exploration
+
+Establish at least two frames to dive between.
+ - real and imaginary starting widths
+    (real -> horiz, imag -> vert, or was it the other way?)
+ - real and imaginary ending widths 
+    (second width is redundant if aspect ratio is constant)
+ - single center point, used as middle of both start and end frames
+
+While exploring, write out points of interest, creating a waypoint list.
+The first and last waypoints are the default endpoints for an 
+editing project.
+
+## Editing
+
+Need at least 2 waypoints to establish the 'outer' parameters for
+a dive animation.  
+
+In terms of an epoch-based definition, these are the starting widths 
+and the appropriate combination of zoom factor, framerate,
+and duration, that achieve the window widths at the start and the end.
+
+In terms of mesh exploration, these are two meshes that define the starting
+and ending widths, and either a framerate or a duration.
+
+## Statistics
+
+For an edit, would like to have statistics gathered across the frames
+that show:
+ - How many pixels are which values (counts and/or hists)
+ - Some kinds of entropy estimates for each frame
+ - Pattern matching results, against a library of shapes
+
+## Sync
+
+Timing of color palette and drawing algorithm type should be flexible
+enough that we can sync it to a music track, or to our defined waypoints.
+
+Sound generations from a sequence might require timing adjustments for
+the animation, to make a rhythm or a tone adjustment more consistent.
+e.g. if we have a good start of a beat, but the 4th bar is messy, we could
+rush or delay that phrase to make it fit more consistently.
+
+Sound sync will probably require a similar kind of adjustment, except the
+satisfaction is backwards, where we take a rhythm or a tone definition,
+and bend the temporal locations that we've targeted as waypoints to
+line up with the beat or tone.
+
+
+
+# Notes on Architecture and Classes
+
+FractalContext
+ - Keeps and manages overall parameters.
+ - Instantiates algorithm and precision-aware subclasses.
+ - Holds the Timeline for the animation.
+ - Responsible for generating images for frame numbers.
+
+
+MediaView
+ - Configures encoding parameters for output rendering.
+
+
+DiveTimeline
+ - Keeps framerate and image/mesh size.
+ - Keeps a sequential list of DiveTimelineSpans, which define animation parameters through their keyframes.
+ - Instantiates an appropriate Algo for every frame requested.
+
+
+Algo
+EscapeAlgo(Algo) and EscapeFrameInfo(FrameInfo)
+JuliaAlgo(EscapeAlgo) and JuliaFrameInfo(FrameInfo)
+ - Chaperone for algorithm-specific intermediates.
+ - Holds both algorithm invocations, and algorithm-specific pre and post processing hooks.
+ - Responsible for intermediates caching.
+
+
+MathSupport
+ - Implementations of library-specific calculations.
+ - Interpolators and multipliers for different numeric types.
+ - Native Python implementation in the base class.
+
+
+# Optimization Discussion
+
+## Inner loop multiplications optimizations at high bit depths
+
+Text taken from:
+https://randomascii.wordpress.com/2011/08/13/faster-fractals-through-algebra/
+
+When using floating-point math the speed of addition, subtraction, and multiplication are generally identical. However floating-point math has only about 52 bits of precision which is woefully inadequate for serious fractal exploration. All of my fractal programs have featured high-precision math routines to allow virtually unlimited zooming. Fractal eXtreme supports up to 7,680 bits of precision. While this isn’t truly unlimited it is close enough because fractal calculation time generally increases as the cube of the zoom level, and even on the fastest PCs available it exceeds most people’s patience between 500 and 2,000 bits of precision.
+
+The speed of squaring and multiplication is critical because in high-precision math they are O(n^2) operations and everything else is O(n). As ‘n’ (the number of 32-bit or 64-bit words) increases, the time to do the multiplications becomes the only thing that matters. Okay, at extremely high precisions multiplication doesn’t have to be O(n^2), but the alternatives are a lot of work for uncertain gain, so I’m ignoring them.
+
+While multiplication and squaring are both O(n^2) it turns out that squaring can be done roughly twice as fast as multiplication. Half of the partial results when squaring are needed twice, but only need to be calculated once.
+
+All the information necessary to find the algebraic optimization is now available.
+
+The observation (which came from somebody I know only through e-mail) was that the Mandelbrot calculations can be done using three squaring operations rather than two squares and a multiply. At the limit this gives approximately a one-third speedup!
+
+The one multiply was being used to calculate zr * zi. If we instead calculate (zr + zi)^2 then we get zr^2 + 2*zr*zi + zi^2. Since we’ve already calculated zr^2 and zi^2 we can just subtract them off and, voila, multiplication through squaring:
+
+
+====
+
+
+

algo.py

@@ -0,0 +1,224 @@
+
+
+import os
+import pickle
+
+import fractalcache as fc
+import fractalpalette as fp
+import divemesh as mesh
+
+from PIL import Image, ImageDraw, ImageFont # EscapeAlgo uses for burn-in
+
+# Should probably import Timeline, but not technically needed, since
+# it's cyclic?
+
+class Algo(object):
+    """
+    Algo is responsible for generating a single, entire, image 
+    frame from a specified mesh.
+
+    This means the algo is responsible for generating and caching its own
+    intermediate results.
+
+    The sequence in run() is the core of the processing architecture.
+
+    Trying something out: presuming that Algos will all write out 
+    intermediate files, without extra parameters to turn that behavior
+    off.  This is only really visible in this base class in mesh_setup(), 
+    in which you can't turn off the mesh pickle file writing.
+
+    In general, because this is a customizable sequence processor,
+    subclasses should NOT call superclass implementations, because they're
+    probably providing more specific behavior.
+    """
+
+    @staticmethod
+    def options_list():
+        """ 
+        Algorithm implementations can fill a dictionary with key/value
+        pairs to pass algorithm-specific parameters from the command line
+        into the per-frame algorithm instantiation
+
+        To do this, the list of available options are returned here, and
+        then iteration over the getopts values is done in parse_options()
+        """
+        return []    
+
+    @staticmethod
+    def load_options_with_math_support(opts, math_support):
+        return {}
+
+    def __init__(self, dive_mesh, frame_number, output_folder_name, extra_params={}):
+        self.algorithm_name = None 
+
+        self.dive_mesh = dive_mesh
+        self.frame_number = frame_number
+        self.output_folder_name = output_folder_name
+
+        # Algo should fill in these values
+        self.mesh_array = None
+        self.mesh_real_array = None
+        self.mesh_imag_array = None
+        self.counts_array = None
+        self.last_values_array = None
+        self.last_values_real_array = None
+        self.last_values_imag_array = None
+        self.processed_array = None
+        self.output_image_file_name = None
+
+    def get_metadata(self):
+        return {'frame_number' : self.frame_number,
+            'fractal_type': self.algorithm_name,
+            'precision_type': self.dive_mesh.mathSupport.precisionType,
+            'mesh_center': str(self.dive_mesh.getCenter()),
+            'complex_real_width' : str(self.dive_mesh.realMeshGenerator.baseWidth),
+            'complex_imag_width' : str(self.dive_mesh.imagMeshGenerator.baseWidth),
+            'mesh_is_uniform' : str(self.dive_mesh.isUniform())}
+
+    def run(self):
+        """ 
+        Frame generation happens in phases, with intermediate hooks available.
+        """
+        self.beginning_hook()
+        self.mesh_setup()
+        self.generate_counts()
+        self.pre_process_hook()
+        self.process_counts()
+        self.pre_image_hook()
+        self.generate_image()
+        self.ending_hook()
+        return self.output_image_file_name 
+
+    def beginning_hook(self):
+        pass
+
+    def mesh_setup(self):
+        """ 
+        Might not be important to split mesh array creation into
+        its owns step, but it does add flexibility to Algo subclasses
+        to handle this differently, or maybe to not save results out
+        to file.
+        """
+        # Originally relied on the directory structure to exist, but after
+        # clearing intermediates a few times, it got annoying, so let's just
+        # cross our fingers that this all ends up where we hoped it would.
+        if not os.path.exists(self.output_folder_name):
+            os.makedirs(output_folder_name)
+
+        mesh_base_name = u"%d.mesh.pik" % self.frame_number
+        mesh_file_name = os.path.join(self.output_folder_name, mesh_base_name)
+        with open(mesh_file_name, 'wb') as mesh_handle:
+            pickle.dump(self.dive_mesh, mesh_handle)
+
+        self.mesh_real_array = self.dive_mesh.generateRealMesh()
+        self.mesh_imag_array = self.dive_mesh.generateImagMesh()
+
+    def generate_counts(self):
+        """ Business end of getting results for a mesh """
+        raise NotImplementedError('Subclass must implement generate_counts()') 
+
+    def pre_process_hook(self):
+        pass
+
+    def process_counts(self):
+        """ 
+        The main thing process_counts() does is make sure that 
+        'processed_array' is loaded with information for generate_image()
+        to run with.
+
+        This base-class implementation just points one to the other.
+        """
+        self.processed_array = self.counts_array
+
+    def pre_image_hook(self):
+        pass
+
+    def generate_image(self):
+        """
+        Apart from actually generating an image, another important thing
+        that generate_image() does is sets output_image_file_name to the
+        appropriate file name that we generated.
+        """
+        raise NotImplementedError('Subclass must implement generate_image()') 
+
+    def ending_hook(self):
+        pass
+
+class EscapeAlgo(Algo):
+
+    @staticmethod
+    def options_list():
+        whole_list = Algo.options_list()
+
+        whole_list.extend(["escape-radius=", 
+            "max-escape-iterations=",
+            "burn",
+            "palette=",
+        ]) 
+
+        return whole_list
+
+    @staticmethod
+    def load_options_with_math_support(opts, math_support):
+        options = Algo.load_options_with_math_support(opts, math_support)
+
+        for opt,arg in opts:
+            if opt in ['--escape-radius']:
+                options['escape_radius'] = float(arg)
+            elif opt in ['--max-escape-iterations']:
+                options['max_escape_iterations'] = int(arg)
+            elif opt in ['--burn']:
+                options['burn_in'] = True
+            elif opt in ['--palette']:
+                options['palette'] = arg
+
+        return options
+
+    def __init__(self, dive_mesh, frame_number, output_folder_name, extra_params={}):
+        super().__init__(dive_mesh, frame_number, output_folder_name, extra_params)
+
+        # Load, with optional default values
+        self.escape_radius = extra_params.get('escape_radius', 10.0)
+        self.max_escape_iterations = int(extra_params.get('max_escape_iterations', 255))
+        self.burn_in = extra_params.get('burn_in', False)
+        self.palette = extra_params.get('palette', fp.FractalPalette())
+
+    def burn_text_to_drawing(self, burn_in_text, drawing):
+        burn_in_location = (10,10)
+        burn_in_margin = 5 
+        burn_in_font = ImageFont.truetype('fonts/cour.ttf', 12)
+        burn_in_size = burn_in_font.getsize_multiline(burn_in_text)
+        drawing.rectangle(((burn_in_location[0] - burn_in_margin, burn_in_location[1] - burn_in_margin), (burn_in_size[0] + burn_in_margin * 2, burn_in_size[1] + burn_in_margin * 2)), fill="black")
+        drawing.text((burn_in_location[0]-2, burn_in_location[1] - 2), burn_in_text, 'white', burn_in_font)
+
+class JuliaAlgo(EscapeAlgo):
+
+    @staticmethod
+    def options_list():
+        whole_list = EscapeAlgo.options_list()
+
+        whole_list.extend(["julia-center="]) 
+
+        return whole_list
+
+    @staticmethod
+    def load_options_with_math_support(opts, math_support):
+        # Considered loading this with default values, but didn't
+        # want to compete with the defaults in __init__()?
+        options = EscapeAlgo.load_options_with_math_support(opts, math_support)
+
+        for opt,arg in opts:
+            if opt in ['--julia-center']:
+                options['julia_center'] = math_support.createComplex(arg)
+
+        return options
+
+    def __init__(self, dive_mesh, frame_number, output_folder_name, extra_params={}):
+        super().__init__(dive_mesh, frame_number, output_folder_name, extra_params)
+
+        #self.julia_center = extra_params.get('julia_center', self.dive_mesh.mathSupport.createComplex(0,0)
+        self.julia_center = extra_params.get('julia_center', self.dive_mesh.mathSupport.createComplex(-.8,.145))
+
+        #print("settled on %s" % self.julia_center)
+
+