NB. A J script to draw a basin of attraction
NB. using steepest descent
NB. Cliff Reiter, April 2002
NB. Based upon 
NB. Clifford A. Reiter, Visualizing Steepest Descent, 
NB. The Visual Computer, 8 1 (1991), 64-67. 

NB. needed for file image i/o
require 'fvj2\raster5'

NB. the sample function, gradient and critical points
f4=:[: +/ 0 0 _4 0 1&p.
gf4=:0 _8 0 4&p.
]CP=:,/>{2#<1 _1*%:2

NB. function to compute unit vectors
unit=:% +/&.:*:

NB. step of steepest descent -- includes
NB. a line search technique
sd_step=:2 : 0
y0=.2{.y.
f0=.u. y0
u=.unit v. y0
f1=.u. y0-u
b=.1
if. f1> f0 do.
  while. (f1>f0)*.b>1e_10 do.
  b=.-:b
  f1=. u. y0-b*u
  end.
  else.
  whilst. (f1>f2)*.b<1e10 do.
  b=.+:b
  f1=. u. y0-b*u
  f2=. u. y0-2*b*u
  end.
end.
f2=. u. y0-2*b*u
t=. b*(f2-4*f1)%2*f2-2*f1
ft=. u. y0-t*u
if. ft<f1 do. z=.y0-t*u else. z=.y0-b*u end.
z,y0,>:{:y.
)

NB. halting conditon for iteration
sd_test=:(0 1&{ (1e_6&<)@:dist@:- 2 3&{)*.{: < 253"_
NB. post processing for iteration
sd_post=:CPb@:(0 1&{), 4&{
NB. pre processing for iteration
sd_pre=:], _:, _: , 0:

NB. steepest descent conjunction
SD =: 2 : 0
([:sd_post[:u. sd_step v.^:sd_test^:_ sd_pre)
)

NB. num_pixels zl_clur3 centerleft,upperright
NB. where the right argument is given as complex numbers
NB. utility to find x-y sample points
zl_clur3=: 4 : 0
w=.-~/9 o.y.
h=.-~/11 o.y.
(h-(+:h)*(i.%<:)>.x.*2*h%w) ,~"0/ ({.y.)+w*(i.%<:)x.
)

NB. a useful palette
]p18q=:0,,/0 2|:<.1 0.72 0.45 */ hue 2 4 0 1 5 3{(i.%])6

NB. euclidean distance
dist=:+/&.:*:"1
NB. Critical point basin
CPb=: 3 : 0
1 i.~ 1e_5 > dist CP-"1 y.
)

NB. [num_pix,cen_left,up_rgt] fct mk_sd gradient 'filename.bmp'
NB. make a steepest descent image
mk_sd=:2 : 0
100 _2 2j2 u. mk_sd v. y.
:
z=.({.x.) zl_clur3 }.x.
Z=:i.0 1 2
for_k. i.#z do.
  Z=:Z,(u. SD v.) :: (0 0 253"_)"1 k{z
  if. 0=10|k do. (": k) (1!:2) 2 end.
  end.
'b k'=.0 1|: Z
(p18q;(253~:k)*1+(3|k)+3*b) writebmp8 y.
)

NB. makes a small image in the J directory
50 _8 8j8 f4 mk_sd gf4 'temp_steepest_descent4.bmp'