resistor-cube 0.0.0.4 → 0.0.1
raw patch · 2 files changed
+57/−69 lines, 2 filesdep +comfort-arraydep +lapackdep −hmatrixdep −utility-ht
Dependencies added: comfort-array, lapack
Dependencies removed: hmatrix, utility-ht
Files
- resistor-cube.cabal +4/−4
- src/Main.hs +53/−65
resistor-cube.cabal view
@@ -1,5 +1,5 @@ Name: resistor-cube-Version: 0.0.0.4+Version: 0.0.1 Synopsis: Compute total resistance of a cube of resistors Description: This is an example of how to compute the total resistance@@ -19,7 +19,7 @@ Cabal-Version: >=1.10 Source-Repository this- Tag: 0.0.0.4+ Tag: 0.0.1 Type: darcs Location: http://hub.darcs.net/thielema/resistor-cube @@ -30,8 +30,8 @@ Executable resistor-cube Main-is: Main.hs Build-Depends:- hmatrix >=0.16 && <0.17,- utility-ht >=0.0.11 && <0.1,+ lapack >=0.2.2 && <0.3,+ comfort-array >=0.3.1 && <0.4, base >=4.5 && <5 Hs-Source-Dirs: src Default-Language: Haskell2010
src/Main.hs view
@@ -2,47 +2,51 @@ Consider a cube of resistors of equal resistance. What is the overall resistance from one corner to the opposite one? -}+{-# LANGUAGE TypeOperators #-} 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 qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Format ((##)) -import Control.Applicative (liftA3)-import Control.Functor.HT (outerProduct)+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable ((!))+import Data.Array.Comfort.Shape ((:+:)((:+:))) -{--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)+type Corner = (Coord,Coord,Coord)+type Edge = (Dim,Coord,Coord) -allCoords :: [Coord]-allCoords = [minBound .. maxBound]+type ShapeEnum = Shape.Enumeration+type CornerShape = (ShapeEnum Coord, ShapeEnum Coord, ShapeEnum Coord)+type EdgeShape = (ShapeEnum Dim, ShapeEnum Coord, ShapeEnum Coord) -allDims :: [Dim]-allDims = [minBound .. maxBound]+type Matrix height width = Matrix.General height width -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+cornerShape :: CornerShape+cornerShape = (Shape.Enumeration, Shape.Enumeration, Shape.Enumeration) +edgeShape :: EdgeShape+edgeShape = (Shape.Enumeration, Shape.Enumeration, Shape.Enumeration) -voltageMatrix :: Matrix Double++{-+We also need a currentMatrix that implements Kirchhoff's current law+but it turns out to equal 'transpose voltageMatrix'.+This makes fullMatrix symmetric.+-}+voltageMatrix :: Matrix EdgeShape CornerShape Double voltageMatrix =- Matrix.fromLists $- outerProduct- (\(Edge ed ex ey) c ->+ Matrix.fromRowMajor $+ Array.sample+ (edgeShape, cornerShape)+ (\((ed,ex,ey), c) -> let ((cx, cy), cz) = selectCornerCoords ed c in if ex==cx && ey==cy then@@ -50,60 +54,44 @@ C0 -> 1 C1 -> -1 else 0)- (liftA3 Edge allDims allCoords allCoords)- (liftA3 Corner allCoords allCoords allCoords) selectCornerCoords :: Dim -> Corner -> ((Coord, Coord), Coord)-selectCornerCoords ed (Corner cx cy cz) =+selectCornerCoords ed (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+sourceCorner = (C0,C0,C0)+destCorner = (C1,C1,C1) --- almost transposed voltageMatrix with some removed rows-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)+resistances :: Vector.Vector EdgeShape Double+resistances = Vector.constant edgeShape 1 -fullMatrix :: Matrix Double+fullMatrix :: Triangular.Symmetric (():+:(EdgeShape:+:CornerShape)) 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)]]+ Triangular.stackSymmetric+ (Triangular.symmetricFromList MatrixShape.RowMajor () [0])+ (Matrix.singleRow MatrixShape.RowMajor $+ Vector.unit (edgeShape:+:cornerShape) (Right sourceCorner)) $+ Triangular.stackSymmetric+ (Triangular.diagonal MatrixShape.RowMajor resistances)+ voltageMatrix+ (Vector.constant+ (MatrixShape.symmetric MatrixShape.RowMajor cornerShape) 0) + 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)- allDims- let cornerPot c = NC.atIndex potentialVec (flattenCornerIndex c)- let totalVoltage = cornerPot destCorner - cornerPot sourceCorner- print $ totalVoltage / totalCurrent+ fullMatrix ## "%2.f"+ let solutionVec =+ Matrix.unliftColumn MatrixShape.ColumnMajor+ (Triangular.solve fullMatrix) $+ Vector.unit (():+:(edgeShape:+:cornerShape)) (Right (Right destCorner))+ print $ - solutionVec ! Right (Right destCorner) {- result: total resistance is 5/(2+2+2)