// Examples from chapter 3 "Data Structures" of the book
// F# for Scientists
//
// (C) Flying Frog Consultancy Ltd., 2007

#light

#r @"C:\Program Files\FSharp-1.9.6.2\bin\FSharp.PowerPack.dll"

/// 3
module DataStructures =

  /// 3.1
  module AlgorithmicComplexity =
    let rec ipow_1 x = function
      | 0 -> 1.
      | n -> x * ipow_1 x (n - 1)

    let rec ipow_2 (x : float) = function
      | 0 -> 1.
      | 1 -> x
      | n ->
          let x2 = ipow_2 x (n / 2) in
          x2 * if n % 2 = 0 then x2 else x * x2

    ipow_1 2. 13
    ipow_2 2. 13

  /// 3.2
  module Arrays =
    module ArrayLiterals =
      let a = [|1;2|]
      let b = [|3;4;5|]
      let c = [|6;7;9|]

    open ArrayLiterals

    /// 3.2.2
    module Indexing =
      b.[1]
      c.[2] <- 8
      c
      c.[3] <- 8

    /// 3.2.3
    module Concatenation =
      Array.append a b

      let e = Array.concat [a; b; c]
      e

    /// 3.2.4
    module Aliasing =
      let d = c

      d.[0] <- 17
      c, d

      c.[0] <- 6
      c, d

    let d = Array.copy c

    d.[0] <- 17
    c, d

    /// 3.2.5
    module Subarrays =
      open Concatenation
      Array.sub e 2 5

    /// 3.2.6
    module Creation =
      let a = Array.init 4 (fun n -> n * n)
      a
    
    /// 3.2.7
    module Iteration =
      Array.iter (printf "%d\n") a

    /// 3.2.8
    module Map =
      let b = Array.map (fun n -> n * n) a
      b

    /// 3.2.9
    module Folds =
      Array.fold_left ( + ) 0 b

      let to_list a = Array.fold_right (fun h t -> h :: t) a []
      to_list [|0; 1; 4; 9|]

      Array.to_list [|0; 1; 4; 9|]

    /// 3.2.10
    module Sorting =
      let a = [|1; 5; 3; 4; 7; 9|]
      
      Array.sort compare a
      
      a


    /// 3.2.11
    module PatternMatching =
      let cross (a : float list) b =
        match a, b with
        | [x1; y1; z1], [x2; y2; z2] ->
            [y1 * z2 - z1 * y2; z1 * x2 - x1 * z2; x1 * y2 - y1 * x2]
        | _ -> invalid_arg "cross"

  /// 3.3
  module Lists =
    [1;2] @ [3;4]

    /// 3.3.1
    module Sorting =
      List.sort compare [1;5;3;4;7;9]

    /// 3.3.2
    module Searching =
      /// 3.3.2.1: Membership
      List.mem 4 [1;5;3;4;7;9]

      List.mem 6 [1;5;3;4;7;9]

      /// 3.3.2.2: Predicate
      List.find (fun i -> (i - 6) * i > 0) [1; 3; 4; 5; 7; 9]

      List.find (fun i -> (i - 6) * i = 0) [1; 3; 4; 5; 7; 9]

      [1; 3; 4; 5; 7; 9]
      |> List.find (fun i -> (i - 6) * i > 0)

      /// 3.3.2.3: Association lists
      let list = List.map (fun n -> n, n * n) [1;2;3;4;5]

      list

      List.assoc 4 list

    /// 3.3.3
    module Filtering =
      List.partition (fun x -> x <= 3) [1;2;3;4;5]

      List.filter (fun x -> x <= 3) [1;2;3;4;5]

    /// 3.3.4
    module MapsAndFolds =
      List.map2
      
      let combine a b = List.map2 (fun a b -> a, b) a b

      combine [1; 2; 3] ["a"; "b"; "c"]

      combine [1; 2; 3] [2; 3; 4; 5]

      let vec_add (a : float list) b =
        List.map2 ( + ) a b
      
      vec_add [1.0; 2.0; 3.0] [2.0; 3.0; 4.0]

      let vec_dot (a : float list) b =
        List.fold_left2 (fun d a b -> d + a * b) 0. a b

      vec_dot [1.0; 2.0; 3.0] [2.0; 3.0; 4.0]

    /// 3.3.5
    module PatternMatching =
      let rec down_sample : float list -> float list = function
        | h1 :: h2 :: t -> 0.5 * (h1 + h2) :: down_sample t
        | [] | [_] -> []
      
      down_sample [5.0]
      
      [0.0; 1.0; 0.0; -1.0; 0.0; 1.0; 0.0; -1.0]
      |> down_sample

      module ViaFold =
        let iter f list =
          List.fold_left (fun () e -> f e) () list

        let map f list =
          List.fold_right (fun h t -> f h :: t) list []
      
      let rec iter f = function
        | [] -> ()
        | h :: t ->
            f h
            iter f t

      let rec map f = function
        | [] -> []
        | h :: t -> f h :: map f t
      
      let rec fold_left f accu list =
        match list with
        | h::t -> fold_left f (f accu h) t
        | [] -> accu
      
      let rec fold_right f list accu =
        match list with
        | h::t -> f h (fold_right f t accu)
        | [] -> accu

  /// 3.4
  module Sets =

    Set.empty

    let set1 = Set.singleton 3
    
    set [1; 2; 3]

    /// 3.4.1
    module Insertion =
      let set2 = Set.add 5 set1

      set1, set2

      let s = set [10; 1; 9; 2; 8; 4; 7; 4; 6; 7; 7]

      Set.to_list s

    open Insertion

    /// 3.4.3
    module Cardinality =
      
      Set.cardinal s

    /// 3.4.4
    module SetTheoreticOperations =
      Set.union (set [1;3;5]) (set [3;5;7])

      set [1;3;5] + set [3;5;7]

      Set.inter (set [1;3;5]) (set [3;5;7])

      Set.diff (set [1;3;5]) (set [3;5;7])

      set [1;3;5] - set [3;5;7]
      
      Set.subset (set [2; 4; 6]) s

    set [1;2;3;4;5] = set [5;4;3;2;1]


  /// 3.5
  module Hashtables =
  
    // The OCaml-compatible Hashtbl module has been relegated to the
    // FSharp.PowerPack.dll library.

    let m =
      [ "Hydrogen", 1.0079; "Carbon", 12.011; "Nitrogen", 14.00674;
        "Oxygen", 15.9994; "Sulphur", 32.06 ]
      |> Hashtbl.of_seq

    /// 3.5.2
    module Searching =
      Hashtbl.find m "Carbon"

      m.["Carbon"]

    /// 3.5.3
    module InsertionReplacementAndRemoval =

      Hashtbl.add m "Tantalum" 180.9

      Hashtbl.replace m "Tantalum" 180.948

      m.["Tantalum"] <- 180.9
      
      Hashtbl.remove m "Oxygen"
    
    Hashtbl.iter (printf "%s -> %f\n") m


  /// 3.6
  module Maps =
  
    /// 3.6.1
    module Creation =
      Map.empty
    
      let m =
        Map.of_list
          [ "Hydrogen", 1.0079; "Carbon", 12.011; "Nitrogen", 14.00674;
            "Oxygen", 15.9994; "Sulphur", 32.06 ]

      let list_of m = Map.fold (fun k v t -> (k, v) :: t) m []

    open Creation

    /// 3.6.2
    module Searching =
      Map.find "Carbon" m

      m.["Carbon"]

      let m2 = Map.add "Tantalum" 180.948 m

      list_of m, list_of m2
      
      Map.remove "Oxygen" m

    /// 3.6.3
    module HigherOrderFunctions =

      Map.iter (printf "%s -> %f\n") m


  /// 3.8
  module Sequences =
  
    let double s =
      Seq.map (( * ) 2) s
    
    double [|1; 2; 3|]
    
    double (set [1; 2; 3])


  /// 3.9
  module HeterogeneousContainers =
    type number =
      | Integer of int
      | Real of float
      | Complex of float * float

    let nums = [Integer 1; Real 2.; Complex(3., 4.)]

    let complex_of_number = function
      | Integer n -> complex (float n) 0.
      | Real x -> complex x 0.
      | Complex(x, y) -> complex x y

    List.map complex_of_number nums



  /// 3.10
  module Trees =
    type tree = Node of tree list

    let rec balanced_tree = function
      | 0 -> Node []
      | n -> Node [balanced_tree(n - 1); balanced_tree(n-1)]

    let example = balanced_tree 3

    let rec leaf_count = function
      | Node [] -> 1
      | Node list -> List.fold_left ( + ) 0 (List.map leaf_count list)

    leaf_count example


    type 'a ptree = PNode of 'a * 'a ptree list

    let rec zero_ptree_of_tree (Node list) =
      PNode(0, List.map zero_ptree_of_tree list)

    zero_ptree_of_tree(Node [Node []; Node []])
    
    let aux1 f h (n, t) =
      let n, h = f h n
      n, h::t
    
    let rec aux2 (Node list) n =
      let n', plist =
        List.fold_right (aux1 aux2) list (n+1, [])
      n', PNode(n, plist)
    
    aux2 (Node[Node []; Node []; Node []]) 1

    let counted_ptree_of_tree tree =
      snd(aux2 tree 1)
    
    counted_ptree_of_tree example

    let rec depth_ptree_of_tree (Node list) =
      let list = List.map depth_ptree_of_tree list in
      let depth_of (PNode(d, _)) = d
      let depth = Seq.fold max (-1) (Seq.map depth_of list)
      PNode(depth + 1, list)

    depth_ptree_of_tree example

    let rec unbalanced_tree n =
      let aux n t = Node [balanced_tree n; t] in
      List.fold_right aux [0 .. n - 1] (Node [])

    let wonky_tree = unbalanced_tree 3
    
    wonky_tree

    depth_ptree_of_tree wonky_tree


    type bin_tree = Leaf | Node of bin_tree * bin_tree

    type 'a pbin_tree = Leaf of 'a | Node of 'a pbin_tree * 'a * 'a pbin_tree

    let rec pbin_tree_of_ptree = function
      | PNode(v, []) -> Leaf v
      | PNode(v, [l; r]) ->
          Node(pbin_tree_of_ptree l, v, pbin_tree_of_ptree r)
      | _ -> invalid_arg "pbin_tree_of_ptree"

    pbin_tree_of_ptree (depth_ptree_of_tree (balanced_tree 3))

    // The syntax "('a, 'b) t" for polymorphic arguments is being deprecated
    type ('a, 'b) t =
      | Leaf of 'a
      | Node of ('a, 'b) t * 'b * ('a, 'b) t



  /// 3.10.2
  module UnbalancedTrees =
    type particle = { m: float; r: float }

    /// The abs_float function from OCaml is being deprecated
    let force2 p1 p2 =
      let d = p2.r - p1.r in
      p1.m * p2.m / (d * abs_float d)

    force2 {m=1.; r=0.1} {m=3.; r=0.8}

    /// 3.10.2.1
    module ArrayBasedForceComputation =
      let array_force p ps =
        Array.fold_left (fun f p2 -> f + force2 p p2) 0. ps

      let rand = new System.Random()

      let random_particle _ =
        {m = rand.NextDouble(); r = rand.NextDouble()}
         
      let random_array n = Array.init n random_particle   

      let origin = random_particle()

      let system = random_array 100000

      array_force origin system

    /// 3.10.2.2
    module TreeBasedForceComputation =
      type partition =
        | Leaf of particle list
        | Node of partition * particle * partition

      type accu = {mp: float; mprp: float; p: partition}

      let make_leaf ps =
        let sum = List.fold_left ( + ) 0.
        let m = List.map (fun p -> p.m) ps |> sum
        let mr = List.map (fun p -> p.m * p.r) ps |> sum
        { mp = m; mprp = mr; p = Leaf ps }

      make_leaf [{m=1.; r= -1.}; {m=3.; r=1.}]

      let make_node a1 a2 =
        let mp = a1.mp + a2.mp
        let mprp = a1.mprp + a2.mprp
        { mp = mp;
          mprp = mprp;
          p = Node(a1.p, {m = mp; r = mprp / mp}, a2.p)}

      type system = { lower: float; tree: partition; upper: float }

      // The F# equivalent of OCaml's epsilon_float value has not yet been ossified
      let rec partition x0 ps x2 =
        match ps with
        | [] | [_] as ps -> make_leaf ps
        | ps when x2 - x0 < epsilon_float -> make_leaf ps
        | ps ->
            let x1 = (x0 + x2) / 2.
            let ps1, ps2 = List.partition (fun p -> p.r < x1) ps
            make_node (partition x0 ps1 x1) (partition x1 ps2 x2)

      let rec make_system lower ps upper =
        { lower = lower;
          tree = (partition lower (List.of_seq ps) upper).p;
          upper = upper }

      let particles =
        [| {m=3.; r=0.1};
           {m=1.; r=0.8};
           {m=1.; r=0.82} |]

      make_system 0. particles 1.

      let metric p pp x0 x2 =
        if x0 <= p.r && p.r < x2 then infinity else
          let sqr x = x * x
          let fmin = p.m * pp.m / sqr(p.r - pp.r)
          let g x y = (pp.r - x) / sqr(p.r - y)
          let fmax = p.m * pp.m / (x2 - x0) * (g x0 x2 - g x2 x0)
          fmax - fmin

      metric {m=3.; r=0.1} {m=2.; r=0.81} 0. 1.
      metric {m=3.; r=0.1} {m=2.; r=0.81} 0.5 1.
      metric {m=3.; r=0.1} {m=2.; r=0.81} 0.75 1.
      metric {m=3.; r=0.1} {m=2.; r=0.81} 0.75 0.875

      let rec force_aux delta p x0 x2 = function
        | Leaf list ->
            Seq.fold (fun f p2 -> f + force2 p p2) 0. list
        | Node(left, pp, right) ->
            if metric p pp x0 x2 < delta then
              force2 p pp
            else
              let x1 = 0.5 * (x0 + x2)
              force_aux delta p x0 x1 left +
                force_aux delta p x1 x2 right

      let force p system delta =
        force_aux delta p system.lower system.upper system.tree

      open ArrayBasedForceComputation

      let system' = make_system 0. system 1.

    /// 3.10.2.3
    module PerformanceComparison =
      open ArrayBasedForceComputation
      open TreeBasedForceComputation
    
      #time

      // Let's use a larger number of particles to get more accurate timings
      let system = random_array 10000000
      let system' = make_system 0. system 1.

      force origin system' 1e-9
      // 2.110464073e+15 in 1.349s

      force origin system' 1e-6
      // 2.110464073e+15 in 0.304s

      force origin system' 1e-3
      // 2.110464073e+15 in 0.059s

      array_force origin system
      // 2.110464073e+15 in 0.228s

      #time



  /// 3.10.3
  module AbstractSyntaxTrees =
  
    /// 3.10.3.1
    module Definition =

      type expr =
        | Int of int
        | Var of string
        | Add of expr * expr
        | Mul of expr * expr

      Mul(Add(Int 1, Var "x"), Int 3)

    /// 3.10.3.2
    module EasierConstruction =

      type expr =
        | Int of int
        | Var of string
        | Add of expr * expr
        | Mul of expr * expr
        with
          static member ( + ) (f, g) = Add(f, g)
          static member ( - ) (f, g) = Add(f, Mul(Int(-1), g))
          static member ( * ) (f, g) = Mul(f, g)

      (Int 1 + Var "x") * Int 3

    /// 3.10.3.3
    module EvaluatingExpressions =
      open EasierConstruction

      let rec eval vars = function
        | Int n -> n
        | Var v -> List.assoc v vars
        | Add(f, g) -> eval vars f + eval vars g
        | Mul(f, g) -> eval vars f * eval vars g

      eval ["x", 2] ((Int 1 + Var "x") * Int 3)

    /// 3.10.3.4
    module TermRewriting =
      open EasierConstruction

      let rec rewrite rule expr =
        let expr' =
          match expr with
          | Int _ | Var _ as f -> rule f
          | Add(f, g) -> rule(rewrite rule f + rewrite rule g)
          | Mul(f, g) -> rule(rewrite rule f * rewrite rule g)
        if expr = expr' then expr else rewrite rule expr'

      let expand = function
        | Mul(Add(f, g), h) | Mul(h, Add(f, g)) -> f * h + g * h
        | f -> f

      rewrite expand ((Int 1 + Var "x") * Int 3)