# partition.rb: segment partitioning functions
# Copyright (c) 2008 by Vincent Fourmond

  
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
  
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details (in the COPYING file).

require 'CTioga/utils'


module CTioga

  Version::register_svn_info('$Revision: 792 $', '$Date: 2008-05-23 00:19:31 +0200 (Fri, 23 May 2008) $')

  # A module for small convenience functions.
  module Utils

    # The following code is stolen from SciYAG/lib/utils.rb,
    # and is used to partition segments in natural-looking
    # subdivisions.

    
    # Our natural way to split decades
    NaturalDistances = Dobjects::Dvector[1, 2, 2.5, 5, 10]

    # Attempts to partition the given segment in at most _nb_
    # segments of equal size. The segments don't necessarily start on
    # the edge of the original segment
    def self.partition_segment(min, max, nb)
      if min > max             
        return partition_segment(max, min, nb)
      elsif min.nan? or max.nan?
        return [0] * nb
      elsif min == max
        return self.partition_segment(min * 0.7, min * 1.3, nb)
      end
      distance = max - min
      min_distance = distance/(nb + 1.0) # Why + 1.0 ? To account
      # for the space that could be left on the side.
      
      # The order of magnitude of the distance:
      order = min_distance.log10.floor
      
      # A distance which is within [1, 10 [ (but the latter is never reached.
      normalized_distance = min_distance * 10**(-order)
      final_distance = NaturalDistances.min_gt(normalized_distance) * 
        10**(order)
      #       puts "Distance: #{distance} in #{nb} : #{normalized_distance} #{final_distance}"
      
      # We're getting closer now: we found the natural distance between
      # ticks.
      
      start = (min/final_distance).ceil * final_distance
      retval = []
      val = start
      while val <= max
        retval << val
        # I use this to avoid potential cumulative addition
        # rounding errors
        val = start + final_distance * retval.size
      end
      return retval

    end

    # Our natural way to split decades - except that all successive
    # element now divide each other.
    NaturalDistancesNonLinear = Dobjects::Dvector[1, 2.5, 5, 10]

    # A class that handles a "natural distance". Create it by giving
    # the initial distance. You can then get its value, lower it,
    # increase it, find the first/last element of an interval that is
    # a multiple if it
    class NaturalDistance

      # Creates the biggest 'natural distance' that is smaller
      # than _distance_
      def initialize(distance)
        @order = distance.log10.floor
        normalized = distance * 10**(-@order)
        @index = 0
        for dist in NaturalDistances
          if dist > normalized
            break
          else
            @index += 1
          end
        end
      end

      # Returns the actual value of the distance
      def value
        return NaturalDistances[@index] * 10**(@order)
      end

      # Goes to the natural distance immediately under the current one
      def decrease
        if @index > 0
          @index -= 1
        else
          @index = NaturalDistances.size - 1
          @order -= 1
        end
      end

      # Goes to the natural distance immediately under the current one
      def increase
        if @index < NaturalDistances.size - 1
          @index += 1
        else
          @index = 0
          @order += 1
        end
      end

      # Find the minimum element inside the given interval
      # that is a multiple of this distance
      def find_minimum(x1, x2)
        x1,x2 = x2, x1 if x1 > x2

        dist = value
        v = (x1/dist).ceil * dist
        if v <= x2
          return v
        else
          return false
        end
      end

      # Returns the value of the next decade (ascending)
      def next_decade(x)
        decade = 10.**(@order + 1)
        if @index == 0          # We stop at 0.5 if the increase is 0.5
          decade *= 0.5
        end
        nb = Float((x/decade).ceil)
        if nb == x/decade
          return (nb + 1)*decade
        else
          return nb * decade
        end
      end

      # Returns the number of elements to #next_decade
      def nb_to_next_decade(x)
        dist = value
        dec = next_decade(x)
        return dec/dist - (x/dist).ceil + 1
      end

      # Returns a list of all the values corresponding to
      # the distance
      def to_next_decade(x)
        next_decade = next_decade(x)
        dist = value
        start = (x/dist).ceil * dist
        retval = []
        while start + retval.size * dist <= next_decade
          retval << start + retval.size * dist
        end
        return retval
      end

    end


    # Attempts to partition the segment image of _min_, _max_
    # by the Proc object _to_ into at most _nb_ elements. The
    # reverse of the transformation, _from_, has to be provided.
    def self.partition_nonlinear(to, from, x1, x2, nb)
      x1, x2 = x2, x1 if x1 > x2
      if x1.nan? or x2.nan?
        return [0] * nb
      elsif x1 == x2            # Nothing to do
        return self.partition_segment(x1 * 0.7, x2 * 1.3, nb)
      end

      xdist = x2 - x1
      xdist_min = xdist/(nb + 1.0) # Why + 1.0 ? To account
      # for the space that could be left on the side.

      y1 = to.call(x1)
      y2 = to.call(x2)
      
      # Make sure y1 < y2
      y1, y2 = y2, y1 if y1 > y2

      # We first need to check if the linear partitioning of
      # the target segment could be enough:
      
      candidate = self.partition_segment(y1, y2, nb)
      candidate_real = candidate.map(&from)

      # We inspect the segment: if one of the length deviates from the
      # average expected by more than 25%, we drop it

      length = []
      p candidate_real, xdist_min
      0.upto(candidate.size - 2) do |i|
        length << (candidate_real[i+1] - candidate_real[i]).abs/(xdist_min)
      end
      p length
      # If everything stays within 25% off, we keep that
      if length.min > 0.75 and length.max < 1.7
        return candidate
      end
      

      # We start with a geometric measure of the distance, that
      # will most likely scale better:
      ydist = y1 * (y2/y1).abs ** (1/(nb + 1.0))

      cur_dist = NaturalDistance.new(ydist)
      
      retval = []
      
      cur_y = y1
      # This flag is necessary to avoid infinite loops
      last_was_decrease = false
      
      distance_unchanged = 0
      last_real_distance = false
      while cur_y < y2
        candidates = cur_dist.to_next_decade(cur_y)
        # We now evaluate the distance in real
        real_distance = (from.call(cur_y) - from.call(candidates.last)).abs/
          candidates.size
        if last_real_distance && (real_distance == last_real_distance)
          distance_unchanged += 1
        else
          distance_unchanged = 0
        end
#         p [:cur_y=, cur_y, :y2=, y2, :real_distance, real_distance, 
#            :distance=, cur_dist, :xdist_min, xdist_min, 
#            :candidates=, *candidates]

        if (real_distance > 1.25 * xdist_min) && 
            (distance_unchanged < 3)
          cur_dist.decrease
          last_was_decrease = true
        elsif real_distance < 0.75 * xdist_min && 
            !last_was_decrease && (distance_unchanged < 3) && 
            candidates.last <= 10 * y2
          cur_dist.increase
          last_was_decrease = false
        else
          retval += candidates
          cur_y = candidates.last
          last_was_decrease = false
        end
        last_real_distance = real_distance
      end

      # We need to select them so
      return retval.select do |y|
        y >= y1 and y <= y2
      end

    end


  end

end