class net ni nh no mu k =
  let sigmoid x = 1. /. (1. +. exp (-.x)) in
  let dsigmoid x = x *. (1. -. x) in 
  (** add bias to both input and hidden layers *)
  let ni         = ni + 1 in 
  let nh         = nh + 1 in
  let rand _ _ = Math.rand (-.1.) 1. in 

object (self)
  val ai      = Array.make ni 1.0;
  val ah      = Array.make nh 0.0;
  val ao      = Array.make no 0.0;
  val wi      = Matrix.init ni nh rand;
  val wo      = Matrix.init nh no rand;  
  val ch      = Array.make nh 0.0;
  val co      = Array.make no 0.0;

  method feedforward input = 

(*    for k = 0 to  *)

    (* copy activations over - skip bias as it is set at 1.0 *)
    for i = 0 to ni - 2 do ai.(i) <- input.(i); done;

    (* hidden activations *)
    for j = 0 to nh - 1 do 
      let sum = ref 0. in 
	for i = 0 to ni - 1 do 
	  sum := !sum +. ai.(i) *. wi.(i).(j);
	done;
	ah.(j) <- sigmoid !sum;
    done;

    (** output activation *)
    for k = 0 to no - 1 do 
      let sum = ref 0. in 
	for j = 0 to nh - 1 do 
	  sum := !sum +. ah.(j) *. wo.(j).(k);
	done;
	ao.(k) <- sigmoid !sum;
	(*ao.(k) <- !sum;*)
    done;
    Array.copy ao;


  method private backprop target = 
    (* calculate error terms for output *)
    for k = 0 to no - 1 do 
      co.(k) <- dsigmoid ao.(k) *. (target.(k) -. ao.(k));
    done;

    (* calc error term for hidden layer *)
    for j = 0 to nh - 1 do 
      let woj = wo.(j) in 
      let sum = ref 0.0 in 
	for k = 0 to no - 1 do 
	  sum := !sum +. co.(k) *. woj.(k);
	done;
	ch.(j) <- dsigmoid ah.(j) *. !sum;
    done;

    (* update hidden weights *)
    for j = 0 to nh - 1 do
      let woj = wo.(j) and ahj = ah.(j) in 
	for k = 0 to no - 1 do 
	  woj.(k) <- woj.(k) +. mu *. co.(k) *. ahj;
	done;
    done;

    (* update input weights *)
    for i = 0 to ni - 1 do 
      let wii = wi.(i) and aii = ai.(i) in 
	for j = 0 to nh - 1 do 
	  wii.(j) <- wii.(j) +. mu *. ch.(j) *. aii;
	done;
    done;

  method train inputs targets epochs =
    for epoch = 0 to epochs - 1 do 
      for i = 0 to Array.length inputs - 1 do 
	let _ = self#feedforward inputs.(i)  in 
	let _ = self#backprop    targets.(i) in ();
      done;
    done;

  method print () = 
    Printf.printf "Input Weights:\n\t";
    (** input layer weights *)
    for i = 0 to ni - 1 do 
      for j = 0 to nh - 1 do 
	Printf.printf "%0.3f " wi.(i).(j);
      done;
      Printf.printf "\n\t";
    done;
    
    Printf.printf "\nHidden Weights:\n\t";    
    (** hidden layer weights *)    
    for j = 0 to nh - 1 do 
      for k = 0 to no - 1 do 
	Printf.printf "%0.3f " wo.(j).(k);
      done;
      Printf.printf "\n\t";
    done;

  method ssq inputs targets = 
    let error = ref 0. in
      for i = 0 to Array.length inputs - 1 do
	let _ = self#feedforward inputs.(i) in 
	  for k = 0 to no - 1 do 
	    error := !error +. 0.5 *. ((targets.(i).(k) -. ao.(k)) ** 2.);
	  done;
      done;
      !error;

  method meanActivity inputs = 
    let outputs = Array.map self#feedforward inputs in 
    let add total out = Array.mapi (fun i _ -> total.(i) +. out.(i)) out in 
    let sum = Array.fold_left add (Array.make (Array.length outputs) 0.) outputs in
    let mean = Array.map (fun v -> v /. float (Array.length inputs)) sum in mean;


  method formalizeBias inputs percent =  
    (** start with hidden layer **)
    for j = 0 to nh - 1 do 
      let a = Array.make (Array.length inputs) 0. in
	for k = 0 to Array.length inputs - 1 do 
	  let input = inputs.(k) in
	  let sum = ref 0. in 
	    (*skip bias *)
	    for i = 0 to ni - 2 do 
	      sum := !sum  +. input.(i) *. wi.(i).(j) 
	    done;
	    (** pre sigmoid! *)
	    a.(k) <- !sum;
	done;
	
	let seq = Math.seq (-.100.) 100. 10000 in 
	  try 
	    for i = 0 to Array.length seq - 1 do 
	      let x = seq.(i) in
	      let s = ref 0. in 
		for j = 0 to Array.length a - 1 do 
		  let y = a.(j) in 
		    s := !s +. (if x > y then 1. else 0.);
		done;
		s := !s /. float (Array.length inputs);
		if !s >= percent then begin
		  wi.(ni-1).(j) <- x;
		  raise Exceptions.Break;
		end;
	    done;
	  with Exceptions.Break -> ();
    done;

    (** output layer *)
    for j = 0 to no - 1 do 
      let a = Array.make (Array.length inputs) 0. in 
	for n = 0 to Array.length inputs - 1 do 
	    let input = inputs.(n) in 

	      (** copy over to include bias *)
	      for i = 0 to ni - 2 do ai.(i) <- input.(i); done;

	      (* skip hidden bias *)
	      for k = 0 to nh - 2 do 
		let sum = ref 0. in 
		  for i = 0 to ni - 1 do sum := !sum +. ai.(i) *. wi.(i).(k); done;
		  ah.(k) <- sigmoid !sum;
	      done;
	      
	      let sum = ref 0. in 
		(** skip bias *)
		for k = 0 to nh - 2 do sum := !sum +. ah.(j) *. wo.(k).(j) done;
		(** pre sigmoid *)
		a.(n) <- !sum;
	  done;
	  
	  let seq = Math.seq (-.100.) 100. 10000 in 
	    try 
	      for i = 0 to Array.length seq - 1 do 
		let x = seq.(i) in
		let s = ref 0. in 
		  for j = 0 to Array.length a - 1 do 
		    let y = a.(j) in 
		      s := !s +. (if x > y then 1. else 0.);
		  done;
		  s := !s /. float (Array.length inputs);
		  if !s >= percent then begin
(*
		    Printf.fprintf stderr "%f %f\n" x !s;
		    flush stderr;
*)
		    wo.(nh-1).(j) <- x;
		    raise Exceptions.Break;
		  end;
	      done;
	    with Exceptions.Break -> ();
      done;

   (* Printf.fprintf stderr "%f\n" (self#meanActivity inputs).(0);*)

  initializer 
  let isum       = Matrix.float_sum wi in 
  let osum       = Matrix.float_sum wo in 
  let norm sum w = w *. (k /. sum) in 
  let _  = Matrix.dmap (norm isum) wi in 
  let _  = Matrix.dmap (norm osum) wo in ();
end;;