[Gzz-commits] gzz/lava/gfx/shopt opt.py smcurv.py

[Tuomas J. Lukka]

CVSROOT:        /cvsroot/gzz
Module name:    gzz
Changes by:     Tuomas J. Lukka <address@hidden>        02/08/28 08:00:44

Modified files:
        lava/gfx/shopt : opt.py smcurv.py 

Log message:
        Some further work; conjugate gradient is better

CVSWeb URLs:
http://savannah.gnu.org/cgi-bin/viewcvs/gzz/gzz/lava/gfx/shopt/opt.py.diff?tr1=1.3&tr2=1.4&r1=text&r2=text
http://savannah.gnu.org/cgi-bin/viewcvs/gzz/gzz/lava/gfx/shopt/smcurv.py.diff?tr1=1.2&tr2=1.3&r1=text&r2=text

Patches:
Index: gzz/lava/gfx/shopt/opt.py
diff -c gzz/lava/gfx/shopt/opt.py:1.3 gzz/lava/gfx/shopt/opt.py:1.4
*** gzz/lava/gfx/shopt/opt.py:1.3       Sat Aug 24 01:51:42 2002
--- gzz/lava/gfx/shopt/opt.py   Wed Aug 28 08:00:44 2002
***************
*** 1,5 ****
--- 1,6 ----
  
  from Numeric import *
+ import RandomArray
  from Tkinter import *
  import math
  import os
***************
*** 40,69 ****
        geom += "%s %s %s\n" % (x[i], y[i], z[i])
      return geom
  
  
  class MomGrad:
      def __init__(self, func, x):
        self.func = func
        self.mom = x * 0
        self.x = x
      def round(self, step, brake):
        (u, g) = self.func(self.x)
        self.u = u
!       self.mom += g
!       l = add.reduce(ravel(self.mom*self.mom))
!       stp = -step / sqrt(l) * self.mom
!       self.x += stp
!       # self.x = maximum(self.x, 0)
        self.mom *= brake
        # print "Mom, l:", self.mom, l
        return self.u
  
  
  
! def slider(obj, name, var, default, from_, to, resolution=0.001):
      setattr(obj, name, var)
      var.set(default)
!     Scale(obj, label=name, orient="horizontal", 
        resolution=resolution,
        from_=from_, to=to, variable=var).pack(
            fill="both", expand=1)
--- 41,145 ----
        geom += "%s %s %s\n" % (x[i], y[i], z[i])
      return geom
  
+ def flatsum(x):
+     return add.reduce(ravel(x))
+ 
+ class LineMin:
+     def __init__(self):
+       self.l = 0.01
+     def linemin(self, func, x, dx):
+       (u, g) = func(x)
+       origu = u
+       origg = g
+       # print "Start linemin at ",u
+       dot = -flatsum(g * dx)
+       if dot < 0: 
+           print "Inverting: ",dot
+           dx = dx * -1
+           dot = -dot
+       l = 0
+       while dot > 0:
+           prevl = l
+           l += self.l
+           l *= 2
+           (uc, gc) = func(x + l * dx)
+           dot = -flatsum(gc * dx)
+           # print "Lextend: ",l,dot
+       curmax = l     # Point where gradient dot was known to be positive
+       curmin = prevl # Point where gradient dot was known to be negative
+       while (curmax-curmin)/(curmax+curmin) > 0.05:
+           testl = 0.5*(curmax+curmin)
+           (uc, gc) = func(x + testl * dx)
+           # print "Lloop: ",curmin,curmax,uc
+           dot = -flatsum(gc * dx)
+           if dot < 0:
+               curmax = testl
+           else:
+               curmin = testl
+       final = (0.5 * (curmax +curmin))
+       (tu, tg) = func(x + final * dx)
+       self.l = final
+       if tu > (u * 1.03): # Allow small jiggle to help out of local minima
+           print "ARGH! LINE SEARCH FAILURE: ",u, tu, self.l, final, g, tg
+           return x, origu, origg
+       return x + final * dx, tu, tg
+ 
+ class ConjGrad:
+     def __init__(self, func, x):
+       self.func = func
+       self.x = x
+       self.lin = LineMin()
+       (u,g) = func(x)
+       self.grad = g
+       self.dir = -g
+     def round(self, step, brake):
+       (self.x, up, gp) = self.lin.linemin(self.func, self.x, self.dir)
+       # Polak-Ribiere
+       beta = flatsum(gp * (gp-self.grad)) / flatsum(self.grad**2);
+       self.dir = -gp + beta * self.dir
+       self.grad = gp
+       return up
+     def kick(self):
+       (u,g) = self.func(self.x)
+       self.x += 0.01 * (RandomArray.random(g.shape)-0.5) * g
+       (u,g) = self.func(self.x)
+       self.grad = g
+       self.dir = -g
+       
+ 
+ 
+       
  
  class MomGrad:
      def __init__(self, func, x):
        self.func = func
        self.mom = x * 0
        self.x = x
+       self.lin = LineMin()
      def round(self, step, brake):
+ 
        (u, g) = self.func(self.x)
        self.u = u
!       self.mom -= g
!       # l = flatsum(self.mom*self.mom)
!       # stp = -step / sqrt(l) * self.mom
!       # self.x += stp
        self.mom *= brake
+ 
+       self.x, foo = self.lin.linemin(self.func, self.x, self.mom)
        # print "Mom, l:", self.mom, l
        return self.u
+     def kick(self):
+       (u,g) = self.func(self.x)
+       self.x += 0.01 * (RandomArray.random(g.shape)-0.5) * g
  
  
  
! def slider(obj, name, var, default, from_, to, resolution=0.0001, 
parent=None):
!     if parent == None: parent = obj
      setattr(obj, name, var)
      var.set(default)
!     Scale(parent, label=name, orient="horizontal", 
        resolution=resolution,
        from_=from_, to=to, variable=var).pack(
            fill="both", expand=1)
***************
*** 75,96 ****
        Button(self, text="RESET", command=self.reset).pack()
        self.pack(expand=1, fill="both")
        Button(self, text="Start", command = self.idle).pack()
        # Button(self, text="Kick", command = self.stretch.kick).pack()
!       slider(self, "dice", IntVar(), 10, 1, 100, resolution=1);
!       slider(self, "niter", IntVar(), 5, 1, 50, resolution=1);
  
        slider(self, "mom_step", DoubleVar(), 0.03, 0, 0.2)
        slider(self, "mom_brake", DoubleVar(), 0.97, 0.8, 1.0)
  
        self.pot = DoubleVar(); 
        Label(self, textvariable=self.pot).pack()
        
        self.reset()
        self.i = 0
  
      def reset(self):
        self.problem.reset(self.dice.get())
!       self.opt = MomGrad(self.problem, self.problem.defaultX())
      def idle(self):
        # print "IDLE"
        u = self.opt.round(step = self.mom_step.get(), brake = 
self.mom_brake.get())
--- 151,178 ----
        Button(self, text="RESET", command=self.reset).pack()
        self.pack(expand=1, fill="both")
        Button(self, text="Start", command = self.idle).pack()
+       Button(self, text="kick", command = self.kick).pack()
        # Button(self, text="Kick", command = self.stretch.kick).pack()
!       slider(self, "dice", IntVar(), 20, 1, 100, resolution=1);
!       slider(self, "niter", IntVar(), 1, 1, 50, resolution=1);
  
        slider(self, "mom_step", DoubleVar(), 0.03, 0, 0.2)
        slider(self, "mom_brake", DoubleVar(), 0.97, 0.8, 1.0)
  
+       self.problem.addWidgets(self)
+ 
        self.pot = DoubleVar(); 
        Label(self, textvariable=self.pot).pack()
+ 
        
        self.reset()
        self.i = 0
  
      def reset(self):
        self.problem.reset(self.dice.get())
!       self.opt = ConjGrad(self.problem, self.problem.defaultX())
!     def kick(self):
!       self.opt.kick()
      def idle(self):
        # print "IDLE"
        u = self.opt.round(step = self.mom_step.get(), brake = 
self.mom_brake.get())
Index: gzz/lava/gfx/shopt/smcurv.py
diff -c gzz/lava/gfx/shopt/smcurv.py:1.2 gzz/lava/gfx/shopt/smcurv.py:1.3
*** gzz/lava/gfx/shopt/smcurv.py:1.2    Sat Aug 24 14:39:31 2002
--- gzz/lava/gfx/shopt/smcurv.py        Wed Aug 28 08:00:44 2002
***************
*** 9,32 ****
  
  class SmoothCurve:
      def addWidgets(self, parent):
!       pass
      def reset(self, dice):
        self.dice = dice
      def defaultX(self):
        step = 1.0 / self.dice
!       rang = arrayrange(0, 1 + step, step)
        self.l = rang[NewAxis, :]
  
!       self.matrix = cos(math.pi * rang * (arange(rang.shape[0]) [:, NewAxis] 
))
! 
!       return concatenate( (0.1*self.l, 0.1 * self.l * self.l) )
      def __call__(self, v):
        # print "V:",v
  
        # def scalarsinglepot(x):
        def scalarneighpot(x1, y1, x2, y2):
            # return 1 * (sqrt((x2-x1)**2 + (y2-y1)**2)/self.dice - 1.0)**2 / 
self.dice
!           return 200 * ((x2-x1)**2 + (y2-y1)**2)
        def scalardneighpot(x1, y1, x2, y2, x3, y3):
            # Calculate angle of turning
            #return 0*x3
--- 9,52 ----
  
  class SmoothCurve:
      def addWidgets(self, parent):
!       frame = Frame(parent)
!       slider(self, "eqdist", DoubleVar(), 0.1, 0, 1, parent=frame)
!       slider(self, "kdist", DoubleVar(), 1, 0, 10, parent=frame)
!       slider(self, "distexp", DoubleVar(), 2, 0, 10, parent=frame, 
resolution=2)
!       slider(self, "kcurv", DoubleVar(), 1, 0, 10, parent=frame)
!       slider(self, "curvexp", DoubleVar(), 1, 0, 10, parent=frame)
!       slider(self, "xoffs", DoubleVar(), 1, 0, 5, parent=frame)
!       frame.pack(expand=1, fill="both")
! 
      def reset(self, dice):
        self.dice = dice
      def defaultX(self):
        step = 1.0 / self.dice
!       rang = arrayrange(0, 1 + self.dice) / (self.dice + 0.0)
        self.l = rang[NewAxis, :]
  
!       x = concatenate( (0.1*self.l, 0.05 + 0.0 * self.l ) )
!       x[1,0] = 0
!       x[1,-1] = 0.1
!       x[0,1] = 0.0
!       x[0,-2] = 0.1
!       return x
      def __call__(self, v):
        # print "V:",v
+       eqdist = self.eqdist.get()
+       kdist = self.kdist.get()
+       kcurv = self.kcurv.get()
+       xoffs = self.xoffs.get()
+       distexp = self.distexp.get()
+       curvexp = self.curvexp.get()
+       v[0,-1] = xoffs
+       v[0,-2] = xoffs
  
        # def scalarsinglepot(x):
        def scalarneighpot(x1, y1, x2, y2):
            # return 1 * (sqrt((x2-x1)**2 + (y2-y1)**2)/self.dice - 1.0)**2 / 
self.dice
!           dist = sqrt((x2-x1)**2 + (y2-y1)**2)
!           return kdist * ((dist*self.dice-eqdist)**distexp)/self.dice
        def scalardneighpot(x1, y1, x2, y2, x3, y3):
            # Calculate angle of turning
            #return 0*x3
***************
*** 34,43 ****
            dx2 = (x3-x2)
            dy1 = (y2-y1)
            dy2 = (y3-y2)
!           l1 = sqrt(dx1**2+dy1**2)
!           l2 = sqrt(dx2**2+dy2**2)
!           dot = (dx1*dx2 + dy1*dy2)/l1/l2
!           return ((1-dot))  
  
        grad = v*0
  
--- 54,65 ----
            dx2 = (x3-x2)
            dy1 = (y2-y1)
            dy2 = (y3-y2)
!           l1sq = (dx1**2+dy1**2)
!           l2sq = (dx2**2+dy2**2)
!           dot = (dx1*dx2 + dy1*dy2)
!           cosangle = dot / (sqrt(l1sq * l2sq+0.00001))
!           return kcurv * (1 - cosangle) ** curvexp
!           # (((l1sq+l2sq-2*dot))) 
  
        grad = v*0
  
***************
*** 78,85 ****
--- 100,116 ----
  
        grad[:,0] = 0
        grad[:,-1] = 0
+       grad[0,1] = 0
+       grad[0,-2] = 0
+       # Then, set the gradients so that the second and second-last points 
move towards
+       # the vertical line
+       # grad[0,1] = 10000 * (v[0,1]-v[0,0])
+ #     grad[0,-2] = 10000 * (v[0,-2]-v[0,-1])
+ #     artifu = 10000 * (0.5* (v[0,1]-v[0,0])**2 +  0.5*(v[0,-2]-v[0,-1])**2)
+       
  
        return (sum(ravel(gs[0]) + sum(ravel(gss[0]))), grad)
+ 
      def plot(self, v):
        writegeom(gvline(v[0,:], v[1,:]), 0)