(*
 * lablai - An ML Artificial Inteligence library
 * Copyright (C) 2006  Till Crueger
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *)

(* File $RCSfile$ *)
(* last edited by $Author: till_crueger $ *)
(* $Date: 2008-01-13 03:46:55 +0100 (So, 13 Jan 2008) $, $Revision: 39 $ *)
 
type transfer = Sigmoid | Linear | Zero;;

type neuron = float array * transfer;;

type layer = neuron array;;

type t = layer list;;

let empty : t = [];;

let weight_init = ref 1.0;;

let sigmoid x = 1.0 /. (1.0+.exp(-.x));;
let linear x = x;;
let zero _ = 0.0;;

let apply transfer =
   match transfer with 
      Sigmoid -> sigmoid
   |  Linear -> linear
   |  Zero -> zero
   ;;
   
(* careful, this returns the derivative in terms of the output of the original function *)
let derive transfer =
   match transfer with 
      Sigmoid -> (fun x -> x *. (1.0 -. x))
   |  Linear -> (fun _ -> 1.0)
   |  Zero -> (fun _ -> 0.0)
   ;;


let make_neuron num_syns transfer =
   Array.init (num_syns+1) 
      (fun _ -> (Random.float (!weight_init*.2.0)) -. !weight_init),
   transfer
   ;;

let make_layer num_in num_out transfer = 
   Array.init num_out ( fun _ -> make_neuron num_in transfer)
   ;;
   
let make_net ins outs hidden transfers =
   let size = ins :: hidden @ [outs] in
   let rec loop size transfers =
      if List.length size -1 != List.length transfers then
         raise (Invalid_argument "Layouts do not match")
      else
         match size with
            n :: m :: size -> (make_layer n m (List.hd transfers)) 
                                      :: (loop (m::size) (List.tl transfers))
         |  _ -> []
   in
   loop size transfers
   ;;
   
let make_classifier ins outs hidden =
   let transfers = Sigmoid :: (List.rev_map (fun _ -> Sigmoid) hidden) in
   make_net ins outs hidden transfers
   ;;
   
let make_approximator ins outs hidden =
   let transfers = List.rev_map (fun _ -> Sigmoid) hidden in
   make_net ins outs hidden (transfers @ [Linear])
   ;;
   
let copy net =
   List.map (Array.map (fun (weights,transfer) -> (Array.copy weights,transfer))) net
   ;;
   
let get_synapses (weights,_) =
   weights
   ;;
   
let get_transfer (_,transfer) =
   transfer
   ;;
   
let get_in_size net =
   let neurons = List.hd net in
   Array.length (get_synapses neurons.(0)) -1
   ;;
   
let get_out_size net =
   let neurons = List.hd (List.rev net) in
   Array.length neurons
   ;;
   
let get_depth =
   List.length
   ;;
   
let conc_in_place net1 net2 =
   if get_out_size net1 != get_in_size net2 then
      raise (Invalid_argument "Net size does not match")
   else 
      net1 @ net2
   ;;
   
let conc net1 net2 =
   conc_in_place (copy net1) (copy net2)
   ;;

 let ( *| ) (v1,transfer) v2 =
   let tmp = ref 0.0 in
   for i = 0 to (Array.length v1) -1 do
      tmp := !tmp +. Array.unsafe_get v1 i *. Array.unsafe_get v2 i
   done;
   (apply transfer) !tmp
   ;;
   
let ( |*| ) mat vec =
   let res = Array.make (Array.length mat) 0.0 in
   for i = 0 to Array.length mat -1 do
      Array.unsafe_set res i ((Array.unsafe_get mat i) *| vec)
   done;
   res
   ;;


let evaluate net input =
   (* helper to do the processing of all layers *)
   let rec loop rep activation =
      (* get the current layer *)
      match rep with
         neurons :: rest_rep ->
            begin
               if Array.length (get_synapses neurons.(0)) != Array.length activation +1 then
                  raise (Invalid_argument "Net size does not match input")
               else ();
               (* add the bias *)
               let activation = Array.append [|1.0|] activation in
               (* calculate the weighted sums *)
               let output = neurons |*| activation in
               (* descend to next layer *)
               loop rest_rep output
            end
      |  [] -> activation (* we have reached the output layer  -> return *)
   in
   (* apply inputs to all layers *)
   loop net input
   ;;
   
let train_in_place alpha ?(decay=0.0) () net (input,example) =
   (* unrolled loop for the training *)
   let loop rep activation =
      (* this loop handles layers 2 to n *)
      let rec loop2 rep activation =
         match rep with
            neurons :: rest_rep ->
               begin
                  if Array.length (get_synapses neurons.(0)) != Array.length activation +1 then
                     raise (Invalid_argument "Net size does not match input")
                  else ();
                  (* add the bias *)
                  let activation = Array.append [|1.0|] activation in
                  (* calculate the outputs*)
                  let output = neurons |*| activation in
                  (* descend to next layer and get weighted sums of deltas from there*)
                  let delta_sums = loop2 rest_rep output in
                  (* calculate deltas from those sums *)
                  let deltas = Array.mapi (fun i x -> (derive (get_transfer neurons.(i)) output.(i))*.x) delta_sums in
                  (* create array to hold new deltas*)
                  let new_deltas = Array.make (Array.length activation) 0.0 in
                  
                  (* calculate the differences*)
                  let n = Array.length neurons in
                  let m = Array.length (get_synapses neurons.(0)) in
                  for i = 0 to n - 1 do
                     let weightsi = get_synapses neurons.(i) in
                     let deltai = deltas.(i) in
                     for j = 0 to m - 1 do
                        let weightij = Array.unsafe_get weightsi j in
                        let change = alpha *. deltai *. Array.unsafe_get activation j -. decay*.weightij in
                        (* get this part of the corresponding new_delta *)
                        Array.unsafe_set new_deltas j (Array.unsafe_get new_deltas j +. deltai *. weightij);
                        Array.unsafe_set weightsi j  (weightij +. change);
                     done
                  done;
                  
                  (* pass the deltas to caller, but throw away the one for the bias *)
                  Array.sub new_deltas 1 (Array.length activation - 1) 
               end
         |   [] -> Array.mapi (fun i x -> example.(i) -. x) activation
      in
      (* Special faster case for the first (and most often largest) layer *)
      match rep with 
         neurons :: rest_rep ->
            begin
               if Array.length (get_synapses neurons.(0)) != Array.length activation +1 then
                  raise (Invalid_argument "Net size does not match input")
               else ();
               (* add the bias *)
               let activation = Array.append [|1.0|] activation in
               (* calculate the weighted sums *)
               let output = neurons |*| activation in
               (* descend to next layer and get weighted sums of deltas from there*)
               let delta_sums = loop2 rest_rep output in
               (* calculate deltas from those sums *)
               let deltas = Array.mapi (fun i x -> (derive (get_transfer neurons.(i)) output.(i))*.x) delta_sums in
                
               (* calculate the differences*)
               let n = Array.length neurons in
               let m = Array.length (get_synapses neurons.(0)) in
               for i = 0 to n - 1 do
                  let weightsi = get_synapses neurons.(i) in
                  let deltai = Array.unsafe_get deltas i in
                  for j = 0 to m - 1 do
                     let weightij = Array.unsafe_get weightsi j in
                     let change = alpha *. deltai *. Array.unsafe_get activation j -. decay*.weightij in
                     (* get this part of the corresponding new_delta *)
                     Array.unsafe_set weightsi j  (weightij +. change);
                  done
               done
               (* from the first layer nothing has to be returned *)
            end
      | [] ->
         ()
   in loop net input
   ;;
   
let train alpha ?decay () net (input,example) =
   let new_net = copy net in
   train_in_place alpha ?decay () new_net (input,example) ;
   new_net
   ;;
 
let squared_error net (input,example) =
   let est = evaluate net input in
   let sum = ref 0.0 in
   for i = 0 to Array.length example -1 do
      sum := !sum +. ((example.(i)-.est.(i))**2.0)
   done;
   !sum +. 0.0
   ;;

(* 
 * $Log$
 * Revision 1.8  2008/01/13 02:46:55  till_crueger
 * - Added more types of special teachers to build more complex teaching tasks
 *
 * Revision 1.7  2008/01/12 17:15:07  till_crueger
 *
 * - Added possibilities for a weight decay term for MLPs
 *
 * Revision 1.5  2008/01/12 14:17:21  till_crueger
 * - Changed teacher structure to add a framework for exams after each epoche
 *
 * Revision 1.4  2008/01/11 14:25:50  till_crueger
 * - Changed teacher Interface to be able to combine different teachers
 *
 * Revision 1.3  2007/12/19 01:17:56  till_crueger
 *
 * - Changed Interface for MLPs
 * - Fixed some small bugs for MLPs as approximators
 * - Added q-Learner that uses an MLP as State-Value-Function
 *
 * Revision 1.2  2007/12/17 23:27:29  till_crueger
 *
 * - Added generic teacher for all kinds of supervised learning mechanism
 *
 * Revision 1.1  2007/12/17 22:10:54  till_crueger
 *
 * - Added module for non-object-oriented Multi-Layer-Perceptrons
 *
 *)