packages feed

linear-circuit (empty) → 0.0

raw patch · 7 files changed

+432/−0 lines, 7 filesdep +QuickCheckdep +basedep +comfort-graphsetup-changed

Dependencies added: QuickCheck, base, comfort-graph, containers, hmatrix, linear-circuit, non-empty, transformers, utility-ht

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2015, Henning Thielemann++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Henning Thielemann nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ linear-circuit.cabal view
@@ -0,0 +1,56 @@+Name:                linear-circuit+Version:             0.0+Synopsis:            Compute resistance of linear electrical circuits+Description:+  Compute resistance of linear electrical circuits.+  .+  For examples see test directory.+Homepage:            http://hub.darcs.net/thielema/linear-circuit+License:             BSD3+License-File:        LICENSE+Author:              Henning Thielemann+Maintainer:          haskell@henning-thielemann.de+Category:            Math+Build-Type:          Simple+Cabal-Version:       >=1.10++Source-Repository this+  Tag:         0.0+  Type:        darcs+  Location:    http://hub.darcs.net/thielema/linear-circuit++Source-Repository head+  Type:        darcs+  Location:    http://hub.darcs.net/thielema/linear-circuit++Library+  Exposed-Modules:+    Math.LinearCircuit+  Build-Depends:+    comfort-graph >=0.0 && <0.1,+    hmatrix >=0.16 && <0.17,+    containers >=0.4 && <0.6,+    utility-ht >=0.0.11 && <0.1,+    base >=4.5 && <5+  Hs-Source-Dirs:      src+  Default-Language:    Haskell2010+  GHC-Options:         -Wall++Test-Suite test-linear-circuit+  Type:                exitcode-stdio-1.0+  Hs-Source-Dirs:      test+  Main-is:             Main.hs+  Other-Modules:+    ResistorCube+    Tree+  Build-Depends:+    linear-circuit,+    QuickCheck >=2 && <3,+    comfort-graph,+    non-empty >0.2 && <0.4,+    transformers >=0.4 && <0.5,+    containers,+    utility-ht,+    base+  Default-Language:    Haskell2010+  GHC-Options:         -Wall
+ src/Math/LinearCircuit.hs view
@@ -0,0 +1,96 @@+module Math.LinearCircuit (resistance) where++import qualified Data.Graph.Comfort as Graph+import Data.Graph.Comfort (Graph)++import qualified Numeric.Container as NC+import qualified Numeric.LinearAlgebra.HMatrix as HMatrix+import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import Numeric.LinearAlgebra.HMatrix (Field, (<\>))+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified Data.Map as Map+import qualified Data.List as List+import Data.Monoid (mconcat)++import Control.Functor.HT (outerProduct)+import Data.Bool.HT (if')+++voltageMatrix ::+   (Graph.Edge edge, Ord node, Field a) =>+   Graph edge node a nodeLabel -> Matrix a+voltageMatrix gr =+   Matrix.fromLists $+   outerProduct+      (\e n ->+         if' (Graph.from e == n) 1 $+         if' (Graph.to   e == n) (-1) $+         0)+      (Graph.edges gr)+      (Graph.nodes gr)++{- |+It is almost currentMatrix = trans voltageMatrix,+except that a row is deleted in currentMatrix.+-}+currentMatrix ::+   (Graph.Edge edge, Ord node, Field a) =>+   Graph edge node a nodeLabel -> node -> node -> Matrix a+currentMatrix gr _src dst =+   Matrix.fromLists $+   outerProduct+      (\n e ->+         if' (Graph.from e == n) 1 $+         if' (Graph.to   e == n) (-1) $+         0)+      (List.delete dst $ Graph.nodes gr)+      (Graph.edges gr)++resistanceMatrix ::+   (Graph.Edge edge, Ord node, Field a) =>+   Graph edge node a nodeLabel -> Matrix a+resistanceMatrix gr =+   HMatrix.diag $ Vector.fromList $+   Map.elems $ Graph.edgeLabels gr++fullMatrix ::+   (Graph.Edge edge, Ord node, Field a) =>+   Graph edge node a nodeLabel -> node -> node -> Matrix a+fullMatrix gr src dst =+   let currents = currentMatrix gr src dst+       voltages = voltageMatrix gr+   in  Matrix.fromBlocks+          [[NC.konst 0 (1, Matrix.cols currents),+               Matrix.asRow $ Vector.fromList $+               map (\n -> if n==src then 1 else 0) $ Graph.nodes gr],+           [resistanceMatrix gr, voltages],+           [currents, NC.konst 0 (Matrix.rows currents, Matrix.cols voltages)]]++rhs ::+   (Graph.Edge edge, Ord node, Field a) =>+   Graph edge node a nodeLabel -> node -> node -> Vector a+rhs gr src dst =+   mconcat+      [NC.konst 0 1,+       NC.konst 0 (length (Graph.edges gr)),+       Vector.fromList $+       map (\n -> if n==src then 1 else 0) $+       List.delete dst $ Graph.nodes gr]+++solution ::+   (Graph.Edge edge, Ord node, Field a) =>+   Graph edge node a nodeLabel -> node -> node -> Vector a+solution gr src dst =+   fullMatrix gr src dst <\> rhs gr src dst++resistance ::+   (Graph.Edge edge, Ord node, Field a) =>+   Graph edge node a nodeLabel -> node -> node -> a+resistance gr src dst =+   solution gr src dst+   `NC.atIndex`+   (length (Graph.edges gr) + length (takeWhile (dst/=) $ Graph.nodes gr))
+ test/Main.hs view
@@ -0,0 +1,22 @@+module Main where++import qualified ResistorCube+import qualified Tree++import qualified Test.QuickCheck as QC+++approx :: Double -> Double -> Bool+approx x y = abs (x-y) < 1e-8++test :: (QC.Testable prop) => String -> prop -> IO ()+test msg prop =+   putStr (msg ++ ": ") >> QC.quickCheck prop++main :: IO ()+main = do+   test "resistor cube" (approx ResistorCube.resistance (5/6))+   test "resistor tree"+      (\x -> approx (Tree.treeResistance x) (Tree.graphResistance x))+   test "orientation of resistors"+      (uncurry approx . Tree.flippedResistances)
+ test/ResistorCube.hs view
@@ -0,0 +1,41 @@+{- |+Consider a cube of resistors of equal resistance.+What is the overall resistance from one corner to the opposite one?+-}+module ResistorCube where++import qualified Math.LinearCircuit as LinearCircuit++import qualified Data.Graph.Comfort as Graph+import Data.Graph.Comfort (Graph)++import Control.Applicative (liftA2, liftA3)+++data Coord = C0 | C1 deriving (Eq, Ord, Show, Enum, Bounded)+data Corner = Corner Coord Coord Coord deriving (Eq, Ord, Show)+++allCoords :: [Coord]+allCoords = [minBound .. maxBound]++dimEdges ::+   (Coord -> Coord -> Coord -> Corner) ->+   [(Graph.UndirEdge Corner, Double)]+dimEdges corner =+   liftA2+      (\a b -> (Graph.undirEdge (corner C0 a b) (corner C1 a b), 1))+      allCoords allCoords++graph :: Graph Graph.UndirEdge Corner Double ()+graph =+   Graph.fromList+      (map (flip (,) ()) $ liftA3 Corner allCoords allCoords allCoords)+      (dimEdges (\x y z -> Corner x y z) +++       dimEdges (\y z x -> Corner x y z) +++       dimEdges (\z x y -> Corner x y z))++resistance :: Double+resistance =+   LinearCircuit.resistance graph+      (Corner C0 C0 C0) (Corner C1 C1 C1)
+ test/Tree.hs view
@@ -0,0 +1,184 @@+{- |+Arrange resistors according to a tree of parallel and serial compositions.+Compare resistance of trees with the general graph resistance computation.+-}+module Tree where++import qualified Math.LinearCircuit as LinearCircuit++import qualified Test.QuickCheck as QC++import qualified Data.Graph.Comfort as Graph+import Data.Graph.Comfort (Graph)++import qualified Control.Monad.Trans.Class as MT+import qualified Control.Monad.Trans.State as MS+import Control.Monad (liftM, liftM2, replicateM)++import qualified Data.Map as Map; import Data.Map (Map)++import qualified Data.NonEmpty.Class as NonEmptyC+import qualified Data.NonEmpty as NonEmpty+import qualified Data.Foldable as Fold+import qualified Data.List.Match as Match+import Data.Functor.Classes (Eq1, Ord1, Show1, eq1, compare1, showsPrec1)+import Data.Monoid (mappend)+import Data.Ord.HT (comparing)+import Data.Eq.HT (equating)+++data T =+     Resistance Double+   | Serial (NonEmpty.T [] T)+   | Parallel (NonEmpty.T [] T)+   deriving (Show)++instance QC.Arbitrary T where+   arbitrary =+      let res = liftM Resistance $ QC.choose (0,1)+          go 0 = res+          go size =+            let subTree n =+                  let x = QC.resize (div size n) QC.arbitrary+                  in  liftM2 NonEmpty.cons x (replicateM (n-1) x)+            in  QC.frequency $+                  (3, res) :+                  (1, liftM Serial   (QC.choose (1,size) >>= subTree)) :+                  (1, liftM Parallel (QC.choose (1,size) >>= subTree)) :+                  []+      in  QC.sized go+   shrink tree =+      case tree of+         Resistance res ->+            let simpleRess = [0,1]+            in  if elem res simpleRess+                  then []+                  else map Resistance simpleRess+         Parallel xs -> NonEmpty.flatten xs ++ map Parallel (QC.shrink xs)+         Serial xs   -> NonEmpty.flatten xs ++ map Serial   (QC.shrink xs)+++parallel2 :: Double -> Double -> Double+parallel2 0 0 = 0+parallel2 x y = x*y / (x+y)++treeResistance :: T -> Double+treeResistance x =+   case x of+      Resistance res -> res+      Serial xs -> Fold.foldl1 (+) $ fmap treeResistance xs+      Parallel xs -> Fold.foldl1 parallel2 $ fmap treeResistance xs+++newtype EdgeId = EdgeId Int+   deriving (Eq, Ord, Show)++instance Enum EdgeId where+   fromEnum (EdgeId n) = n+   toEnum = EdgeId++newEdgeId :: (Monad m) => MS.StateT EdgeId m EdgeId+newEdgeId = do+   n <- MS.get+   MS.put $ succ n+   return n++data Edge a =+   Edge {+      edgeId :: EdgeId,+      edgeFrom, edgeTo :: a+   }+   deriving (Show)++instance Eq (Edge a) where (==) = equating edgeId+instance Ord (Edge a) where compare = comparing edgeId++instance Eq1 Edge where eq1 = (==)+instance Ord1 Edge where compare1 = compare+instance Show1 Edge where showsPrec1 = showsPrec++instance Fold.Foldable Edge where+   foldMap f (Edge _ x y) = mappend (f x) (f y)++instance Graph.Edge Edge where+   from (Edge _ n _) = n+   to (Edge _ _ n) = n++instance Graph.Reverse Edge where+   reverseEdge (Edge n from to) = Edge n to from+++newtype Node = Node Int+   deriving (Eq, Ord, Show)++instance Enum Node where+   fromEnum (Node n) = n+   toEnum = Node++newNode :: (Monad m) => MS.StateT Node m Node+newNode = do+   n <- MS.get+   MS.put $ succ n+   return n++edgesFromTree ::+   T -> (Node, Node) ->+   MS.StateT EdgeId (MS.State Node) (Map (Edge Node) Double)+edgesFromTree tree (from, to) =+   case tree of+      Resistance res -> do+         e <- newEdgeId+         return $ Map.singleton (Edge e from to) res+      Serial xs -> do+         ns <- sequence $ Match.replicate (NonEmpty.tail xs) $ MT.lift newNode+         fmap Map.unions $ sequence $+            NonEmpty.flatten $+            NonEmptyC.zipWith edgesFromTree xs $+            NonEmpty.mapAdjacent (,) $+            NonEmpty.cons from $ NonEmpty.snoc ns to+      Parallel xs -> do+         fmap Map.unions $ mapM (flip edgesFromTree (from,to)) $+            NonEmpty.flatten xs++graphFromTree :: T -> (Graph Edge Node Double (), (Node, Node))+graphFromTree tree =+   let ((edgeMap, globalEnds), lastNode) =+         flip MS.runState (Node 0) $ flip MS.evalStateT (EdgeId 0) $ do+            ends <- MT.lift $ liftM2 (,) newNode newNode+            edges <- edgesFromTree tree ends+            return (edges, ends)+   in  (Graph.fromMap+          (Map.fromList $ map (flip (,) ()) [Node 0 .. pred lastNode])+          edgeMap,+        globalEnds)++graphResistance :: T -> Double+graphResistance =+   uncurry (uncurry . LinearCircuit.resistance) . graphFromTree++++data+   FlippedGraph =+      FlippedGraph (Graph Edge Node (Double, Bool) ()) (Node, Node)+   deriving (Show)++instance QC.Arbitrary FlippedGraph where+   arbitrary = do+      (graph, ends) <- fmap graphFromTree QC.arbitrary+      flpGraph <-+         Graph.traverseEdge (\res -> liftM ((,) res) QC.arbitrary) graph+      return $ FlippedGraph flpGraph ends++flippedResistances :: FlippedGraph -> (Double, Double)+flippedResistances (FlippedGraph graph ends) =+   let flippedGraph =+          Graph.fromMap+             (Graph.nodeLabels graph)+             (Map.fromList $+              map+                 (\(e, (res, flp)) ->+                    (if flp then Graph.reverseEdge e else e, res)) $+              Map.toList $ Graph.edgeLabels graph)+   in  (uncurry (LinearCircuit.resistance (Graph.mapEdge fst graph)) ends,+        uncurry (LinearCircuit.resistance flippedGraph) ends)