10
\$\begingroup\$

Follow up on Plotting polynomial roots

I rewrote the program to C++ and am now using GNU GSL scientific library and gnuplot and also used OpenMP for multithreading. The result is that it takes about 7 seconds compared to 14 hours when I used python with numpy and matplotlib.

I added the functionality to plot the image with different colours based on the density of the points with a radius of 0.1. The algorithm for this is O(n^2) which is very slow. I basically brute force to check the distance for each point to every other point in the system and check if it is <= 0.1 and if it is I add to the neighbour count for that point. In the end I map the neighbour count to a scalar between 0 and 1.5 which represents the colour.

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_poly.h>
#include <cmath>

#define maxd(a,b) (((a) > (b)) ? (a) : (b))
#define mind(a,b) (((a) < (b)) ? (a) : (b))

unsigned long long cur_root;
unsigned long long capacity;

void
polynomial_recursive(double **roots, double min_x, double max_x, int min_degree, int max_degree, int cur_degree, double *coeff)
{
    if (cur_degree <= max_degree)
      {    
        for (double c = -4; c <= 4; c++)
          {
            double *new_coeff;

            new_coeff = (double *) malloc (sizeof (double) * (cur_degree + 1));
            if (new_coeff == NULL)
              {
                perror ("malloc failed");
              }
            memcpy (new_coeff, coeff, cur_degree * sizeof (double));
            new_coeff[cur_degree] = c;
            if (cur_degree >= min_degree && new_coeff[cur_degree] != 0)
              {
                double z[cur_degree*2];

                gsl_poly_complex_workspace *w
                    = gsl_poly_complex_workspace_alloc (cur_degree + 1);

                gsl_poly_complex_solve (new_coeff, cur_degree + 1, w, z);

                gsl_poly_complex_workspace_free (w);

                for (int i = 0; i < cur_degree; i++)
                  {
                    if (cur_root + 2 >= capacity)
                      {
                        capacity *= 2;
                        double *tmp = (double *) realloc (*roots, capacity * sizeof (double));
                        if (tmp == NULL)
                          {
                            perror ("realloc failed");
                          }
                        else
                          {
                            *roots = tmp;
                          }
                      }
                    (*roots) [cur_root++] = z[2*i];
                    (*roots) [cur_root++] = z[2*i+1];
                  }
              }
            #pragma omp task shared(new_coeff)
            polynomial_recursive (roots, min_x, max_x, min_degree, max_degree, cur_degree + 1, new_coeff);
            #pragma omp taskwait
            free (new_coeff);
          }
      }
}

unsigned long long getDensityForPoint(double x, double y, double *points, int size, double radius)
{
    unsigned long long count = 0;
    for (int i = 0; i < size; i += 2)
      {
        double x2 = points[i];
        double y2 = points[i+1];
        if (x != x2 || y != y2)
          {
            double dist = fabs(x2 - x) + fabs(y2 - y);;
            if (dist <= radius)
              {
                count++;
              }
          }
      }
    return count;
}

double
calculateDensityScalarFromDensity(unsigned long min, unsigned long max, unsigned long density)
{
    return 0.0 + ((1.5 - 0.0) / (double)(max - min)) * (density - min);
} 

int
main(void)
{
    double *roots;
    unsigned long long *densities, min, max;

    capacity = 100;
    cur_root = 0;
    roots = (double *) malloc (capacity * sizeof (double));
    if (roots == NULL)
      {
        perror ("malloc failed");
      }

    #pragma omp parallel
    {
        #pragma omp single 
        {
            polynomial_recursive (&roots, -4, 4, 2, 5, 0, NULL);
        }
    }

    densities = (unsigned long long *) malloc (cur_root / 2 * sizeof (unsigned long long));
    max = 0;
    min = ~0;

    #pragma omp parallel for
    for (int i = 0; i < cur_root; i += 2)
      {
        densities[i / 2] = getDensityForPoint(roots[i], roots[i+1], roots, cur_root, 0.1);
        max = maxd(densities[i / 2], max);
        min = mind(densities[i / 2], min);
      }

    printf("max %llu, min %llu\n", max, min);
    FILE *gnuplotPipe = popen ("gnuplot -persistent", "w");
    FILE *temp = fopen("data.temp", "w");

    #pragma omp parallel for
    for (int i = 0; i < cur_root; i += 2)
      {
        fprintf (temp, "%lf %lf %lf\n", roots[i], roots[i+1], calculateDensityScalarFromDensity(min, max, densities[i / 2])); 
      }

    free(densities);
    free(roots);

    fprintf (gnuplotPipe, "set zrange [0.0:1.5]\nset cbrange [0.0:1.5]\nset term pngcairo\nset terminal png size 1920,1080\nset output \"plot.png\"\nset palette defined (0 \"black\", 0.5 \"red\", 1 \"yellow\", 1.5 \"white\")\nplot 'data.temp' with points pointtype 5 pointsize 0.001 palette\n");
}

It is this part that is taking so long:

for (int i = 0; i < cur_root; i += 2)
  {
    densities[i / 2] = getDensityForPoint(roots[i], roots[i+1], roots, cur_root, 0.1);
    max = maxd(densities[i / 2], max);
    min = mind(densities[i / 2], min);
  }

I am using the manhattan distance because that is the fastest I could think of. The problem is that I'm doing this for very complex polynomials. I set min_degree to 20 and max_degree to 22 and instead of using between -4 and 4 I use only either -PI or +PI so nothing in between.

Going up to polynomials of 22 order and only using two different coefficients I calculated that there would be about 300 million roots and so with a O(n^2) time complexity it would take around 75 days to complete which is insane. I've heard about k-d trees but haven't successfully implemented them to gain any major speed advantage. Another idea I had is to downscale the data.

How do I speed up the density calculation?

Here is another image that I generated with max_degree 5 and coefficients between -4 and 4 and the density function enabled (black is low density, white is high density):

enter image description here

Another awesome image: enter image description here

\$\endgroup\$
  • 2
    \$\begingroup\$ A better algorithm would be preferable, but there are certainly ways to micro-optimize this code to squeeze more performance out of it. First things first, is your compiler actually taking the #pragma omp parallel for hint to parallelize this loop? (BTW, parallelizing a for loop that calls fprintf is somewhere between silly and a bad idea…). \$\endgroup\$ – Cody Gray Jun 2 '17 at 15:27
  • 2
    \$\begingroup\$ As you are writing a bitmap eventually, it may/will be a lot faster to blit the results onto the bitmap right away. With some alpha blending you will get approximately the same result. Or take a grid of integers (i.e. buckets) to count values and transform that into an image. The easy fix at this point is to stop collecting the actual roots, as their number explodes. \$\endgroup\$ – mvds Jun 2 '17 at 15:31
  • 1
    \$\begingroup\$ Also, for the last coefficient you only need to loop over positive values, due to symmetries in your problem. It's an easy performance gain of 50%. \$\endgroup\$ – mvds Jun 2 '17 at 15:33
  • 1
    \$\begingroup\$ Also, when either collecting roots or collecting a bitmap, you can limit yourself to the right upper quadrant, as the resulting plot can be proven to be symmetric about the x and y axis. \$\endgroup\$ – mvds Jun 2 '17 at 15:34
  • \$\begingroup\$ I think this should be tag C and not C++. \$\endgroup\$ – coincoin Jun 2 '17 at 15:52
7
\$\begingroup\$

Keeping all else untouched, this gives a 1000x performance boost:

    int w=1920,h=1080;

    uint64_t *buckets = (uint64_t*)calloc(w*h,sizeof(buckets));
    double min_x,max_x,min_y,max_y;
    min_x = max_x = roots[0];
    min_y = max_y = roots[1];

    for ( int i = 0; i < cur_root; i+=2 )
    {
            min_x = roots[i]<min_x?roots[i]:min_x;
            max_x = roots[i]>max_x?roots[i]:max_x;
            min_y = roots[i+1]<min_y?roots[i+1]:min_y;
            max_y = roots[i+1]>max_y?roots[i+1]:max_y;
    }

    printf("range: (%f,%f)-(%f,%f)\n",min_x,min_y,max_x,max_y);

    double pad_x = (max_x-min_x)/100;
    double pad_y = (max_y-min_y)/100;
    min_x -= pad_x;
    max_x += pad_x;
    min_y -= pad_y;
    max_y += pad_y;

    for ( int i = 0; i < cur_root; i+=2 )
    {
            int x = (roots[i]-min_x)/(max_x-min_x)*w;
            int y = (roots[i+1]-min_y)/(max_y-min_y)*h;
            buckets[x+w*y]++;
    }

    uint64_t v_min,v_max;
    v_max = v_min = buckets[0];
    for ( int i = 0; i < w*h; i++ )
    {
            v_min = buckets[i]<v_min?buckets[i]:v_min;
            v_max = buckets[i]>v_max?buckets[i]:v_max;
    }

    printf("max %llu, min %llu\n", v_max, v_min);
    FILE *gnuplotPipe = popen ("gnuplot -persistent", "w");
    FILE *temp = fopen("data.temp", "w");

    for ( int i = 0; i < w*h; i++ )
    {
            uint64_t val = buckets[i];
            if ( !val ) continue;
            int x = i%w;
            int y = i/w;
            double fx = min_x + (max_x-min_x) * x / w;
            double fy = min_y + (max_y-min_y) * y / h;
            fprintf (temp, "%lf %lf %lf\n", fx, fy, calculateDensityScalarFromDensity(v_min, v_max, val));
    }

The approach to take the bins and turn them into a list of points is a little backward, but it leaves the gnuplot command untouched. You chould look into a way to plot the data directly as a bitmap (which might be done using gnuplot, but I don't know gnuplot that well).

Also, as I tried to explain in the comments, your coefficient space is redundant; the roots for [c0,c1,c2,c3] are the same as the roots for [a*c0,a*c1,a*c2,a*c3]. Also, if the last coefficient is 0, you're basically solving a polynomial of a lower order. So you could loop the last coefficient starting from 1:

void polynomial_recursive(double **roots, double min_x, double max_x, int min_degree, int max_degree, int cur_degree, double *coeff)
{
    for (double c = cur_degree==max_degree?1:min_x; c <= max_x; c++)
    {
            double *new_coeff;

            new_coeff = (double *) malloc (sizeof (double) * (cur_degree + 1));
            if (new_coeff == NULL)
            {
                    perror ("malloc failed");
            }
            memcpy (new_coeff, coeff, cur_degree * sizeof (double));
            new_coeff[cur_degree] = c;

            if ( cur_degree == max_degree )
            {
                    double z[cur_degree*2];

                    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (cur_degree + 1);
                    gsl_poly_complex_solve (new_coeff, cur_degree + 1, w, z);
                    gsl_poly_complex_workspace_free (w);

                    for (int i = 0; i < cur_degree; i++)
                    {
                            if (cur_root + 2 >= capacity)
                            {
                                    capacity *= 2;
                                    double *tmp = (double *) realloc (*roots, capacity * sizeof (double));
                                    if (tmp == NULL)
                                    {
                                            perror ("realloc failed");
                                    }
                                    else
                                    {
                                            *roots = tmp;
                                    }
                            }
                            (*roots) [cur_root++] = z[2*i];
                            (*roots) [cur_root++] = z[2*i+1];
                    }
            } else {
                    #pragma omp task shared(new_coeff)
                    polynomial_recursive (roots, min_x, max_x, min_degree, max_degree, cur_degree + 1, new_coeff);
                    #pragma omp taskwait
            }
            free (new_coeff);
    }
}

This shaves another 50% of the runtime.

I think you can also rearrange the code to do less initialisations and mallocs, but that will not save as much.

Turns out gnuplot does support this, the command looks like:

plot 'data.temp' binary array=%dx%d format=\"%%lu\" with image

This code will be a lot faster in writing the pixel data and give some stunning images: (losing the actual axis values as we are plotting pixels from 0..3000)

    FILE *gnuplotPipe = popen ("gnuplot -persistent", "w");
    FILE *temp = fopen("data.temp", "w");
    fwrite(buckets,w*h,sizeof(*buckets),temp);
    fclose(temp);

    free(roots);

    printf("Plotting...\n");

    fprintf (gnuplotPipe,
            "set cbrange [0.0:30]\n"
            "set xrange [0:%d]\n"
            "set yrange [0:%d]\n"
            "set term pngcairo\n"
            "set terminal png size %d,%d\n"
            "set output \"plot.png\"\n"
            "set palette defined (0 \"white\", 1 \"black\", 100 \"red\", 200 \"yellow\", 300 \"green\", 400 \"blue\")\n"
            "plot 'data.temp' binary array=%dx%d format=\"%%lu\" with image\n",w,h,w+218,h+60,w,h);

roots degree 6 9mp

\$\endgroup\$
  • \$\begingroup\$ Thanks for the tips, for some reason when I ran the code with the bucket implementation with max_degree being set to 5 I got a black image. Here is a pastebin: pastebin.com/HcXfKqvT and this is what it looks like: imgur.com/a/eLAnr \$\endgroup\$ – Linus Jun 2 '17 at 18:07
  • \$\begingroup\$ You have a few points with extremely high value (e.g. 0,0). You should either hardcode some low v_max or use a logarithmic scale. \$\endgroup\$ – mvds Jun 2 '17 at 18:18
  • \$\begingroup\$ That works, but I lose a lot of detail: imgur.com/a/Nx61B \$\endgroup\$ – Linus Jun 2 '17 at 18:23
  • \$\begingroup\$ You should look into a way to plot the bitmap directly \$\endgroup\$ – mvds Jun 2 '17 at 18:27
  • \$\begingroup\$ Not sure if gnuplot supports that but I'll try. What do you mean by alpha blending? \$\endgroup\$ – Linus Jun 2 '17 at 18:33
4
\$\begingroup\$

Use what C++ provides

If you use the C++ versions of the headers:

#include <cstdlib>
#include <cstdio> // or better, <iostream>
#include <cstring>

You'll get the names safely in the std namespace. I see you include <cmath> but assume its names are exported to the global namespace - that's an error. You can't depend on that happening (it's just an artefact of your implementation).

Also, prefer new/delete and new[]/delete[] over malloc()/free(). You might be able to avoid memory management completely by using the standard containers; that's a Good Thing.

std::max() and std::min() are safer choices compared to your maxd() and mind() macros.

Avoid I/O in parallel code

These prints may get interleaved in unpredictable ways:

#pragma omp parallel for
for (int i = 0; i < cur_root; i += 2)
  {
    fprintf (temp, "%lf %lf %lf\n",
             roots[i], roots[i+1],
             calculateDensityScalarFromDensity(min, max, densities[i / 2])); 
  }

You might do something like

#pragma omp parallel for
for (int i = 0; i < cur_root; i += 2)
  {
    auto density = calculateDensityScalarFromDensity(min, max, densities[i / 2])
#pragma omp critical
    fprintf (temp, "%lf %lf %lf\n",
             roots[i], roots[i+1],
             density); 
  }

Bin your results

I don't mean to throw everything away! The output is a count of results at each pixel location (possibly with some small blurring? I'm not quite sure). Having obtained values, the only use for them is to increment the count of the corresponding bin - it's like obtaining a histogram. That reduces your storage requirements, and your output is almost formatted for you.

Use symmetry

As noted in comments, we can get away with recording only the positive roots, and produce the rest of the output by reflection.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.