diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,10 @@ +2012-04-25 John D. Ramsdell + + * cmu.cabal (Version): Released as version 1.6. + + * src/Algebra/CommutativeMonoid/LinDiaphEq.hs: Added code that + solves inhomogeneous equations. + 2012-02-29 John D. Ramsdell * cmu.cabal (Version): Released as version 1.5. diff --git a/cmu.cabal b/cmu.cabal --- a/cmu.cabal +++ b/cmu.cabal @@ -1,5 +1,5 @@ Name: cmu -Version: 1.5 +Version: 1.6 Maintainer: ramsdell@mitre.org Cabal-Version: >= 1.2 License: GPL @@ -19,7 +19,7 @@ Library Build-Depends: base >= 3 && < 5, containers >= 0.3, array Exposed-Modules: Algebra.CommutativeMonoid.Unification - Algebra.CommutativeMonoid.HomLinDiaphEq + Algebra.CommutativeMonoid.LinDiaphEq Hs-Source-Dirs: src GHC-Options: -Wall -fno-warn-name-shadowing -fwarn-unused-imports @@ -28,7 +28,7 @@ Main-Is: Algebra/CommutativeMonoid/Main.hs Build-Depends: base >= 3 && < 5, containers >= 0.3, array Other-Modules: Algebra.CommutativeMonoid.Unification - Algebra.CommutativeMonoid.HomLinDiaphEq + Algebra.CommutativeMonoid.LinDiaphEq Hs-Source-Dirs: src GHC-Options: -Wall -fno-warn-name-shadowing -fwarn-unused-imports diff --git a/src/Algebra/CommutativeMonoid/HomLinDiaphEq.hs b/src/Algebra/CommutativeMonoid/HomLinDiaphEq.hs deleted file mode 100644 --- a/src/Algebra/CommutativeMonoid/HomLinDiaphEq.hs +++ /dev/null @@ -1,143 +0,0 @@ --- Homogeneous Linear Diaphantine Equation solver --- --- Copyright (C) 2009 John D. Ramsdell --- --- This program is free software: you can redistribute it and/or modify --- it under the terms of the GNU General Public License as published by --- the Free Software Foundation, either version 3 of the License, or --- (at your option) any later version. - --- This program 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 General Public License for more details. - --- You should have received a copy of the GNU General Public License --- along with this program. If not, see . - --- | --- Module : Algebra.CommutativeMonoid.HomLinDiaphEq --- Copyright : (C) 2009 John D. Ramsdell --- License : GPL --- --- Homogeneous Linear Diaphantine Equation solver. --- --- The solver uses the algorithm of Contejean and Devie as specified --- by David Papp and Bela Vizari in \"Effective Solutions of Linear --- Diophantine Equation Systems with an Application to Chemistry\", --- Rutcor Research Report RRR 28-2004, September, 2004, --- , after --- modification so as to ensure every basis vector is considered. --- --- The algorithm for systems of homogeneous linear Diophantine --- equations follows. Let e[k] be the kth basis vector for 1 <= k <= --- n. To find the minimal, non-negative solutions M to the system of --- equations sum(i=1,n,a[i]*v[i]) = 0, the algorithm of Contejean and --- Devie is: --- --- 1. [init] A := {e[k] | 1 <= k <= n}; M := {} --- --- 2. [new minimal results] M := M + {a in A | a is a solution} --- --- 3. [unnecessary branches] A := {a in A | all m in M : some --- 1 <= k <= n : m[k] < a[k]} --- --- 4. [test] If A = {}, stop --- --- 5. [breadth-first search] A := {a + e[k] | a in A, 1 <= k <= n, --- \ \< 0}; go to step 2 - -module Algebra.CommutativeMonoid.HomLinDiaphEq (homLinDiaphEq) where - -import Data.Array -import Data.Set (Set) -import qualified Data.Set as S - -{-- Debugging hack -import System.IO.Unsafe - -z :: Show a => a -> b -> b -z x y = seq (unsafePerformIO (print x)) y ---} - -type Vector a = Array Int a - -vector :: Int -> [a] -> Vector a -vector n elems = - listArray (0, n - 1) elems - --- | The 'homLinDiaphEq' function takes a list of integers that --- specifies a homogeneous linear Diophantine equation, and returns --- the equation's minimal, non-negative solutions. -homLinDiaphEq :: [Int] -> [[Int]] -homLinDiaphEq [] = [] -homLinDiaphEq v = - newMinimalResults (vector n v) (basis n) S.empty - where n = length v - --- Construct the basis vectors for an n-dimensional space -basis :: Int -> Set (Vector Int) -basis n = - foldl (flip S.insert) S.empty - [ z // [(k, 1)] | - k <- indices z ] - where z = vector n $ replicate n 0 - --- The main loop has been reorganized to ensure every basis vector is --- considered. The breadth-first search step is now the last step. - --- Add elements of a that solve the equation to m and the output -newMinimalResults :: Vector Int -> Set (Vector Int) -> - Set (Vector Int) -> [[Int]] -newMinimalResults v a m = - loop m (S.toList a) -- Test each element in a - where - loop m [] = - nextSearch v a m -- Generate new a and try again - loop m (x:xs) - | prod v x == 0 && S.notMember x m = - elems x:loop (S.insert x m) xs -- Answer found - | otherwise = - loop m xs - --- Generate the next set of test vectors--if there aren't any, your done -nextSearch :: Vector Int -> Set (Vector Int) -> - Set (Vector Int) -> [[Int]] -nextSearch v a m = - if S.null a' then - [] - else - newMinimalResults v (breadthFirstSearch v a') m - where - a' = unnecessaryBranches a m - --- Remove unnecessary branches. A test vector is not necessary if all --- of its elements are greater than or equal to the elements of some --- minimal solution. -unnecessaryBranches :: Set (Vector Int) -> Set (Vector Int) -> Set (Vector Int) -unnecessaryBranches a m = - S.filter f a - where - f x = all (g x) (S.toList m) - g x y = not (lessEq y x) - --- Compare vectors element-wise. -lessEq :: Vector Int -> Vector Int -> Bool -lessEq x y = - all (\i-> x!i <= y!i) (indices x) - --- Breadth-first search using the algorithm of Contejean and Devie -breadthFirstSearch :: Vector Int -> Set (Vector Int) -> Set (Vector Int) -breadthFirstSearch v a = - S.fold f S.empty a - where - f x acc = - foldl (flip S.insert) acc - [ x // [(k, x!k + 1)] | - k <- indices x, - prod v x * v!k < 0 ] -- Contejean-Devie contribution - --- Inner product -prod :: Vector Int -> Vector Int -> Int -prod x y = - sum [ x!i * y!i | i <- indices x ] diff --git a/src/Algebra/CommutativeMonoid/LinDiaphEq.hs b/src/Algebra/CommutativeMonoid/LinDiaphEq.hs new file mode 100644 --- /dev/null +++ b/src/Algebra/CommutativeMonoid/LinDiaphEq.hs @@ -0,0 +1,182 @@ +-- Linear Diaphantine Equation solver +-- +-- Copyright (c) 2009 The MITRE Corporation +-- +-- This program is free software: you can redistribute it and/or modify +-- it under the terms of the GNU General Public License as published by +-- the Free Software Foundation, either version 3 of the License, or +-- (at your option) any later version. + +-- This program 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 General Public License for more details. + +-- You should have received a copy of the GNU General Public License +-- along with this program. If not, see . + +-- | +-- Module : Algebra.CommutativeMonoid.LinDiaphEq +-- Copyright : (C) 2009 John D. Ramsdell +-- License : GPL +-- +-- Linear Diaphantine Equation solver. +-- +-- The solver uses the algorithm of Contejean and Devie as specified +-- in \"An Efficient Incremental Algorithm for Solving Systems of +-- Linear Diophantine Equations\", Information and Computation +-- Vol. 113, pp. 143-174, 1994 after a modification explained below. +-- +-- The algorithm for systems of homogeneous linear Diophantine +-- equations follows. Let e[k] be the kth basis vector for 1 <= k <= +-- n. To find the minimal, non-negative solutions M to the system of +-- equations sum(i=1,n,a[i]*v[i]) = 0, the modified algorithm of +-- Contejean and Devie is: +-- +-- 1. [init] A := {e[k] | 1 <= k <= n}; M := {} +-- +-- 2. [new minimal results] M := M + {a in A | a is a solution} +-- +-- 3. [breadth-first search] A := {a + e[k] | a in A, 1 <= k <= n, +-- \ \< 0} +-- +-- 4. [unnecessary branches] A := {a in A | all m in M : some +-- 1 <= k <= n : m[k] < a[k]} +-- +-- 5. [test] If A = {}, stop, else go to 2. +-- +-- The original algorithm reversed steps 3 and 4. +-- +-- This module provides a solver for a single linear Diophantine +-- equation a*v = b, where a and v are vectors, not matrices. +-- Conceptually, it uses the homogeneous solver after appending -b as +-- the last element of v and by appending 1 to a at each step in the +-- computation. The extra 1 is omitted when an answer is produced. +-- +-- Steps 3 and 4 were switched because the use of the original +-- algorithm for the problem 2x + y - z = 2 produces a non-minimal +-- solution. linDiaphEq [2,1,-1] 2 = [[1,0,0],[0,2,0]], but the +-- original algorithm produces [[1,0,0],[0,2,0],[1,1,1]]. +-- +-- The algorithm is likely to be Fortenbacher's algorithm, the one +-- generalized to systems of equations by Contejean and Devie, but I +-- have not been able to verified this fact. I learned how to extend +-- Contejean and Devie's results to an inhomogeneous equation by +-- reading \"Effective Solutions of Linear Diophantine Equation +-- Systems with an Application to Chemistry\" by David Papp and Bela +-- Vizari, Rutcor Research Report RRR 28-2004, September, 2004, +-- . +-- +-- The example that shows a problem with the original algorithm +-- follows. For the problem linDiaphEq [2,1,-1] 2, the value of a and +-- m at the beginning of the loop is: +-- +-- @ +-- a m +-- [[0, 0, 1], [0, 1, 0], [1, 0, 0]] [] +-- [[0, 1, 1], [0, 2, 0]] [[1, 0, 0]] +-- [] [[1, 0, 0], [0, 2, 0]] +-- @ +-- +-- Consider [0, 1, 1] in a. If you remove unnecessary branches first, +-- the element will stay in a. After performing breadth-first search, +-- a will contain [1, 1, 1], which is the unwanted, non-minimal +-- solution. + +module Algebra.CommutativeMonoid.LinDiaphEq (linDiaphEq) where + +import Data.Array +import Data.Set (Set) +import qualified Data.Set as S + +{-- Debugging hack +import System.IO.Unsafe + +z :: Show a => a -> b -> b +z x y = seq (unsafePerformIO (print x)) y + +zz :: Show a => a -> a +zz x = z x x + +pr :: Set (Vector Int) -> [[Int]] +pr s = map elems $ S.toList s + +zzz :: Set (Vector Int) -> Set (Vector Int) +zzz s = z (pr s) s +--} + +type Vector a = Array Int a + +vector :: Int -> [a] -> Vector a +vector n elems = + listArray (0, n - 1) elems + +-- | The 'linDiaphEq' function takes a list of integers that specifies +-- the coefficients of linear Diophantine equation and a constant, +-- and returns the equation's minimal, non-negative solutions. When +-- solving an inhomogeneous equation, solve the related homogeneous +-- equation and add in those solutions. +linDiaphEq :: [Int] -> Int -> [[Int]] +linDiaphEq [] _ = [] +linDiaphEq v c = + newMinimalResults (vector n v) c (basis n) S.empty + where n = length v + +-- Construct the basis vectors for an n-dimensional space +basis :: Int -> Set (Vector Int) +basis n = + S.fromList [ z // [(k, 1)] | k <- indices z ] + where z = vector n $ replicate n 0 + +-- This is the main loop. + +-- Add elements of a that solve the equation to m and the output +newMinimalResults :: Vector Int -> Int -> Set (Vector Int) -> + Set (Vector Int) -> [[Int]] +newMinimalResults _ _ a _ | S.null a = [] +newMinimalResults v c a m = + loop m (S.toList a) -- Test each element in a + where + loop m [] = -- When done, prepare for next iteration + let a' = breadthFirstSearch v c a -- Step 3 + a'' = unnecessaryBranches a' m in -- Step 4 +-- The original algorithm reverses these two steps. +-- let a' = unnecessaryBranches a m +-- a'' = breadthFirstSearch v c a' in + newMinimalResults v c a'' m + loop m (x:xs) + | prod v x == c && S.notMember x m = + elems x:loop (S.insert x m) xs -- Answer found + | otherwise = + loop m xs + +-- Breadth-first search using the algorithm of Contejean and Devie +breadthFirstSearch :: Vector Int -> Int -> Set (Vector Int) -> Set (Vector Int) +breadthFirstSearch v c a = + S.fold f S.empty a + where + f x acc = + foldl (flip S.insert) acc + [ x // [(k, x!k + 1)] | + k <- indices x, + (prod v x - c) * v!k < 0 ] -- Fortenbacher contribution + +-- Inner product +prod :: Vector Int -> Vector Int -> Int +prod x y = + sum [ x!i * y!i | i <- indices x ] + +-- Remove unnecessary branches. A test vector is not necessary if all +-- of its elements are greater than or equal to the elements of some +-- minimal solution. +unnecessaryBranches :: Set (Vector Int) -> Set (Vector Int) -> Set (Vector Int) +unnecessaryBranches a m = + S.filter f a + where + f x = all (g x) (S.toList m) + g x y = not (lessEq y x) + +-- Compare vectors element-wise. +lessEq :: Vector Int -> Vector Int -> Bool +lessEq x y = + all (\i-> x!i <= y!i) (indices x) diff --git a/src/Algebra/CommutativeMonoid/Unification.hs b/src/Algebra/CommutativeMonoid/Unification.hs --- a/src/Algebra/CommutativeMonoid/Unification.hs +++ b/src/Algebra/CommutativeMonoid/Unification.hs @@ -65,7 +65,7 @@ import Data.List (transpose) import Data.Map (Map) import qualified Data.Map as Map -import Algebra.CommutativeMonoid.HomLinDiaphEq +import Algebra.CommutativeMonoid.LinDiaphEq -- Chapter 8, Section 5 of the Handbook of Automated Reasoning by -- Franz Baader and Wayne Snyder describes unification in @@ -207,7 +207,7 @@ case assocs (add t0 (neg t1)) of [] -> Substitution Map.empty t -> - let basis = homLinDiaphEq (map snd t) in + let basis = linDiaphEq (map snd t) 0 in mgu (map fst t) basis -- Construct a most general unifier the minimal non-negative solutions @@ -265,7 +265,7 @@ -- -- To compute a most general unifier, the set of minimal non-negative -- integer solutions to a linear equation must be found. See module --- Algebra.CommutativeMonoid.HomLinDiaphEq. +-- Algebra.CommutativeMonoid.LinDiaphEq. -- Input and Output