# Python implementation of gradient walker

from math import sqrt;

def distance(u,v):
  return sqrt((u[1]-v[1])**2+(u[2]-v[2])**2+(u[0]-v[0])**2);

t0 = sqrt(3)-1;
r0 = 2-sqrt(3);

T0 = [t0,t0,t0];
O = [0,0,0];

def grad(f,X):
  h = global_h;
  gr = []
  for t in range(0,3):
    z = X[:];
    z[t] = z[t] + h;
    gr.append((f(z)-f(X))/h);
  return gr;

def vlen(V):
  return sqrt(sum(map(lambda x:x*x, V)))

def f(X):
  x=X[0]; y=X[1];z=X[2];
  if (x<0 or y<0 or z<0):
    return -1.0;
  a = [1-x, 1-y, 1-z, distance(X,O)-1, distance(X,T0)-r0];
  return min(a)

def refine(func, X, scalA):
  poggle = 10;
  scal = scalA/poggle;
  ssf = func(X); bsf = X;
  for x in range(-poggle,1+poggle):
    for y in range(-poggle,1+poggle):
      for z in range(-poggle,1+poggle):
        U = [X[0]+x*scal, X[1]+y*scal, X[2]+z*scal];
        if (func(U) > ssf):
          bsf = U; ssf = func(U);
  return bsf;

def increfine(func, X, s0):
  v = s0;
  while (v > 1.0e-15):
    X = refine(func, X, v);
    v = v * 0.7;
    print X,f(X),grad(f,X);
  return X;


prec = 10;
scal = 1.0/prec;
global_h = 0.001;
list = [];
coords = [scal*x for x in range(0,1+prec)];
for x in coords:
  for y in coords:
    for z in coords:
      if (f([x,y,z])>0):
        print [x,y,z];
        print "refines to ",increfine(f,[x,y,z],scal)
#        list.append([x,y,z]);

for trial in range(20):
  speed = 0.001 * pow(0.9, trial);
  global_h = speed
  print "Trial", trial," Minimum increment tested ",speed
  list2 = []
  for a in list:
    old = f(a);
    v = grad(f,a);
    len = sqrt(sum(map(lambda x:x*x, v)));
    if(len==0):
# grad=0 implies local maximum
      print "Grad=0 at ",a, f(a);
    else:
      v = map(lambda x:x/len, v);
      bsf = a; ssf = f(bsf);
      scal = 1; w=a;
      while(f(w)>=ssf):
#       print scal,w,f(w)
        bsf=w; ssf=f(w);
        w = [a[i] + scal*speed*v[i] for i in range(0,3)];
#       print scal,w,f(w)
        scal = scal*2;
#    print w, ssf
    list2.append(bsf)
  list = list2[:]
  print "Maximum f-value :",max(map(f,list))
  print "Maximum grad value :",max(map(lambda x:vlen(grad(f,x)), list))

list = map(lambda x : increfine(f,x,0.0001), list)