*** thuban.orig/Thuban/Model/classgen.py	2004-11-24 18:33:49.000000000 +0100
--- thuban/Thuban/Model/classgen.py	2004-11-24 18:38:50.000000000 +0100
***************
*** 17,24 ****
  # $Id: classgen.py,v 1.19 2004/10/07 14:28:51 jan Exp $
  
  import operator
  
! from color import Color, Transparent
  from range import Range
  from classification import Classification, ClassGroupSingleton, \
      ClassGroupRange, ClassGroupProperties
--- 17,26 ----
  # $Id: classgen.py,v 1.19 2004/10/07 14:28:51 jan Exp $
  
  import operator
+ from math import sqrt
+ import random
  
! from color import Color, Transparent, interpolate_rgb, interpolate_hsv
  from range import Range
  from classification import Classification, ClassGroupSingleton, \
      ClassGroupRange, ClassGroupProperties
***************
*** 155,160 ****
--- 157,185 ----
  
      return (adjusted, clazz)
  
+ def generate_ranges_from_list( _list, ramp):
+     """Return a classification based on the breaks indicated in the list.
+ 
+     The list *must* include min and max values. The ranges generated will
+     be of the form [a,b[ (lower bound included, upper bound excluded) except
+     for the greater one, which will include the upper bound, i.e. the max.
+ 
+     n-1 ranges will be generated for a list of n elements.
+     """
+     clazz = Classification()
+     numGroups = len( _list )
+     if numGroups == 0: return clazz
+     n = numGroups - 1
+     end = "["
+     for i in range(n):
+         prop = ramp.GetProperties( float(i)/(n-1) )
+         if i == n - 1:
+             end = "]"
+         clazz.AppendGroup( ClassGroupRange(
+             Range( ("[", _list[i], _list[i+1],end) ),
+             prop))
+     return clazz
+ 
  
  def calculate_quantiles(_list, percents, _range):
      """Calculate quantiles for the given _list of percents from the
***************
*** 313,324 ****
                  [(q, (q - minIndex+1) / float(numValues)) \
                   for q in quantiles])
  
  class CustomRamp:
  
!     def __init__(self, prop1, prop2):
          """Create a ramp between prop1 and prop2."""
          self.prop1 = prop1
          self.prop2 = prop2
  
      def GetRamp(self):
          """Return this ramp."""
--- 338,775 ----
                  [(q, (q - minIndex+1) / float(numValues)) \
                   for q in quantiles])
  
+ 
+ # TODO: implement Jenks & Caspall completely.
+ 
+ if not  __builtins__.has_key("sum"): # Python < 2.3
+     sum = lambda x: reduce( lambda y,z: y+z, x )
+ 
+ def lineq( x1, y1, x2, y2 ):
+     """Return (a,b) such that y_i = a * x_i + b; (None,b) for x_i = b """
+     if x1 == x2:
+         return (None, y1)
+     else:
+         a = (y2-y1) / (x2-x1)
+         return (a, y1 - a*x1)
+ 
+ 
+ def classify_discont( data, maxclasses = None, resultlog = None ):
+      """Generate "natural" breaks in the sense that they approximate very well
+      the original distribution.
+  
+      Based on original Fortran code 'discont.f' by Jean-Pierre Grimmeau,
+      courtesy of Moritz Lennert.
+      
+      The distribution is approximated recursively by straight lines. 
+  
+      The procedure has no parameters: the tolerance for the stopping criterion
+      is embedded. It's the chi^2 value for alpha=95% with 1 d.o.f, used in
+      comparing a segment with its possible best subdivision.
+  
+      If 'maxclasses' is specified, the algorithm will stop after finding
+      the specified number of classes.
+  
+      'resultlog', if specified, must be a reference to a list
+      that will gather information about the processing, namely tuples of
+  
+      (number_of_classes,
+      class_i,
+      number_of_points_in_class_i,
+      min_value_for_class_i,
+      max_value_for_class_i,
+      mean_value_for_class_i,
+      standard_deviation_for_class_i).
+  
+      This data will be appended to whatever is the value of resultlog at
+      call time.
+      """
+ 
+      from math import sqrt
+ 
+      x = [ float(d) for d in data if d is not None]
+      x.sort()
+      xmin = x[0]
+      xmax = x[-1]
+      xrng = xmax - xmin
+      N = len(x)
+  
+      xinc = min([ x[i+1] - x[i] for i in range(N-1) if x[i+1] != x[i] ])
+      x = [0.0]+[ (x_i - xmin) / xrng for x_i in x ]
+      xn = [0.0]+[ float(i+1) / N for i in range(N) ]
+  
+      xlim = xinc / xrng
+      xinc /= 2.0
+  
+      nclass = [N] # indices for breaks
+      brk = [] # break values
+      no = [] # class count
+      co = [] # indicator for classification of border cases
+      nstart = 0
+      nmax = 0
+      while(1):
+          nend = 0
+          dmax = 0
+          nmax = 0
+          for j in range(len(nclass)):
+              nstart = nend
+              nend = nclass[j]
+              co.append(1e38)
+              (a,b) = lineq( x[nstart], xn[nstart], x[nend], xn[nend] )
+              if a is not None:
+                   den = sqrt( 1.0 + a*a )
+              for k in range(nstart+1, nend+1):
+                  if a is None:
+                      d = abs(x[k]-b)
+                  else:
+                      d = abs( xn[k] - (a*x[k] + b) ) / den
+                      # den is not needed to calculate the max here: it's constant
+                  if d > dmax and x[nend] - x[k] > xlim:
+                      dmax = d
+                      nmax = k
+  
+              if x[nend] != x[nstart]:
+                   co[j] = ( xn[nend] - xn[nstart] ) / (x[nend] - x[nstart] )
+          # Find break values (brk) and class count (no)
+          i = len(nclass) - 1
+          no = [0] * (i+1)
+          brk = [0] * i
+          for j in range(i):
+              no[j] = nclass[j]
+              brk[j] = x[nclass[j]] * xrng + xmin
+              if co[j] <= co[j+1]:
+                  brk[j] = brk[j] - xinc
+                  no[j] = no[j] - 1
+              else:
+                  brk[j] = brk[j] + xinc
+          no[ i ] = N
+          no.insert(0,0)
+          no = [ no[k+1] - no[k] for k in range(len(no)-1) ]
+  
+          if nmax==0:
+              break
+          nclass.append(nmax)
+          nclass.sort()
+  
+          if resultlog is not None and len(brk) > 0:
+              for i in [-1] + range(len(brk)):
+                  if i == -1:
+                      classdata = [ d * xrng + xmin for d in x
+                                    if d * xrng + xmin < brk[0] ]
+                      startpoint = xmin
+                      endpoint = brk[0]
+                  elif i < len(brk) - 1:
+                      classdata = [ d * xrng + xmin for d in x
+                                    if brk[i] < d * xrng + xmin < brk[i+1] ]
+                      startpoint = brk[i]
+                      endpoint = brk[i+1]
+                  else:
+                      classdata = [ d * xrng + xmin for d in x
+                                    if brk[i] < d * xrng + xmin ]
+                      startpoint = brk[i]
+                      endpoint = xmax
+                  ndata = len( classdata )
+                  if ndata > 1:
+                      meanvalue = sum( classdata ) / ndata
+                      stdev = sum( [ ( d - meanvalue ) ** 2
+                                     for d in classdata ] ) / (ndata - 1)
+                  else:
+                      stdev = None
+                      if ndata == 1:
+                          meanvalue = classdata[0]
+                      else:
+                          meanvalue = None
+                  resultlog.append(
+                      ( len(nclass) - 1, i + 2, ndata,
+                        startpoint, endpoint,
+                        meanvalue, stdev ))
+  
+          if maxclasses is not None and len(nclass) == maxclasses+1:
+              break
+  
+          j = nclass.index(nmax)
+          xnhere = xn[nclass[j]]
+          xnpost = xn[nclass[j+1]]
+          xpre  = x[nclass[j-1]]
+          xhere = x[nclass[j]]
+          xpost = x[nclass[j+1]]
+          xnpre = xn[nclass[j-1]]
+          no1 = N * ( xnhere - xnpre )
+          no2 = N * ( xnpost - xnhere )
+          if xpost != xpre:
+              f = N * ( xnpost - xnpre ) / ( xpost - xpre )
+          else:
+              f = 1e38
+          xt1 = f*( xhere - xpre )
+          xt2 = f*( xpost - xhere )
+          if abs( xt1*xt2 ) < 1e-15:
+               if abs(xt2)<1e-15:
+                    xt2 = xinc/2 * f/xrng
+                    xt1 = xt1 - xt2
+               else:
+                    xt1 = xinc/2 * f/xrng
+                    xt2 = xt2 - xt1
+          if xt1 + xt2 > 1e-15:
+               chi2 = ( (no1-no2) - (xt1-xt2) ) ** 2 / (xt1+xt2)
+               if chi2 < 3.84146:
+                    break
+              
+      out = [xmin] + brk + [xmax]
+      out.sort()
+      return out
+ 
+ 
+ def diffabs(x, y):
+     """A measure of distance between two values: absolute value of difference.
+ 
+     May be used bu classify_kmeans, classify_annealing.
+     """
+     return abs(x-y)
+ 
+ def classify_kmeans( data, k = 5, distance = diffabs,
+                      init = None, output = "breaks" ):
+     """k-means clustering with distance function 'distance', initial positions
+     'init' for 'k' centroids. If 'output' is "breaks", the breakpoints are
+     returned, otherwise the centroids are returned.
+     
+     If no initial positions are specified, they are chosen at random within
+     the data range.
+ 
+     Beware: it's common for k-means to return less than 'k' classes in
+     degenerate cases.
+     """
+     x = [ float(d) for d in data if d is not None]
+     x.sort()
+     N = len(x)
+     xinc = min([ x[i+1] - x[i] for i in range(N-1) if x[i+1] != x[i] ])
+     threshold = N * xinc / 2
+     mvt = 1e38
+     if init is None:
+         k -= 1 # k classes, so k-1 medoids
+         r = random.Random()
+         r.seed()
+         medoid = [ r.random() * ( max(x)-min(x) ) + min(x) for i in range(k) ]
+     else:
+         medoid = init
+         k = len(init)
+         if k < 2:
+             return []
+     while( mvt > threshold ):
+         old_medoid = medoid[:]
+         sums = [0]*k
+         n = [0]*k
+         for i in range(N):
+             dist = max(x) + 1
+             nearest = None
+             for m in range(k):
+                 dmxi = distance(medoid[m], x[i])
+                 if dmxi < dist:
+                     nearest = m
+                     dist = dmxi
+             sums[nearest] += x[i]
+             n[nearest] += 1
+         medoid = [ sums[i] / ( n[i] + 1 ) for i in range(k) ]
+         mvt = sum([ abs(mnew-mold) for (mnew,mold) in zip(medoid, old_medoid)])
+ 
+     # uniqueify medoids
+     unique = {}
+     for m in medoid:
+         unique[m] = 1
+     um = unique.keys()
+     um.sort()
+     if output != "breaks":
+         return um
+     if um[0] > min(x):
+         um = [ min(x) ]+um
+     if um[-1] == max(x):
+         um = um[:-1]
+     # find break values ("medians" in the sense of the 'distance' function)
+     out = [min(x)]
+     mi = 1
+     xim1 = min(x)
+     for xi in x:
+         if mi == len(um): break
+         if distance(xi, um[mi]) < distance(xi, um[mi-1]):
+             out.append( (xi+xim1) / 2.0 )
+             mi += 1
+         xim1 = xi
+     return out + [max(x)]
+ 
+ 
+ def classification_quality( breaks, data ):
+     """ Classification quality in the "Table" sense of Jenks & Caspall (1971).
+     
+     That is: a classification is better if the classes it produces fit the
+     original distribution in terms of sum of errors within each class with
+     respect to the mean value for the class. The original distribution has
+     zero error. Maximal error is obtained for a one-class classification
+     assigned the mean value of all observations.
+ 
+     This function can be used to evaluate classification quality in the "Area"
+     sense of the same authors if the data are weighted with the (projected)
+     area of the corresponding polygons.
+ 
+     We have adapted the measure to penalize one-observation classes.
+     """
+     data.sort()
+     cvalues = []
+     thrdown = min(data) - 1
+     ss = 0
+     nn = 0
+     for thrup in breaks:
+         s = 0
+         n = 0
+         for d in data:
+            if thrdown < d <= thrup:
+                s += d
+                n += 1
+         if n == 0:
+             cvalues.append(0)
+         else:
+             cvalues.append( s/float(n) )
+         ss += s
+         nn += n
+         thrdown = thrup
+     overall_mean = float(ss) / nn
+     thrdown = min(data) - 1
+     i = 0
+     evalues = []
+     for thrup in breaks:
+         s = 0
+         n = 0
+         for d in data:
+            if thrdown < d <= thrup:
+                s += abs( d - cvalues[i] )
+                n += 1
+         if n <= 1: # give less weight to one-case classes
+             evalues.append(overall_mean)
+         else:
+             try:
+                 evalues.append( s/float(n) / cvalues[i] * overall_mean )
+             except ZeroDivisionError:
+                 evalues.append( s/float(n) )
+         i += 1
+         thrdown = thrup
+     if overall_mean != 0:
+         return 1.0 - ( sum(evalues) / len(breaks) ) / overall_mean
+     else: 
+         return 1.0 - ( sum(evalues) / len(breaks) )
+ 
+ 
+ def classify_annealing( data, k=5, output="breaks",
+                         decay=.999, trials=30, distance=diffabs,
+                         objective=classification_quality, init=None ):
+     """A simulated annealing classifier.
+ 
+     A vector of centroids is annealed by decreasing-amplitude random
+     perturbation of their position. Centroids never move beyond the median
+     positions with respect to neighbours, so initial conditions are very
+     sensitive.
+ 
+     Annealing parameters are:
+ 
+     - 'k' -- the number of classes considered
+ 
+     - 'decay' -- the decay rate of temperature (amplitude of movement) for each
+     generation 
+ 
+     - 'trials' -- random perturbation is tried 'trials' times for each
+     temperature; only the best movement is kept
+ 
+     - 'distance' -- a two-parameter function defining the distance between two
+     centroids
+ 
+     - 'objective' -- the objective function to be maximized; parameters are
+     candidate vector and data vector
+ 
+     - 'init' -- initial candidate vector
+ 
+     With output="breaks", the function returns the breakpoints for
+     classification; otherwise it returns the best candidate vector of centroids
+     at the end of the annealing process.
+     """
+     r = random.Random()
+     r.seed()
+     x = [ float(d) for d in data if d is not None]
+     x.sort()
+     xmin = x[0]
+     xmax = x[-1]
+     N = len(x)
+     xinc = min([ x[i+1] - x[i] for i in range(N-1) if x[i+1] != x[i] ])
+     k -= 1
+     if init is None:
+         current = [ float(i)/k * (xmax-xmin) + xmin for i in range(k) ]
+     else:
+         current = init
+         k = len(init)
+         if k < 2:
+             return [] # XXX Should raise something...
+     current.insert(0, x[0])
+     spanminus = [ ( current[i+1] - current[i] ) / 2 for i in range(k-1) ]
+     spanplus = spanminus[:]
+     spanminus.insert(0, current[0] - xmin)
+     spanplus.append( xmax-current[-1] )
+ 
+     bestquality = 0
+     mvtrange = 1
+     while(1):
+         previous_bestquality = bestquality
+         bestquality = 0
+         best = None
+         for i in range(trials):
+             candidate = []
+             for j in range(k):
+                 if r.random() < .5:
+                     candidate.append( current[j] +
+                                       r.random() * spanplus[j] * mvtrange )
+                 else:
+                     candidate.append( current[j] -
+                                       r.random() * spanminus[j] * mvtrange )
+             quality = objective( candidate, x )
+             if quality > bestquality:
+                 bestquality = quality
+                 best = candidate
+         current = best
+         spanminus = [ ( current[i+1] - current[i] ) / 2 for i in range(k-1) ]
+         spanplus = spanminus[:]
+         spanminus.insert(0, current[0] - xmin)
+         spanplus.append( xmax - current[-1] )
+         mvtrange *= decay
+         
+         if mvtrange < xinc or bestquality == previous_bestquality:
+             break
+ 
+     out = []
+     unique = {}
+     for m in current:
+         unique[m] = 1
+         um = unique.keys()
+     um.sort()
+     if output == "medoids":
+         return um
+     if um[0] > min(x):
+         um = [ min(x) ] + um
+     if um[-1] == max(x):
+         um = um[:-1]
+     # find break values ("medians" in the sense of the 'distance' function)
+     out = [ min(x) ]
+     mi = 1
+     xim1 = min(x)
+     for xi in x:
+         if mi == len(um): break
+         if distance(xi, um[mi]) < distance(xi, um[mi-1]):
+             out.append( (xi+xim1) / 2.0 )
+             mi += 1
+         xim1 = xi
+     return out + [ max(x) ]
+ 
+ 
  class CustomRamp:
  
!     def __init__(self, prop1, prop2, interpolate_colors = interpolate_rgb):
          """Create a ramp between prop1 and prop2."""
          self.prop1 = prop1
          self.prop2 = prop2
+         self._interpolate = interpolate_colors
  
      def GetRamp(self):
          """Return this ramp."""
***************
*** 336,348 ****
              raise ValueError(_("invalid index"))
  
          newProps = ClassGroupProperties()
! 
!         self.__SetProperty(self.prop1.GetLineColor(),
!                            self.prop2.GetLineColor(),
!                            index, newProps.SetLineColor)
!         self.__SetProperty(self.prop1.GetFill(), self.prop2.GetFill(),
!                            index, newProps.SetFill)
! 
          w = (self.prop2.GetLineWidth() - self.prop1.GetLineWidth()) \
              * index \
              + self.prop1.GetLineWidth()
--- 787,798 ----
              raise ValueError(_("invalid index"))
  
          newProps = ClassGroupProperties()
!         newProps.SetLineColor( self._interpolate( self.prop1.GetLineColor(),
!                                                   self.prop2.GetLineColor(),
!                                                   index ) )
!         newProps.SetFill( self._interpolate( self.prop1.GetFill(),
!                                              self.prop2.GetFill(),
!                                              index ) )
          w = (self.prop2.GetLineWidth() - self.prop1.GetLineWidth()) \
              * index \
              + self.prop1.GetLineWidth()
***************
*** 355,388 ****
  
          return newProps
  
-     def __SetProperty(self, color1, color2, index, setf):
-         """Use setf to set the appropriate property for the point
-         index percent between color1 and color2. setf is a function
-         to call that accepts a Color object or Transparent.
-         """
- 
-         if color1 is Transparent and color2 is Transparent:
-             setf(Transparent)
-         elif color1 is Transparent:
-             setf(Color(
-                  color2.red   * index,
-                  color2.green * index,
-                  color2.blue  * index))
-         elif color2 is Transparent:
-             setf(Color(
-                  color1.red   * index,
-                  color1.green * index,
-                  color1.blue  * index))
-         else:
-             setf(Color(
-                 (color2.red   - color1.red)   * index + color1.red,
-                 (color2.green - color1.green) * index + color1.green,
-                 (color2.blue  - color1.blue)  * index + color1.blue))
- 
  class MonochromaticRamp(CustomRamp):
      """Helper class to make ramps between two colors."""
  
!     def __init__(self, start, end):
          """Create a Monochromatic Ramp.
  
          start -- starting Color
--- 805,814 ----
  
          return newProps
  
  class MonochromaticRamp(CustomRamp):
      """Helper class to make ramps between two colors."""
  
!     def __init__(self, start, end, interpolate_colors = interpolate_rgb):
          """Create a Monochromatic Ramp.
  
          start -- starting Color
***************
*** 397,403 ****
          ep.SetLineColor(end)
          ep.SetFill(end)
  
!         CustomRamp.__init__(self, sp, ep)
  
  grey_ramp         = MonochromaticRamp(Color(1, 1, 1),  Color(0, 0, 0))
  red_ramp          = MonochromaticRamp(Color(1, 1, 1),  Color(.8, 0, 0))
--- 823,829 ----
          ep.SetLineColor(end)
          ep.SetFill(end)
  
!         CustomRamp.__init__(self, sp, ep, interpolate_colors)
  
  grey_ramp         = MonochromaticRamp(Color(1, 1, 1),  Color(0, 0, 0))
  red_ramp          = MonochromaticRamp(Color(1, 1, 1),  Color(.8, 0, 0))



*** thuban/Thuban/UI/classgen.py	2003-08-08 08:01:30.000000000 -0500
--- thuban-dcabkp/Thuban/UI/classgen.py	2004-02-26 18:49:33.000000000 -0500
***************
*** 1,6 ****
--- 1,7 ----
  # Copyright (c) 2003 by Intevation GmbH
  # Authors:
  # Jonathan Coles <jonathan@intevation.de>
+ # Daniel Calvelo <dca@users.sourceforge.net>
  #
  # This program is free software under the GPL (>=v2)
  # Read the file COPYING coming with Thuban for details.
***************
*** 18,24 ****
  
  from Thuban import _
  
! from Thuban.Model.classification import ClassGroupProperties
  
  from Thuban.Model.table import FIELDTYPE_INT, FIELDTYPE_DOUBLE, \
       FIELDTYPE_STRING
--- 19,25 ----
  
  from Thuban import _
  
! from Thuban.Model.classification import ClassGroupProperties, ClassGroupRange, ClassGroupSingleton
  
  from Thuban.Model.table import FIELDTYPE_INT, FIELDTYPE_DOUBLE, \
       FIELDTYPE_STRING
***************
*** 27,36 ****
  from Thuban.Model.range import Range
  from Thuban.UI.common import ThubanBeginBusyCursor, ThubanEndBusyCursor
  
  import classifier, resource
  
! from Thuban.Model.classgen import \
!     generate_uniform_distribution, generate_singletons, generate_quantiles, \
      CustomRamp, grey_ramp, red_ramp, green_ramp, blue_ramp, green_to_red_ramp, \
      HotToColdRamp, FixedRamp
  
--- 28,43 ----
  from Thuban.Model.range import Range
  from Thuban.UI.common import ThubanBeginBusyCursor, ThubanEndBusyCursor
  
+ from Thuban.Model.histogram import bincount, histcount
+ from Thuban.UI.graph import BinHistogram, PointDensityGram
+ from Thuban.UI.sliders import Sliders
+ 
  import classifier, resource
  
! from Thuban.Model.color import interpolate_rgb, interpolate_hsv
! from Thuban.Model.classgen import generate_ranges_from_list,\
!     generate_uniform_distribution, generate_singletons, generate_quantiles,\
!     classify_discont, classify_kmeans, classify_annealing, \
      CustomRamp, grey_ramp, red_ramp, green_ramp, blue_ramp, green_to_red_ramp, \
      HotToColdRamp, FixedRamp
  
***************
*** 43,48 ****
--- 50,59 ----
  GENCOMBOSTR_UNIFORM = _("Uniform Distribution")
  GENCOMBOSTR_UNIQUE = _("Unique Values")
  GENCOMBOSTR_QUANTILES = _("Quantiles from Table")
+ GENCOMBOSTR_CUSTOM = _("Custom Breaks")
+ GENCOMBOSTR_DISCONT = _("Distribution Discontinuities")
+ GENCOMBOSTR_KMEANS = _("Distribution Classes")
+ GENCOMBOSTR_KANNEAL = _("Optimized Distribution Classes (slow)")
  
  PROPCOMBOSTR_CUSTOM     = _("Custom Ramp")
  PROPCOMBOSTR_GREY       = _("Grey Ramp")
***************
*** 97,102 ****
--- 108,117 ----
          if self.type in (FIELDTYPE_INT, FIELDTYPE_DOUBLE):
              self.genpanels.append((GENCOMBOSTR_UNIFORM, GenUniformPanel))
              self.genpanels.append((GENCOMBOSTR_QUANTILES, GenQuantilesPanel))
+             self.genpanels.append((GENCOMBOSTR_CUSTOM, GenCustomPanel))
+             self.genpanels.append((GENCOMBOSTR_DISCONT, GenDiscontPanel))
+             self.genpanels.append((GENCOMBOSTR_KMEANS, GenKmeansPanel))
+             self.genpanels.append((GENCOMBOSTR_KANNEAL, GenKannealPanel))
  
          for name, clazz in self.genpanels:
              self.genChoice.Append(name, [clazz, None])
***************
*** 220,226 ****
  
          if genSel in (GENCOMBOSTR_UNIFORM,          \
                        GENCOMBOSTR_UNIQUE,           \
!                       GENCOMBOSTR_QUANTILES):
  
              numGroups = genPanel.GetNumGroups()
  
--- 235,245 ----
  
          if genSel in (GENCOMBOSTR_UNIFORM,          \
                        GENCOMBOSTR_UNIQUE,           \
!                       GENCOMBOSTR_QUANTILES,        \
!                       GENCOMBOSTR_CUSTOM,           \
!                       GENCOMBOSTR_DISCONT,          \
!                       GENCOMBOSTR_KMEANS,           \
!                       GENCOMBOSTR_KANNEAL):
  
              numGroups = genPanel.GetNumGroups()
  
***************
*** 282,287 ****
--- 301,353 ----
                  else:
                      self.parent._SetClassification(self.clazz)
  
+             elif genSel == GENCOMBOSTR_CUSTOM:
+ 
+                 amin = genPanel.GetMin()
+                 amax = genPanel.GetMax()
+                 precision = genPanel.GetPrecision()
+                 raw_list = [ p * (amax - amin) + amin
+                              for p in genPanel.sliders.GetList()]
+                 _list = [ precision * int( x / precision ) for x in raw_list ]
+                 if _list[0]>amin:
+                     _list = [amin]+_list
+                 if _list[-1]<amax:
+                     _list.append( amax )
+                 self.clazz = generate_ranges_from_list( _list, ramp)
+                 self.parent._SetClassification(self.clazz)
+ 
+             elif genSel == GENCOMBOSTR_DISCONT:
+                 
+                 genPanel.breaks = classify_discont(
+                     genPanel._list,
+                     maxclasses = genPanel.spin_nclasses.GetValue()
+                     )
+                 self.clazz = generate_ranges_from_list(
+                     genPanel.breaks,
+                     ramp)
+                 self.parent._SetClassification(self.clazz)
+                 
+             elif genSel == GENCOMBOSTR_KMEANS:
+ 
+                 self.clazz = generate_ranges_from_list(
+                     classify_kmeans( genPanel.GetList(),
+                                      k = genPanel.kSpinner.GetValue()
+                                      ),
+                     ramp)
+                 self.parent._SetClassification(self.clazz)
+ 
+             elif genSel == GENCOMBOSTR_KANNEAL:
+                 self.clazz = generate_ranges_from_list(
+                     classify_annealing(
+                     genPanel.GetList(),
+                     k = genPanel.kSpinner.GetValue(),
+                     init = classify_kmeans( genPanel.GetList(),
+                                             k = genPanel.kSpinner.GetValue(),
+                                             output = "medoids")
+                     ), ramp)
+                 self.parent._SetClassification(self.clazz)
+                 
+ 
      def OnCancel(self, event):
          self.Close()
  
***************
*** 906,917 ****
--- 972,1308 ----
              finally:
                  ThubanEndBusyCursor()
  
+ 
+ class GenCustomPanel(wxPanel):
+ 
+     def __init__(self, parent, layer, fieldName, fieldType):
+         wxPanel.__init__(self, parent, -1)
+ 
+         default_nbins = 10
+         self.parent = parent
+         self.layer = layer
+         self.fieldName = fieldName
+         self.fieldType = fieldType
+ 
+         topSizer = wxStaticBoxSizer(wxStaticBox(self, -1, _("Frequency count")), 
+                                     wxVERTICAL)
+ 
+         spinSizer = wxBoxSizer( wxHORIZONTAL )
+         spinSizer.Add( wxStaticText(self, -1, _("Number of bins:")), 0, wxALL )
+         self.spin_numBins = wxSpinCtrl(self, -1, style=wxTE_RIGHT)
+         self.spin_numBins.SetRange(2, sys.maxint)
+         self.spin_numBins.SetValue( default_nbins )
+         EVT_SPINCTRL(self, self.spin_numBins.GetId(), self.OnChangeNBins)
+ 
+         spinSizer.Add( self.spin_numBins, 0, wxALL )
+         topSizer.Add( spinSizer, 0, wxALIGN_RIGHT|wxALIGN_TOP)
+ 
+         curdata = self.GetList()
+         self._min = min(curdata)
+         self._max = max(curdata)
+         
+         self.hist = BinHistogram(self, data=curdata, bins=default_nbins )
+         topSizer.Add( self.hist, 0, wxALIGN_TOP|wxALIGN_CENTER_HORIZONTAL )
+ 
+         self.datashow = PointDensityGram(self, data=curdata)
+         topSizer.Add(self.datashow, 0, wxALIGN_TOP|wxALIGN_CENTER_HORIZONTAL )
+ 
+         ThubanBeginBusyCursor()
+         # FIXME Ugly: this goes through parent's objects. Interfaces should be cleaned up a bit
+         # FIXME 2 : putting this in __init__ prevents update when another classification has been
+         #           generated but the dialog is not closed and reopened. This should be triggered
+         #           by the OnOk handler or the _OnGenType handler. Interfaces should be cleaned up...
+         current_breaks = []
+         current_clazz = self.parent.parent.classGrid.GetTable().GetClassification()
+         for i in range(current_clazz.GetNumGroups()):
+             aGroup = current_clazz.GetGroup( i )
+             if isinstance( aGroup, ClassGroupRange ):
+                 current_breaks.append( aGroup.GetMin() )
+                 current_breaks.append( aGroup.GetMax() )
+             elif isinstance( aGroup, ClassGroupSingleton ):
+                 current_breaks.append( aGroup.GetValue() )
+             else:
+                 pass
+         if len(current_breaks) == 0 or len(current_breaks) > 25:
+             # custom breaks become unusable after that...
+             initial_breaks = [ float(i+1) / 6 for i in range(5) ]
+         else:
+             unique = {}
+             for cb in current_breaks:
+                 if self._max > cb > self._min:
+                     unique[cb] = 1
+             initial_unscaled_breaks = unique.keys()
+             initial_unscaled_breaks.sort()
+             initial_breaks = [ ( b - self._min ) / ( self._max - self._min )
+                                for b in initial_unscaled_breaks ]
+             if initial_breaks[0] < 1e-4:
+                 initial_breaks = initial_breaks[1:]
+             if initial_breaks[-1] > 1 - 1e-4:
+                 initial_breaks = initial_breaks[:-1]
+         
+         self.sliders = Sliders(self, pos = initial_breaks)
+ 
+         topSizer.Add(self.sliders, 0, wxALIGN_TOP|wxALIGN_CENTER_HORIZONTAL )
+ 
+         self.SetSizer(topSizer)
+         topSizer.SetSizeHints( self )
+         self.SetAutoLayout(True)
+ 
+         # FIXME: this can be problematic for large tables.
+         # An option would be to calculate the precision from min(int(log10(x)))
+         # for all data x taking due care of FP precision problems.
+         l = self.GetList()
+         l.sort()
+         self._precision = min( [ l[i+1] - l[i]
+                                  for i in range(len(l) - 1)
+                                  if l[i+1] != l[i] ] )
+         ThubanEndBusyCursor()
+ 
+     def GetNumGroups(self):
+         return self.sliders.getNum()+1 #sliders are breaks, giving +1 intervals
+ 
+     def GetMin(self):
+         return self._min
+ 
+     def GetMax(self):
+         return self._max
+ 
+     def GetPrecision(self):
+         return self._precision
+ 
+     def GetList(self):
+         _list = []
+         table = self.layer.ShapeStore().Table()
+         if table is not None:
+             ThubanBeginBusyCursor()
+             for i in range(table.NumRows()):
+                 _list.append(table.ReadValue(i, self.fieldName))
+             ThubanEndBusyCursor()
+         return _list
+ 
+     def OnChangeNBins(self, event):
+         self.hist.nbins = self.spin_numBins.GetValue()
+         self.hist.OnPaint( None )
+ 
+ ID_DISCONT_BUTTON = 4100
+ ID_DISCONT_SPINNER = 4101
+ 
+ class GenDiscontPanel(wxPanel):
+     def __init__(self, parent, layer, fieldName, fieldType):
+         wxPanel.__init__(self, parent, -1)
+         self.parent = parent
+         self.layer = layer
+         ThubanBeginBusyCursor()
+         _table = self.layer.ShapeStore().Table()
+         _list = [ _table.ReadValue(i, fieldName)
+                   for i in range( _table.NumRows() ) ]
+         self._min = min(_list)
+         self._max = max(_list)
+         _list.sort()
+         self._precision = min( [ _list[i+1] - _list[i]
+                                  for i in range(len(_list) - 1)
+                                  if _list[i+1] != _list[i] ] )
+         self.resultlog = []
+         self.nclasses = None
+         self.breaks = classify_discont( _list, resultlog = self.resultlog )
+         self._list = _list
+         theSizer = wxStaticBoxSizer(wxStaticBox(
+             self, -1,
+             _("Got %d classes in the distribution of %s")
+             % (len(self.breaks) - 1, _table.NumRows())), 
+                                     wxVERTICAL)
+         theSizer.Add( wxStaticText(
+             self, -1,
+             _("\nThe breakpoints have been obtained\n"
+               "by recursively partitioning the data\n"
+               "distributions.\n\n"
+               "These breakpoints reproduce faithfully\n"
+               "the original distribution of the tabular data.\n"
+               "\nBased on code by J-P. Grimmeau, IGEAT-ULB "
+ 	      "(grimmeau@ulb.ac.be).\n\n")), 0, wxALL)
+         subSizer = wxBoxSizer( wxHORIZONTAL )
+         detailsButton = wxButton( self, ID_DISCONT_BUTTON, _("See Details") )
+         EVT_BUTTON( self, ID_DISCONT_BUTTON, self.OnGetDetails )
+         self.spin_nclasses = wxSpinCtrl(
+             self, ID_DISCONT_SPINNER, style = wxTE_RIGHT )
+         self.spin_nclasses.SetRange(2, len(self.breaks)-1 )
+         self.spin_nclasses.SetValue( len(self.breaks)-1 )
+         
+         subSizer.Add( detailsButton, 0, wxALL )
+         subSizer.Add( wxStaticText(
+             self, -1, _("  Fix number of classes ") ), 0, wxALIGN_RIGHT )
+         subSizer.Add( self.spin_nclasses, 0, wxALIGN_RIGHT )
+         theSizer.Add( subSizer, 0, wxALL )
+         self.SetSizer( theSizer )
+         self.SetAutoLayout(True)
+         ThubanEndBusyCursor()
+ 
+     def OnGetDetails(self, event):
+         # first, calculate format string 
+         from math import log
+         ndec = 1
+         f = self._precision * 10
+         while ndec < 10:
+             if abs( int( f ) - f ) < 1e-10 :
+                 break
+             ndec += 1
+             f *= 10
+         nlen = int( log( abs(self._max) ) / log(10) ) + 1 #integer part digits
+         nlen += 1 # decimal dot
+         nlen += ndec # fraction digits
+         fmt = "%"+repr(nlen)+"."+repr(ndec)+"f"
+         # second, generate HTML table for the above
+         linefmt = ("%4d ","%4d ","%6d ", fmt+" " , fmt+" ", fmt+" ", fmt+" ")
+         fmtlen = (5,5,7,nlen,nlen,nlen,nlen)
+         headers = (_("total # classes"),_("class"),_("# data points"),
+                    _("lower bound"),_("upper bound"),
+                    _("mean value"),_("standard deviation") )
+         title = _("Details from Distribution Discontinuities")
+         html = "<html><head><title>"
+         html += title
+         html += "</title></head>\n<body>\n"
+         html += "<table align=\"right\">\n"
+         html += "<tr>\n"
+         for h in headers:
+             html += "  <th><font color=\"#000055\">"+h+"</font></th>"
+         html+= "</tr>\n"
+         # third, open the dialog to display it all.
+         refnclass = self.resultlog[0][0]
+         for det in self.resultlog:
+             if det[0] != refnclass:
+                 html += "<tr><td>&nbsp;</td></tr>"
+                 refnclass = det[0]
+             html += "<tr>\n"
+             for i in range(len(det)):
+                 if det[i] is not None:
+                     html += "<td>"+( linefmt[i] % det[i] )+"</td>\n"
+                 else:
+                     html += "<td>&nbsp;</td>\n"
+             html += "</tr>\n"
+         html += "</table></body></html>\n"
+         
+         from wxPython.html import wxHtmlWindow
+         self.msgwindow = wxDialog( self, -1, title, style = wxDIALOG_NO_PARENT|wxDEFAULT_DIALOG_STYLE )
+         self.msgwindow.SetSize( (600, 400) )
+         htmlSizer = wxBoxSizer( wxVERTICAL )
+         self.msg = wxHtmlWindow( self.msgwindow, -1, style = wxVSCROLL|wxHSCROLL )
+         self.msg.SetPage( html )
+         self.msg.SetSize( (600, 400) )
+         self.msg.Show()
+         htmlSizer.Add( self.msg, 0, wxALL|wxGROW )
+         self.msgwindow.SetSizer( htmlSizer )
+         self.msgwindow.SetAutoLayout( True )
+         self.msgwindow.Show()
+         self.msgwindow.Raise()
+         EVT_CLOSE( self.msgwindow, self._DestroyDetailsDialog )
+ 
+     def _DestroyDetailsDialog( self, event ):
+         self.msg.Close()
+         self.msgwindow.Destroy()
+         
+     def GetNumGroups(self):
+         return len(self.breaks)
+         
+     def GetMin(self):
+         return self._min
+ 
+     def GetMax(self):
+         return self._max
+ 
+     def GetList(self):
+         return self._list
+ 
+ class GenKmeansPanel(wxPanel):
+     def __init__(self, parent, layer, fieldName, fieldType):
+         wxPanel.__init__(self, parent, -1)
+         self.parent = parent
+         self.layer = layer
+         self.fieldName = fieldName
+         ThubanBeginBusyCursor()
+         _table = self.layer.ShapeStore().Table()
+         topSizer = wxStaticBoxSizer(wxStaticBox(
+             self, -1, _("Distribution Classes")), 
+                                     wxVERTICAL)
+ 
+         spinSizer = wxBoxSizer( wxHORIZONTAL )
+         spinSizer.Add( wxStaticText(
+             self, -1, _("Initial number of Classes:")), 0, wxALL )
+         self.kSpinner = wxSpinCtrl(self, -1, style=wxTE_CENTER)
+         self.kSpinner.SetRange(2, sys.maxint)
+         self.kSpinner.SetValue( 5 )
+ 
+         spinSizer.Add( self.kSpinner, 0, wxALL )
+         topSizer.Add( wxStaticText(
+             self, -1,
+             _("The classes will be generated by clustering\n"
+               "the tabular data using k-means.\n\n"
+               "Bear in mind that this does not take\n"
+               "into account the layer geometry.\n\n"
+               "It's likely you obtain fewer classes than\n"
+               "the initial number you specified.\n"
+               "Furthermore, the algorithm relies on random\n"
+               "initial conditions, thus you may obtain different\n"
+               "results on different runs\n\n")), 0, wxALL) 
+         topSizer.Add( spinSizer, 0, wxALIGN_RIGHT|wxALIGN_TOP)
+         
+         self.SetSizer( topSizer )
+         self.SetAutoLayout(True)
+         ThubanEndBusyCursor()
+ 
+     def GetNumGroups(self):
+         return self.kSpinner.GetValue()
+ 
+     def GetList(self):
+         _list = []
+         table = self.layer.ShapeStore().Table()
+         if table is not None:
+             ThubanBeginBusyCursor()
+             for i in range(table.NumRows()):
+                 _list.append(table.ReadValue(i, self.fieldName))
+             ThubanEndBusyCursor()
+         return _list
+ 
+ class GenKannealPanel(GenKmeansPanel):
+     def __init__(self, parent, layer, fieldName, fieldType):
+         wxPanel.__init__(self, parent, -1)
+         self.parent = parent
+         self.layer = layer
+         self.fieldName = fieldName
+         ThubanBeginBusyCursor()
+         _table = self.layer.ShapeStore().Table()
+         topSizer = wxStaticBoxSizer(wxStaticBox(
+             self, -1, _("Optimized Distribution Classes")), 
+                                     wxVERTICAL)
+ 
+         spinSizer = wxBoxSizer( wxHORIZONTAL )
+         spinSizer.Add( wxStaticText(
+             self, -1, _("Initial number of Classes:")), 0, wxALL )
+         self.kSpinner = wxSpinCtrl(self, -1, style=wxTE_CENTER)
+         self.kSpinner.SetRange(2, sys.maxint)
+         self.kSpinner.SetValue( 5 )
+         
+         spinSizer.Add( self.kSpinner, 0, wxALL )
+         topSizer.Add(wxStaticText(
+             self, -1,
+             _("You will initially obtain a k-means clustering of\n"
+               "your tabular data, from which an extra optimization\n"
+               "will be performed.\n\n"
+               "It's likely you will obtain fewer classes than\n"
+               "specified here.\n\n")), 0, wxALL) 
+         topSizer.Add( spinSizer, 0, wxALIGN_RIGHT|wxALIGN_TOP)
+         
+         self.SetSizer( topSizer )
+         self.SetAutoLayout(True)
+         ThubanEndBusyCursor()
+ 
+ 
  ID_CUSTOMRAMP_COPYSTART = 4001
  ID_CUSTOMRAMP_COPYEND = 4002
  ID_CUSTOMRAMP_EDITSTART = 4003
  ID_CUSTOMRAMP_EDITEND = 4004
  ID_CUSTOMRAMP_SPROP = 4005
  ID_CUSTOMRAMP_EPROP = 4006
+ ID_CUSTOMRAMP_MODECHOICE = 4007
  
  class CustomRampPanel(wxPanel):
  
***************
*** 919,925 ****
          wxPanel.__init__(self, parent, -1)
  
          topSizer = wxStaticBoxSizer(wxStaticBox(self, -1, ""), wxHORIZONTAL)
! 
          bsizer = wxBoxSizer(wxVERTICAL)
          bsizer.Add(wxStaticText(self, -1, _("Start:")), 0, wxALL | wxCENTER, 4)
          self.startPropCtrl = classifier.ClassGroupPropertiesCtrl(
--- 1310,1322 ----
          wxPanel.__init__(self, parent, -1)
  
          topSizer = wxStaticBoxSizer(wxStaticBox(self, -1, ""), wxHORIZONTAL)
! 	
! 	self.modeChoice = wxRadioBox( self, ID_CUSTOMRAMP_MODECHOICE, 
! 	                         _("Mode"),wxDefaultPosition,wxDefaultSize, 
! 				 [_("HSV Ramp"),_("RGB Ramp")],
! 				 1, wxRA_SPECIFY_COLS)
!         topSizer.Add( self.modeChoice, 0, wxLEFT|wxALIGN_CENTER_VERTICAL )
!       
          bsizer = wxBoxSizer(wxVERTICAL)
          bsizer.Add(wxStaticText(self, -1, _("Start:")), 0, wxALL | wxCENTER, 4)
          self.startPropCtrl = classifier.ClassGroupPropertiesCtrl(
***************
*** 968,974 ****
                        | wxALIGN_CENTER_HORIZONTAL \
                        | wxALIGN_CENTER_VERTICAL, 
                     4)
- 
          EVT_BUTTON(self, ID_CUSTOMRAMP_COPYSTART, self._OnCopyStart)
          EVT_BUTTON(self, ID_CUSTOMRAMP_COPYEND, self._OnCopyEnd)
          EVT_BUTTON(self, ID_CUSTOMRAMP_EDITSTART, self._OnEditStart)
--- 1365,1370 ----
***************
*** 979,986 ****
          topSizer.SetSizeHints(self)
  
      def GetRamp(self):
          return CustomRamp(self.startPropCtrl.GetProperties(),
!                           self.endPropCtrl.GetProperties())
  
      def _OnCopyStart(self, event):
          self.endPropCtrl.SetProperties(self.startPropCtrl.GetProperties())
--- 1375,1387 ----
          topSizer.SetSizeHints(self)
  
      def GetRamp(self):
+     	if self.modeChoice.GetStringSelection()==_("HSV Ramp"):
+ 	    interpolate_function = interpolate_hsv
+ 	else:
+ 	    interpolate_function = interpolate_rgb
          return CustomRamp(self.startPropCtrl.GetProperties(),
!                           self.endPropCtrl.GetProperties(),
!                           interpolate_function)
  
      def _OnCopyStart(self, event):
          self.endPropCtrl.SetProperties(self.startPropCtrl.GetProperties())



*** thuban/Thuban/Model/color.py	2003-08-05 07:38:52.000000000 -0500
--- thuban-dcabkp/Thuban/Model/color.py	2004-02-13 02:37:08.000000000 -0500
***************
*** 6,11 ****
--- 6,13 ----
  # This program is free software under the GPL (>=v2)
  # Read the file COPYING coming with Thuban for details.
  
+ from colorsys import rgb_to_hsv, hsv_to_rgb
+ 
  __version__ = "$Revision: 1.10 $"
  
  class Color:
***************
*** 58,60 ****
--- 60,163 ----
  Transparent = _Transparent()
  Black = Color(0, 0, 0)
  
+ def interpolate_rgb(color1, color2, weight):
+     """Return a weighted average of color1 and color2 in RGB colorspace.
+ 
+     color1 and color2 must be instances of Color.
+     For components R,G,B in [0,1], the weighted average is
+ 
+     color1 + weight * ( color2 - color1 )
+ 
+     so for weight==0, you have color1 and for weight==1 color2.
+     """
+     if color1 is Transparent and color2 is Transparent:
+         return Transparent
+     elif color1 is Transparent:
+         return Color(
+             color2.red   * weight,
+             color2.green * weight,
+             color2.blue  * weight)
+     elif color2 is Transparent:
+         return Color(
+             color1.red   * weight,
+             color1.green * weight,
+             color1.blue  * weight)
+     else:
+         return Color(
+             (color2.red   - color1.red)   * weight + color1.red,
+             (color2.green - color1.green) * weight + color1.green,
+             (color2.blue  - color1.blue)  * weight + color1.blue)
+ 
+ def interpolate_hsv( color1, color2, weight ):
+     """Return a weighted average of color1 and color2 in HSV colorspace.
+ 
+     For each C from H,S,V in [0,1], the weighted average is
+ 
+     C1 + weight * ( C2 - C1 )
+ 
+     so for weight==0, you have color1 and for weight==1 color2.
+     """
+     if color1 is Transparent and color2 is Transparent:
+         return Transparent
+     elif color1 is Transparent:
+         return Color(
+             color2.red   * weight,
+             color2.green * weight,
+             color2.blue  * weight)
+     elif color2 is Transparent:
+         return Color(
+             color1.red   * weight,
+             color1.green * weight,
+             color1.blue  * weight)
+     else:
+         #
+         # This is the interesting part
+         #
+         h1,s1,v1 = rgb_to_hsv( color1.red, color1.green, color1.blue )
+         h2,s2,v2 = rgb_to_hsv( color2.red, color2.green, color2.blue )
+         #
+         # Special case 1: grays are given red hue by default, although
+         #                 they do not have hue proper.
+         #                 This can trick the interpolation: going from
+         #                 white to blue would give pale purple shades
+         #
+         saturation_threshold = 1e-5 # Let's have some flexibility here...
+         if s1 <= saturation_threshold and s2 > saturation_threshold:
+             # Color 2 really has a hue
+             h1 = h2
+         if s2<=saturation_threshold and s1 > saturation_threshold:
+             h2 = h1
+             # All other cases are Ok, including s1==s2==0
+             
+         #
+         # Special case 2: hue is a ring, you don't go from yellow to
+         #                 violet through green, you move directly from
+         #                 yellow through orange and red to violet.
+         #                 Or at least it seems more intuitive that way.
+         #
+         if abs(h2 - h1) > .5: # we go "the other way around"
+             if h1 > h2:
+                 h2 = h2 + 1.0
+             else:
+                 h1 = h1 + 1.0
+             hh = h1 + weight * (h2 - h1)
+             hh = hh - int(hh) # equivalent to hh mod 1
+         else:
+             hh = h1 + weight * (h2 - h1)
+ 
+         #
+         # Special case 3: black is unsaturated by default, so going from
+         #                 black to green goes through greenish dark grays
+         #                 instead of dark greens
+         #
+         value_threshold = 1e-5 
+         if v1 <= value_threshold and v2 > value_threshold: # Color 1 is black
+             s1 = s2
+         if v2 <= value_threshold and v1 > value_threshold:
+             s2 = s1
+             
+         r,g,b = hsv_to_rgb( hh ,
+                             s1 + weight * (s2 - s1) ,
+                             v1 + weight * (v2 - v1) )
+         return Color(r,g,b)
+ 