{- |
Consider a cube of resistors of equal resistance.
What is the overall resistance from one corner to the opposite one?
-}
module Main where

import qualified Numeric.Container as NC
import qualified Data.Packed.Matrix as Matrix
import qualified Data.Packed.Vector as Vector
import qualified Numeric.LinearAlgebra.HMatrix as HMatrix
import Data.Packed.Matrix (Matrix)

import Control.Applicative (liftA3)
import Control.Functor.HT (outerProduct)


{-
Set resistance of a primitive resistor to 1.
This way, voltage equals current.
-}
data Coord = C0 | C1 deriving (Eq, Ord, Show, Enum, Bounded)
data Dim = D0 | D1 | D2 deriving (Eq, Ord, Show, Enum, Bounded)
data Corner = Corner Coord Coord Coord deriving (Eq, Ord, Show)
data Edge = Edge Dim Coord Coord deriving (Eq, Ord, Show)

allCoords :: [Coord]
allCoords = [minBound .. maxBound]

allDims :: [Dim]
allDims = [minBound .. maxBound]

flattenCornerIndex :: Corner -> Int
flattenCornerIndex (Corner x y z) =
   (fromEnum x * 2 + fromEnum y) * 2 + fromEnum z

flattenEdgeIndex :: Edge -> Int
flattenEdgeIndex (Edge d x y) =
   (fromEnum d * 2 + fromEnum x) * 2 + fromEnum y


voltageMatrix :: Matrix Double
voltageMatrix =
   Matrix.fromLists $
   outerProduct
      (\(Edge ed ex ey) c ->
         let ((cx, cy), cz) = selectCornerCoords ed c
         in  if ex==cx && ey==cy
               then
                  case cz of
                     C0 ->  1
                     C1 -> -1
               else 0)
      (liftA3 Edge allDims allCoords allCoords)
      (liftA3 Corner allCoords allCoords allCoords)


-- ToDo: How about cyclic arrangement of dimensions?
selectCornerCoords :: Dim -> Corner -> ((Coord, Coord), Coord)
selectCornerCoords ed (Corner cx cy cz) =
   case ed of
      D0 -> ((cy, cz), cx)
      D1 -> ((cx, cz), cy)
      D2 -> ((cx, cy), cz)

sourceCorner, destCorner :: Corner
sourceCorner = Corner C0 C0 C0
destCorner = Corner C1 C1 C1

currentMatrix :: Matrix Double
currentMatrix =
   Matrix.fromLists $
   outerProduct
      (\c@(Corner _ _ _) (Edge ed ex ey) ->
         let ((cx, cy), cz) = selectCornerCoords ed c
         in  if ex==cx && ey==cy
               then
                  case cz of
                     C0 -> 1
                     C1 -> -1
               else 0)
      (filter (/= sourceCorner) $
       filter (/= destCorner) $
       liftA3 Corner allCoords allCoords allCoords)
      (liftA3 Edge allDims allCoords allCoords)


fullMatrix :: Matrix Double
fullMatrix =
   Matrix.fromBlocks
      [[NC.konst 0 (1,12), Matrix.asRow $ Vector.fromList $ 1 : replicate 7 0],
       [HMatrix.ident 12, voltageMatrix],
       [currentMatrix, NC.konst 0 (6,8)]]

main :: IO ()
main = do
   print fullMatrix
   let [currentVec, potentialVec] =
         Vector.takesV [12,8] $ HMatrix.null1 fullMatrix
   let totalCurrent =
         sum $
         map
            (\d -> NC.atIndex currentVec $ flattenEdgeIndex $ Edge d C0 C0)
            [minBound .. maxBound]
   let cornerPot c = NC.atIndex potentialVec (flattenCornerIndex c)
   let totalVoltage = cornerPot destCorner - cornerPot sourceCorner
   print $ totalVoltage / totalCurrent

{-
result: total resistance is 5/(2+2+2)
-}