(* 
   compile: ocamlc str.cma hopfield_demo.ml -o hopfield_demo
   execute: ./hopfield 10 testdata1
*)

open Printf

let acoef = 1.0
let bcoef = 1.0
let dcoef = 2.0

let set_pos num pos file = 
  let input = open_in file in
    for i=1 to num do
      let strlist = Str.split_delim (Str.regexp ",") (input_line input) in
      let a = (float_of_string (List.nth strlist 0)) in
      let b = (float_of_string (List.nth strlist 1)) in
	pos.(i-1) <- (a,b)
    done;
    close_in input

let set_distance num pos distance = 
  let distance_total = ref 0.0 in
    for i=1 to num do
      for j=1 to num do
	if i=j then
	  distance.(i-1).(j-1) <- 0.0
	else(
	  let a_x = fst pos.(i-1) in
	  let a_y = snd pos.(i-1) in
	  let b_x = fst pos.(j-1) in
	  let b_y = snd pos.(j-1) in
	  let c = (a_x -. b_x) *. (a_x -. b_x) in
	  let d = (a_y -. b_y) *. (a_y -. b_y) in
	  let e = sqrt (c +. d) in
	    distance.(i-1).(j-1) <- e;
	    distance_total := !distance_total +. distance.(i-1).(j-1));
      done;
    done;
    for i=0 to (num-1) do
      for j=0 to (num-1) do
	distance.(i).(j) <- 10.0 *. distance.(i).(j) /. !distance_total
      done;
    done

let update_state num unit distance =
  for i=0 to (num-1) do
    for j=0 to (num-1) do
      let aterm = ref 0.0 and bterm = ref 0.0 and dterm = ref 0.0 in
      let jm = ref 0 and jp = ref 0 and unitin = ref 0.0 in
	for n=1 to num do
	  aterm := !aterm +. unit.(i).(n-1)
	done;
	aterm := ~-.acoef *. (!aterm -. unit.(i).(j));
	
	for n=1 to num do
	  bterm := !bterm +. unit.(n-1).(j)
	done;
	bterm := ~-.bcoef *. (!bterm -. unit.(i).(j));

	if (j-1) = -1 then jm := num-1 else jm := j-1;
	if (j+1) = num then jp := 0 else jp := j+1;
	
	for n=1 to num do
	  dterm := !dterm +. distance.(i).(n-1) *. (unit.(n-1).(!jp) +. unit.(n-1).(!jm));
	done;
	dterm := ~-.dcoef *. !dterm;
	unitin := !aterm +. !bterm +. !dterm +. acoef +. bcoef;
	unit.(i).(j) <- 0.5 *. (1.0 +. (tanh (!unitin /. 0.5)));
    done;
  done

let energy num unit distance =
  let term1 = ref 0.0 and term2 = ref 0.0 and term3 = ref 0.0 in
  let jm = ref 0 and jp = ref 0 in
    for i=0 to (num-1) do 
      for j=0 to (num-1) do
	for m=0 to (num-1) do
	  if m!=j then term1 := !term1 +. unit.(i).(j) *. unit.(i).(m)
	done;
      done;
    done;
    term1 := 0.5 *. acoef *. !term1;

    for j=0 to (num-1) do
      for i=0 to (num-1) do
	for m=0 to (num-1) do
	  if m!=i then term2 := !term2 +. unit.(i).(j) *. unit.(m).(j)
	done;
      done;
    done;
    term2 := 0.5 *. bcoef *. !term2;

    for i=0 to (num-1) do
      for j=0 to (num-1) do
	if (j-1) = -1 then jm := num-1 else jm := j-1;
	if (j+1) = num then jp := 0 else jp := j+1;
	for m=0 to (num-1) do
	  term3 := !term3 +. distance.(i).(m) *. unit.(i).(j) *. (unit.(m).(!jp) +. unit.(m).(!jm));
	done;
      done;
    done;
    term3 := 0.5 *. dcoef *. !term3;
    
    !term1 +. !term2 +. !term3

let display_unit unit num =
  for i=0 to (num-1) do
    for j=0 to (num-1) do
      printf "%f , " unit.(i).(j);
    done;
    print_newline();
  done

let set_unit num unit =
  for i=0 to (num-1) do
    for j=0 to (num-1) do
      unit.(i).(j) <- Random.float 1.0
    done;
  done


let _ =
  let num = (int_of_string Sys.argv.(1)) in
  let file = Sys.argv.(2) in
  let pos = Array.make num (0.0,0.0) in
  let unit = Array.make_matrix num num 0.0 in
  let distance = Array.make_matrix num num 0.0 in
    Random.self_init ();
    set_pos num pos file;
    set_distance num pos distance;
    set_unit num unit;
    display_unit unit num;
    for i=1 to 100 do
    (*  display_unit unit num; *)
      update_state num unit distance;
      let a = (energy num unit distance) in
	printf "number:%d , Energy=%f\n" i a;
    done;
    display_unit unit num;
    printf "**** end ****\n"