sf_graphs.cpp

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

//  (C) Copyright John Maddock 2008.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifdef _MSC_VER
#  pragma warning (disable : 4127) // conditional expression is constant
#  pragma warning (disable : 4180) // qualifier applied to function type has no meaning; ignored
#  pragma warning (disable : 4503) // decorated name length exceeded, name was truncated
#  pragma warning (disable : 4512) // assignment operator could not be generated
#  pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'function_ptr' was previously defined as a type
#endif

#include <boost/math/special_functions.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/function.hpp>
#include <boost/bind.hpp>

#include <list>
#include <map>
#include <string>
#include <boost/svg_plot/svg_2d_plot.hpp>
#include <boost/svg_plot/show_2d_settings.hpp>

class function_arity1_plotter
{
public:
   function_arity1_plotter() : m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0), m_has_legend(false) {}

   void add(boost::function<double(double)> f, double a, double b, const std::string& name)
   {
      if(name.size())
         m_has_legend = true;
      //
      // Now set our x-axis limits:
      //
      if(m_max_x == m_min_x)
      {
         m_max_x = b;
         m_min_x = a;
      }
      else
      {
         if(a < m_min_x) 
            m_min_x = a;
         if(b > m_max_x)
            m_max_x = b;
      }
      m_points.push_back(std::pair<std::string, std::map<double,double> >(name, std::map<double,double>()));
      std::map<double,double>& points = m_points.rbegin()->second;
      double interval = (b - a) / 200;
      for(double x = a; x <= b; x += interval)
      {
         //
         // Evaluate the function, set the Y axis limits 
         // if needed and then store the pair of points:
         //
         double y = f(x);
         if((m_min_y == m_max_y) && (m_min_y == 0))
            m_min_y = m_max_y = y;
         if(m_min_y > y)
            m_min_y = y;
         if(m_max_y < y)
            m_max_y = y;
         points[x] = y;
      }
   }

   void add(const std::map<double, double>& m, const std::string& name)
   {
      if(name.size())
         m_has_legend = true;
      m_points.push_back(std::pair<std::string, std::map<double,double> >(name, m));
      std::map<double, double>::const_iterator i = m.begin();

      while(i != m.end())
      {
         if((m_min_x == m_min_y) && (m_min_y == 0))
         {
            m_min_x = m_max_x = i->first;
         }
         if(i->first < m_min_x)
         {
            m_min_x = i->first;
         }
         if(i->first > m_max_x)
         {
            m_max_x = i->first;
         }

         if((m_min_y == m_max_y) && (m_min_y == 0))
         {
            m_min_y = m_max_y = i->second;
         }
         if(i->second < m_min_y)
         {
            m_min_y = i->second;
         }
         if(i->second > m_max_y)
         {
            m_max_y = i->second;
         }

         ++i;
      }
   }

   void plot(const std::string& title, const std::string& file,
      const std::string& x_lable = std::string(),
      const std::string& y_lable = std::string())
   {
      using namespace boost::svg;

      static const svg_color colors[5] = 
      {
         darkblue,
         darkred,
         darkgreen,
         darkorange,
         chartreuse
      };

      svg_2d_plot plot;
      plot.image_size(600, 400);
      plot.copyright_holder("John Maddock").copyright_date("2008").boost_license_on(true);
      plot.coord_precision(4); // Could be 3 for smaller plots?
      plot.title(title).title_font_size(20).title_on(true);
      plot.legend_on(m_has_legend);

      double x_delta = (m_max_x - m_min_x) / 50;
      double y_delta = (m_max_y - m_min_y) / 50;
      plot.x_range(m_min_x, m_max_x + x_delta)
         .y_range(m_min_y, m_max_y + y_delta);
      plot.x_label_on(true).x_label(x_lable);
      plot.y_label_on(true).y_label(y_lable);
      plot.y_major_grid_on(false).x_major_grid_on(false);
      plot.x_num_minor_ticks(3);
      plot.y_num_minor_ticks(3);
      //
      // Work out axis tick intervals:
      //
      double l = std::floor(std::log10((m_max_x - m_min_x) / 10) + 0.5);
      double interval = std::pow(10.0, (int)l);
      if(((m_max_x - m_min_x) / interval) > 10)
         interval *= 5;
      plot.x_major_interval(interval);
      l = std::floor(std::log10((m_max_y - m_min_y) / 10) + 0.5);
      interval = std::pow(10.0, (int)l);
      if(((m_max_y - m_min_y) / interval) > 10)
         interval *= 5;
      plot.y_major_interval(interval);
      plot.plot_window_on(true);
      plot.plot_border_color(lightslategray)
          .background_border_color(lightslategray)
          .legend_border_color(lightslategray)
          .legend_background_color(white);

      int color_index = 0;

      for(std::list<std::pair<std::string, std::map<double,double> > >::const_iterator i = m_points.begin();
         i != m_points.end(); ++i)
      {
         plot.plot(i->second, i->first)
            .line_on(true)
            .line_color(colors[color_index])
            .line_width(1.)
            .shape(none);
         if(i->first.size())
            ++color_index;
         color_index = color_index % (sizeof(colors)/sizeof(colors[0]));
      }
      plot.write(file);
   }

   void clear()
   {
      m_points.clear();
      m_min_x = m_min_y = m_max_x = m_max_y = 0;
      m_has_legend = false;
   }

private:
   std::list<std::pair<std::string, std::map<double, double> > > m_points;
   double m_min_x, m_max_x, m_min_y, m_max_y;
   bool m_has_legend;
};

template <class F>
struct location_finder
{
   location_finder(F _f, double t, double x0) : f(_f), target(t), x_off(x0){}

   double operator()(double x)
   {
      try
      {
         return f(x + x_off) - target;
      }
      catch(const std::overflow_error&)
      {
         return boost::math::tools::max_value<double>();
      }
      catch(const std::domain_error&)
      {
         if(x + x_off == x_off)
            return f(x_off + boost::math::tools::epsilon<double>() * x_off);
         throw;
      }
   }

private:
   F f;
   double target;
   double x_off;
};

template <class F>
double find_end_point(F f, double x0, double target, bool rising, double x_off = 0)
{
   boost::math::tools::eps_tolerance<double> tol(50);
   boost::uintmax_t max_iter = 1000;
   return boost::math::tools::bracket_and_solve_root(
      location_finder<F>(f, target, x_off), 
      x0, 
      double(1.5), 
      rising, 
      tol, 
      max_iter).first;
}

double sqrt1pm1(double x)
{
   return boost::math::sqrt1pm1(x);
}

double lbeta(double a, double b)
{
   return std::log(boost::math::beta(a, b));
}

int main()
{
   function_arity1_plotter plot;
   double (*f)(double);
   double (*f2)(double, double);
   double (*f2u)(unsigned, double);
   double (*f2i)(int, double);
   double (*f3)(double, double, double);
   double (*f4)(double, double, double, double);
   
   f = boost::math::zeta;
   plot.add(f, 1 + find_end_point(f, 0.1, 40.0, false, 1.0), 10, "");
   plot.add(f, -20, 1 + find_end_point(f, -0.1, -40.0, false, 1.0), "");
   plot.plot("Zeta Function Over [-20,10]", "zeta1.svg", "z", "zeta(z)");

   plot.clear();
   plot.add(f, -14, 0, "");
   plot.plot("Zeta Function Over [-14,0]", "zeta2.svg", "z", "zeta(z)");

   f = boost::math::tgamma;
   double max_val = f(6);
   plot.clear();
   plot.add(f, find_end_point(f, 0.1, max_val, false), 6, "");
   plot.add(f, -1 + find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, -max_val, false), "");
   plot.add(f, -2 + find_end_point(f, 0.1, max_val, false, -2), -1 + find_end_point(f, -0.1, max_val, true, -1), "");
   plot.add(f, -3 + find_end_point(f, 0.1, -max_val, true, -3), -2 + find_end_point(f, -0.1, -max_val, false, -2), "");
   plot.add(f, -4 + find_end_point(f, 0.1, max_val, false, -4), -3 + find_end_point(f, -0.1, max_val, true, -3), "");
   plot.plot("tgamma", "tgamma.svg", "z", "tgamma(z)");

   f = boost::math::lgamma;
   max_val = f(10);
   plot.clear();
   plot.add(f, find_end_point(f, 0.1, max_val, false), 10, "");
   plot.add(f, -1 + find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), "");
   plot.add(f, -2 + find_end_point(f, 0.1, max_val, false, -2), -1 + find_end_point(f, -0.1, max_val, true, -1), "");
   plot.add(f, -3 + find_end_point(f, 0.1, max_val, false, -3), -2 + find_end_point(f, -0.1, max_val, true, -2), "");
   plot.add(f, -4 + find_end_point(f, 0.1, max_val, false, -4), -3 + find_end_point(f, -0.1, max_val, true, -3), "");
   plot.add(f, -5 + find_end_point(f, 0.1, max_val, false, -5), -4 + find_end_point(f, -0.1, max_val, true, -4), "");
   plot.plot("lgamma", "lgamma.svg", "z", "lgamma(z)");

   f = boost::math::digamma;
   max_val = 10;
   plot.clear();
   plot.add(f, find_end_point(f, 0.1, -max_val, true), 10, "");
   plot.add(f, -1 + find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, max_val, true), "");
   plot.add(f, -2 + find_end_point(f, 0.1, -max_val, true, -2), -1 + find_end_point(f, -0.1, max_val, true, -1), "");
   plot.add(f, -3 + find_end_point(f, 0.1, -max_val, true, -3), -2 + find_end_point(f, -0.1, max_val, true, -2), "");