diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/apps/ExperimentMain.hs b/apps/ExperimentMain.hs new file mode 100644 --- /dev/null +++ b/apps/ExperimentMain.hs @@ -0,0 +1,34 @@ +module Main where + + +import Options.Applicative + +import qualified Hopfield.Experiments.Experiment +import qualified Hopfield.Experiments.SmallExperiments +import qualified Hopfield.Experiments.ClusterExperiments +import qualified Hopfield.Experiments.Experiment2SuperAttractors + + +data ExperimentArgs = ExperimentArgs { experimentName :: String + , experimentArgs :: [String] + } deriving (Show) + + +argParser :: ParserInfo ExperimentArgs +argParser = info (helper <*> options) ( fullDesc <> header "Runs Hopfield experiments" ) + where + options = ExperimentArgs <$> argument str ( metavar "EXPERIMENT" <> help "the name of the experiment to run" ) + <*> arguments str ( metavar "EXPERIMENT_ARGS" <> help "experiment parameters" ) + + +main :: IO () +main = do + + ExperimentArgs expName args <- execParser argParser + + case expName of + "experiment" -> Hopfield.Experiments.Experiment.main + "small" -> Hopfield.Experiments.SmallExperiments.main + "cluster" -> Hopfield.Experiments.ClusterExperiments.run args + "super" -> Hopfield.Experiments.Experiment2SuperAttractors.main + _ -> error "unknown experiment" diff --git a/apps/Recognize.hs b/apps/Recognize.hs new file mode 100644 --- /dev/null +++ b/apps/Recognize.hs @@ -0,0 +1,181 @@ +{-# LANGUAGE NamedFieldPuns #-} + +module Main where + +import Codec.Picture +import Control.Monad +import Data.Vector ((!)) +import qualified Data.Vector as V +import Options.Applicative +import Text.Printf +import System.Directory + +import Hopfield.Common +import Hopfield.Hopfield +import Hopfield.Images.ConvertImage +import Hopfield.Boltzmann.RestrictedBoltzmannMachine +import Hopfield.Boltzmann.ClassificationBoltzmannMachine +import Hopfield.Benchmark +import Hopfield.Util + + +-- TODO niklas make --fixseed command line option for deterministic results + +-- Function used to transform binary bits (0 or 1) from images to values +-- stored in the network we are using +transformFunction :: Method -> (Int -> Int) +transformFunction Hopfield = (\x -> 2 * x - 1) +transformFunction _ = id + +toPattern :: Method -> CBinaryPattern -> Pattern +toPattern m (CBinaryPattern { cPattern = pat }) = V.fromList $ map (transformFunction m . fromIntegral) $ pat + + +-- | Generates a black-white pixel value from the given pattern. +-- Returns: 'maxBound' if > 0, otherwise 'minBound' for any numeric output type +-- (e.g. 0/255 for Word8). +genPixelBW :: (Bounded a) => Pattern -> Int -> Int -> Int -> a +genPixelBW pattern x y width | pattern ! (y + x * width) > 0 = maxBound + | otherwise = minBound + + +-- | Converts a 'Pattern' to a 8-bit black-white image. +patternToBwImage :: Pattern -> Int -> Int -> Image Pixel8 +patternToBwImage pattern width height = generateImage (genPixelBW pattern width) width height + + +-- | @recPic method (width, height) imgPaths queryImgPath@ recognises a +-- an image given by @queryImgPath@ by using a network of type @method@ which +-- has been trained using @imgPaths@. +-- The images are rescaled accorging to @width@ and @heigth@ before training +-- the network. +recPic :: Method -> (Int, Int) -> [FilePath] -> FilePath -> IO (Either (Image Pixel8) FilePath) +recPic method (width, height) imgPaths queryImgPath = do + l@(_queryImg:_imgs) <- forM (queryImgPath:imgPaths) (\path -> loadPicture path width height) + let queryPat:imgPats = map (toPattern method) l + + result <- case method of + Hopfield -> matchPattern (buildHopfieldData Storkey imgPats) queryPat + Boltzmann -> do d <- buildBoltzmannData imgPats + Right <$> matchPatternBoltzmann d queryPat + CBoltzmann -> do d <- buildCBoltzmannData imgPats + return . Right $ matchPatternCBoltzmann d queryPat + + return $ case result of + -- TODO apply heuristic instead of returning pattern as image (only required for Hopfield) + Left pattern -> Left $ patternToBwImage pattern width height + Right i -> Right $ imgPaths !! i + +-- @saveChain method (width, height) imgPaths queryImgPath@ uses @method@ to train +-- the netwwork using @imgPaths@. Writes to disk all the intermediate images +-- which were produced in the process of mathching @queryImgPath@. +saveChain :: Method -> (Int, Int) -> [FilePath] -> FilePath -> IO () +saveChain method (width, height) imgPaths queryImgPath = do + l@(_queryImg:_imgs) <- forM (queryImgPath:imgPaths) (\path -> loadPicture path width height) + let queryPat:imgPats = map (toPattern method) l + + case method of + Hopfield -> do chain <- updateChain (buildHopfieldData Storkey imgPats) queryPat + mapM_ (putStrLn . patternToAsciiArt width) chain + cleanupDir + mapM_ save $ zip [(0::Int)..] chain + m -> error $ "Method" ++ show m ++ "does not use a chain of images for recognition" + + where + save (number, pattern) = do let filename = printf "converged-images/%.6d.bmp" number + + createDirectoryIfMissing True "converged-images" + writeBitmap filename (patternToBwImage pattern width height) + + cleanupDir = removeDirectoryRecursive "converged-images" + + +data RecognizeArgs = RunOptions + { method :: String + , width :: Int + , height :: Int + , queryPath :: String + , filePaths :: [String] + , saveAllPatterns :: Bool + } + | BenchmarkOptions + { benchmarkPaths :: [String] + } + | InbuiltBenchmarkOptions + { benchmarkName :: String + } + deriving (Show) + + +runOptions :: Parser RecognizeArgs +runOptions = RunOptions <$> argument str ( metavar "METHOD" <> help "hopfield, boltzmann or cboltzmann" ) + <*> argument auto ( metavar "WIDTH" <> help "width images are resized to" ) + <*> argument auto ( metavar "HEIGHT" <> help "height images are resized to" ) + <*> argument str ( metavar "QUERY_PATH" <> help "image to match" ) + <*> arguments str ( metavar "FILE_PATHS" <> help "images to match against (training set)" ) + <*> switch ( long "save-all-patterns" <> help "save all intermediate patterns to harddisk" ) + + +benchmarkOptions :: Parser RecognizeArgs +benchmarkOptions = BenchmarkOptions <$> arguments str ( metavar "FILE_PATHS" <> help "Target for the greeting" ) + + +inbuiltBenchmarkOptions :: Parser RecognizeArgs +inbuiltBenchmarkOptions = InbuiltBenchmarkOptions <$> argument str ( metavar "NAME" <> help "Name of the inbuilt benchmark" ) + + +recognizeOptions :: Parser RecognizeArgs +recognizeOptions = subparser + ( command "run" ( info (helper <*> runOptions) + ( progDesc "Add a file to the repository" )) + <> command "bench" (info (helper <*> benchmarkOptions) + ( progDesc "run benchmark" )) + <> command "inbuiltbench" (info (helper <*> inbuiltBenchmarkOptions) + ( progDesc "run inbuilt benchmark" )) + ) + + +recognizeArgParser :: ParserInfo RecognizeArgs +recognizeArgParser = info (helper <*> recognizeOptions) + ( fullDesc <> header "Performs Hopfield/Boltzmann recognition" + <> progDesc "To see help on individual commands, run --help on them, e.g. recognize run --help." ) + + +main :: IO () +main = do + recArgs <- execParser recognizeArgParser + + case recArgs of + + RunOptions { method, width, height, queryPath, filePaths, saveAllPatterns } + | width < 1 -> error "width must be > 1" + | height < 1 -> error "height must be > 1" + | queryPath == "" -> error "empty query path" + | filePaths == [] -> error "empty query path" + | otherwise -> do + + let recMethod = case method of + "hopfield" -> Hopfield + "boltzmann" -> Boltzmann + "cboltzmann" -> CBoltzmann + _ -> error "unrecognized method" + + if saveAllPatterns + then + saveChain recMethod (width, height) filePaths queryPath + else do + foundPathOrImage <- recPic recMethod (width, height) filePaths queryPath + case foundPathOrImage of + Right path -> putStrLn path + Left image -> do let convergedPath = "converged.bmp" + -- TODO handle return + _ <- writeBitmap convergedPath image + putStrLn $ "no pattern found, wrote coverged image to " ++ convergedPath + + + BenchmarkOptions { benchmarkPaths = _bp } -> error "benchmark not implemented" + + InbuiltBenchmarkOptions { benchmarkName } -> case benchmarkName of + "bench1" -> bench1 + "bench2" -> bench2 + _ -> error "unknown benchmark name" diff --git a/hopfield.cabal b/hopfield.cabal new file mode 100644 --- /dev/null +++ b/hopfield.cabal @@ -0,0 +1,154 @@ +name: hopfield +version: 0.1.0.0 +license: MIT +author: Mihaela Rosca, Lukasz Severyn, Niklas Hambuechen, Razvan Marinescu, Wael Al Jisihi +copyright: Copyright: (c) 2012 Mihaela Rosca, Lukasz Severyn, Niklas Hambuechen, Razvan Marinescu, Wael Al Jisihi +maintainer: Niklas Hambüchen +category: AI, Machine Learning +stability: experimental +synopsis: Hopfield Networks, Boltzmann Machines and Clusters +description: Attractor Neural Networks for Modelling Associative Memory + . + Report: + . + A third year group project at Imperial College London, + supervised by Prof. Abbas Edalat. + . + This projects implements: + . + * Hopfield Networks + . + * Clusters and Super Attractors + . + * The Restricted Boltzmann Machine + . + * A Boltzmann Machine for classification + . + and comes with a range of experiments to evaluate their properties. + +homepage: https://github.com/imperialhopfield/hopfield +bug-Reports: https://github.com/imperialhopfield/hopfield/issues + +build-type: Simple +cabal-version: >= 1.10 + +source-repository head + type: git + location: git://github.com/imperialhopfield/hopfield.git + + +library + default-language: Haskell2010 + exposed-modules: + Hopfield.Hopfield + , Hopfield.Analysis + , Hopfield.Benchmark + , Hopfield.Boltzmann.ClassificationBoltzmannMachine + , Hopfield.Boltzmann.RestrictedBoltzmannMachine + , Hopfield.Clusters + , Hopfield.Common + , Hopfield.Experiments.ClusterExperiments + , Hopfield.Experiments.Experiment + , Hopfield.Experiments.Experiment2SuperAttractors + , Hopfield.Experiments.ExperimentUtil + , Hopfield.Experiments.SmallExperiments + , Hopfield.Images.ConvertImage + , Hopfield.Measurement + , Hopfield.SuperAttractors + , Hopfield.TestUtil + , Hopfield.Util + other-modules: + hs-source-dirs: + src + build-tools: hsc2hs + build-depends: + base >= 4 && <= 5 + , parallel >= 3.1.0.1 + , array >= 0.4.0.0 + , erf >= 2.0.0.0 + , exact-combinatorics >= 0.2.0.4 + , hmatrix >= 0.11.0.4 + , MonadRandom >= 0.1.8 + , probability >= 0.2.4 + , random >= 1.0.1.1 + , random-fu >= 0.2.3.1 + , rvar >= 0.2.0.1 + , vector >= 0.9.1 + , QuickCheck >= 2.4.2 + , deepseq >= 1.3.0.0 + , monad-loops >= 0.3.3.0 + , split >= 0.2.1.1 + c-sources: + src/Hopfield/Images/convertImage.c + include-dirs: + src/Hopfield + , /usr/include/ImageMagick + includes: + wand/magick_wand.h + cc-options: + -- Can't use "-Wextra -Werror" here due to hsc2hs generating unused main() parameters + -g -std=c99 -O0 -Wall -Wextra + ghc-options: + -Wall -fwarn-unused-imports -auto-all + extra-libraries: + MagickWand MagickCore + + +executable experiment + default-language: Haskell2010 + hs-source-dirs: + apps + main-is: + ExperimentMain.hs + other-modules: + build-depends: + base >= 4 && <= 5 + , hopfield + , optparse-applicative >= 0.5.0.0 + ghc-options: + -Wall -fwarn-unused-imports -auto-all -caf-all -rtsopts -threaded + + + +executable recognize + default-language: Haskell2010 + hs-source-dirs: + apps + main-is: + Recognize.hs + other-modules: + build-depends: + base >= 4 && <= 5 + , hopfield + , random >= 1.0.1.1 + , MonadRandom >= 0.1.8 + , vector >= 0.9.1 + , optparse-applicative >= 0.5.0.0 + , JuicyPixels >= 2.0.0 + , directory >= 1.1.0.2 + ghc-options: + + -Wall -fwarn-unused-imports -auto-all -caf-all -rtsopts -threaded + + +Test-Suite tests + default-language: Haskell2010 + type: exitcode-stdio-1.0 + hs-source-dirs: + test + main-is: + Main.hs + build-depends: + base >= 4 + , hopfield + , erf >= 2.0.0.0 + , hspec >= 1.3.0.1 + , HUnit >= 1.2.4.2 + , QuickCheck >= 2.4.2 + , vector >= 0.9.1 + , MonadRandom >= 0.1.8 + , random >= 1.0.1.1 + , exact-combinatorics >= 0.2.0.4 + , parallel >= 3.1.0.1 + ghc-options: + -Wall -fwarn-unused-imports -auto-all -caf-all -rtsopts diff --git a/src/Hopfield/Analysis.hs b/src/Hopfield/Analysis.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Analysis.hs @@ -0,0 +1,76 @@ +module Hopfield.Analysis where + +-- Module for computing the theoretical error of a network +-- Uses the error derivations for independent patterns and super attractors +-- which can be found at the appendix of the report + +import Data.List +import Data.Number.Erf +import qualified Data.Vector as V + + +import Hopfield.Hopfield +import Hopfield.Util + + + +-- | Computes the probability of error for one element given a hopfield data +-- structure. Note that I claim that the actuall error of probability depends +-- on this, but is not the whole term +-- The assumption is that the patterns which were used to train the network +-- are independent. +computeErrorIndependentPats :: HopfieldData -> Double +computeErrorIndependentPats hopfield = computeErrorIndependentPatsNumbers p n + where pats = patterns hopfield + n = V.length $ pats !! 0 + p = length pats + + +-- | computes the error of a super attractor of a hopfield network. The assumption +-- is that there is only one super attractor and the other patterns are independent. +computeErrorSuperAttractor :: HopfieldData -> Double +computeErrorSuperAttractor hopfield = computeErrorSuperAttractorNumbers d n p + where pats = patterns hopfield + n = V.length $ pats !! 0 + p = length pats + d = snd $ maximumBy (compareBy snd) (getElemOccurrences pats) + + +computeErrorIndependentPatsNumbers :: Int -> Int -> Double +computeErrorIndependentPatsNumbers p n + = 1.0 / 2.0 * (1 - (erf $ sqrt $ n ./. (2 * p))) + + +-- | @computeErrorSuperAttractorNumbers d p n@ +-- Computes the probability of error for a super attractor with degree @d@, in +-- a Hopfield network with @n@ neurons, which has been trained with @p@ patterns. +-- The assumption is that the other patterns are independent +-- for mathematical derivation of the equation, see report. +computeErrorSuperAttractorNumbers :: Int -> Int -> Int -> Double +computeErrorSuperAttractorNumbers d p n + = 1.0 / 2.0 * (1.0 - (erf $ (sqrt (n ./. (2 * (p - d)) )))) + + +-- @patternsToNeuronsRatioFromError err@. Given that the err we accept is @err@, +-- returns the maximum ratio between the number of patterns and the number of +-- neurons which can be used to ensure that the probability of error is just @err@. +-- if p/n is grater than @patternsToNeuronsRatioFromError err@ then the error +-- of a Hopfield network will be greater than err. This method is used to compute +-- the minimum number of neurons given the number of training patterns and the +-- maximum error accepted error. +patternsToNeuronsRatioFromError :: Double -> Double +patternsToNeuronsRatioFromError err = 1.0 / (2.0 * (inverf (1.0 - 2.0 * err)) ^ (2 :: Int)) + + + +-- @minNumberOfNeurons p err@ Given the number of patterns used to train a Hopfield +-- network and the maximum error accepted, returns the minimum number of neurons +-- required for the network. +minNumberOfNeurons :: Int -> Double -> Int +minNumberOfNeurons p err + = 1 + floor (p ./ (patternsToNeuronsRatioFromError err)) + + +maxNumberOfPatterns :: Int -> Double -> Int +maxNumberOfPatterns n err + = floor (patternsToNeuronsRatioFromError err *. n) diff --git a/src/Hopfield/Benchmark.hs b/src/Hopfield/Benchmark.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Benchmark.hs @@ -0,0 +1,12 @@ +module Hopfield.Benchmark where + + +import Hopfield.Hopfield +import Hopfield.Clusters +import Hopfield.Experiments.ClusterExperiments + +bench1 :: IO () +bench1 = print =<< experimentUsingT1NoAvg Hebbian 10 10 + +bench2 :: IO () +bench2 = print =<< performAndPrint1 T2 Hebbian 20 5 0.0 0.5 0.5 1 diff --git a/src/Hopfield/Boltzmann/ClassificationBoltzmannMachine.hs b/src/Hopfield/Boltzmann/ClassificationBoltzmannMachine.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Boltzmann/ClassificationBoltzmannMachine.hs @@ -0,0 +1,312 @@ +{-# LANGUAGE PatternGuards, ScopedTypeVariables #-} + +-- | Base Restricted Boltzmann machine. +module Hopfield.Boltzmann.ClassificationBoltzmannMachine where + +-- http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine + +-- Using RBM for recognition +-- http://uai.sis.pitt.edu/papers/11/p463-louradour.pdf +-- http://www.dmi.usherb.ca/~larocheh/publications/drbm-mitacs-poster.pdf + +import Data.Maybe +import Control.Monad +import Control.Monad.Random +import Data.List +import Data.Vector ((!)) +import qualified Data.Vector as V +import qualified Numeric.Container as NC + +import Hopfield.Common +import Hopfield.Util + +-- In the case of the Boltzmann Machine the weight matrix establishes the +-- weights between visible and hidden neurons +-- w i j - connection between visible neuron i and hidden neuron j + +-- | determines the rate in which the weights are changed in the training phase. +-- http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine#Training_algorithm +learningRate :: Double +learningRate = 0.001 + + +data Mode = Hidden | Visible | Classification + deriving(Eq, Show) + + +data BoltzmannData = BoltzmannData { + weightsB :: Weights -- ^ the weights of the network + , classificationWeights :: Weights -- weights for classification + , biasB :: Bias + , biasC :: Bias + , biasD :: Bias + , patternsB :: [Pattern] -- ^ the patterns which were used to train it + -- can be decuded from weights, maybe should be remove now + , hiddenCount :: Int -- ^ number of neurons in the hidden layer + , pattern_to_class :: [(Pattern, Int)] -- the class of the given pattern + -- classes have to be in consecutive order, from 0 +} + deriving(Show) + + +-- | Retrieves the dimension of the weights matrix corresponding to the given mode. +-- For hidden, it is the width of the matrix, and for visible it is the height. +-- One has to ensure that the appropriate weights matrix is passed with this function. +getDimension :: Mode -> Weights -> Int +getDimension Hidden ws = V.length $ ws +getDimension Visible ws = V.length $ ws ! 0 +getDimension Classification ws = V.length $ ws ! 0 + + +-- | @buildCBoltzmannData patterns@ trains a boltzmann network with @patterns@. +-- The number of hidden neurons is set to the number of visible neurons. +buildCBoltzmannData :: MonadRandom m => [Pattern] -> m BoltzmannData +buildCBoltzmannData [] = error "Train patterns are empty" +buildCBoltzmannData pats = + buildCBoltzmannData' pats nr_visible + where nr_visible = V.length (head pats) + + +-- | @buildCBoltzmannData' patterns nrHidden@: Takes a list of patterns and +-- builds a Boltzmann network (by training) in which these patterns are +-- stable states. The result of this function can be used to run a pattern +-- against the network, by using 'matchPatternBoltzmann'. +buildCBoltzmannData' :: MonadRandom m => [Pattern] -> Int -> m BoltzmannData +buildCBoltzmannData' [] _ = error "Train patterns are empty" +buildCBoltzmannData' pats nrHidden + | first_len == 0 + = error "Cannot have empty patterns" + | any (\x -> V.length x /= first_len) pats + = error "All training patterns must have the same length" + | otherwise = trainBoltzmann pats nrHidden + where + first_len = V.length $ head pats + + + +-- | @getActivationProbability ws bias pat index@ +-- can be used to compute the activation probability for a neuron in the +-- visible layer, or for parts of the sums requires for +-- the probability of the classifications +getActivationSum :: Weights -> Bias -> Pattern -> Int -> Double +getActivationSum ws bias pat index + = bias ! index + dotProduct (columnVector ws index) (toDouble pat) + + +-- | @getActivationProbabilityVisible ws bias h index@ returns the activation +-- probability for a neuron @index@ in a visible pattern, given the weights +-- matrix @ws@, the vector of biases @bias@. Applies the activation function +-- to the activation sum, in order to obtain the probability. +getActivationProbabilityVisible :: Weights -> Bias -> Pattern -> Int -> Double +getActivationProbabilityVisible ws bias h index + = activation $ getActivationSum ws bias h index + + +-- | @getActivationSumHidden ws bias h index@ returns the activation +-- sum for a neuron @index@ in a hidden pattern, given the weights +-- matrix @ws@, the vector of biases @bias@. +getActivationSumHidden :: Weights -> Weights -> Bias -> Pattern -> Pattern -> Int -> Double +getActivationSumHidden ws u c v y index + | Just e <- validPattern Visible ws v = error e + | Just e <- validPattern Classification u y = error e + | otherwise = c ! index + dotProduct (ws ! index) (toDouble v) + dotProduct (u ! index) (toDouble y) + +-- | @getActivationSumHidden ws bias h index@ returns the activation +-- sum for all neurons in the hidden pattern, given the weights +-- matrix @ws@, the vector of biases @bias@. +getHiddenSums :: Weights -> Weights -> Bias -> Pattern -> Pattern -> V.Vector Double +getHiddenSums ws u c v y + = V.fromList [getActivationSumHidden ws u c v y i | i <- [0 .. (V.length ws) - 1] ] + + +-- | @getActivationProbabilityVisible ws u bias v index@ returns the activation +-- probability for a neuron @index@ in a hidden pattern, given the weights +-- matrices @ws@ and @u@, the vector of biases @bias@. Applies the activation function +-- to the activation sum, in order to obtain the probability. +getActivationProbabilityHidden :: Weights -> Weights -> Bias -> Pattern -> Pattern -> Int -> Double +getActivationProbabilityHidden ws u c v y index + = activation $ getActivationSumHidden ws u c v y index + + +-- | @updateNeuronVisible ws bias h index@ updates a neuron in the visible layer by using gibbsSampling, according +-- to the activation probability +updateNeuronVisible :: MonadRandom m => Weights -> Bias -> Pattern -> Int -> m Int +updateNeuronVisible ws bias h index + = gibbsSampling $ getActivationProbabilityVisible ws bias h index + + +-- | Updates a neuron in the hidden layer by using gibbsSampling, according +-- to the activation probability +updateNeuronHidden :: MonadRandom m => Weights -> Weights -> Bias -> Pattern -> Pattern -> Int -> m Int +updateNeuronHidden ws u c v y index + = gibbsSampling $ getActivationProbabilityHidden ws u c v y index + + +-- | Updates the entire visible layer by using gibbsSampling, according +-- to the activation probability +updateVisible :: MonadRandom m => Weights -> Bias -> Pattern -> m Pattern +updateVisible ws bias h + | Just e <- validPattern Hidden ws h = error e + | otherwise = V.fromList `liftM` mapM (updateNeuronVisible ws bias h) updatedIndices + where + updatedIndices = [0 .. (V.length $ ws ! 0) - 1] + + +-- | Updates the entire visible layer by using gibbsSampling, according +-- to the activation probability +updateHidden :: MonadRandom m => Weights -> Weights -> Bias -> Pattern -> Pattern -> m Pattern +updateHidden ws u c v y + | Just e <- validPattern Visible ws v = error e + | otherwise = V.fromList `liftM` mapM (updateNeuronHidden ws u c v y) updatedIndices + where + updatedIndices = [0 .. (V.length ws) - 1 ] + + +-- | Updates a classification vector given the current state of the network ( +-- the u matrix and the vector of biases d, together with a hidden vector h) +updateClassification :: Weights -> Bias -> Pattern -> Pattern +updateClassification u d h + = V.fromList [ if n == newClass then 1 else 0 | n <- allClasses] + where + -- TODO replace with actual sampling using inverse method (with cdf list) + expActivation = exp . (getActivationSum u d h) + newClass = maximumBy (compareBy expActivation) allClasses + allClasses = [0 .. nrClasses - 1] + nrClasses = V.length d + + +-- @getClassificationVector pat_to_classes pat@ returns the classification +-- vector of @pat@, by looking up in @pat@ in @pat_to_classes@ to obtain the +-- class of the pattern. The classification vector is obtained by +-- creating vector with all 0s and only 1 in the position of the class. +-- The length of all classification vectors is the number of classes. +getClassificationVector :: [(Pattern, Int)] -> Pattern -> Pattern +getClassificationVector pat_classes pat + = V.fromList [ if n == pat_class then 1 else 0 | n <- map snd pat_classes] + where pat_class = fromJust $ lookup pat pat_classes + + + +-- | One step which updates the weights in the CD-n training process. +-- The weights are changed according to one of the training patterns. +-- http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine#Training_algorithm +-- @oneTrainingStep bm visible@ updates the parameters of @bm@ (the 2 weight +-- matrices and the biases) according to the training instance @v@ +-- and its classification, obtained by looking in the map kept in @bm@ +oneTrainingStep :: MonadRandom m => BoltzmannData -> Pattern -> m BoltzmannData +oneTrainingStep (BoltzmannData ws u b c d pats nr_h pat_to_class) v = do + let y = getClassificationVector pat_to_class v + h_sum = getHiddenSums ws u c v y + h <- updateHidden ws u c v y + v' <- updateVisible ws b h + let y' = updateClassification u d h + (h_sum' :: V.Vector Double) = getHiddenSums ws u c v' y' + getOuterProduct v1 v2 = NC.toLists $ (fromDataVector v1) `NC.outer` (fromDataVector $ toDouble v2) + getDelta pos neg = map (map (* learningRate)) $ combine (-) pos neg + updateWeights w d_w = vector2D $ combine (+) (list2D w) d_w + deltaBias v1 v2 = V.map ((* learningRate) . fromIntegral) (combineVectors (-) v1 v2) + deltaBiasC v1 v2 = V.map (* learningRate) (combineVectors (-) v1 v2) + updateBias bias delta_bias = combineVectors (+) bias delta_bias + pos_ws = getOuterProduct h_sum v -- "positive gradient for ws" + neg_ws = getOuterProduct h_sum' v' -- "negative gradient for ws" + pos_u = getOuterProduct h_sum y -- "positive gradient for u" + neg_u = getOuterProduct h_sum' y' -- "negative gradient for u" + d_ws = getDelta pos_ws neg_ws -- "delta ws" + new_ws = updateWeights ws d_ws + d_u = getDelta pos_u neg_u -- "delta u" + new_u = updateWeights u d_u + new_b = updateBias b (deltaBias v v') + new_c = updateBias c (deltaBiasC h_sum h_sum') + new_d = updateBias d (deltaBias y y') + return $ BoltzmannData new_ws new_u new_b new_c new_d pats nr_h pat_to_class + + +-- | The training function for the Boltzmann Machine. +-- We are using the contrastive divergence algorithm CD-1 +-- TODO see if making the vis +-- (we could extend to CD-n, but "In pratice, CD-1 has been shown to work surprisingly well." +-- @trainBoltzmann pats nrHidden@ where @pats@ are the training patterns +-- and @nrHidden@ is the number of neurons to be created in the hidden layer. +-- http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine#Training_algorithm +trainBoltzmann :: MonadRandom m => [Pattern] -> Int -> m BoltzmannData +trainBoltzmann pats nr_h = do + ws <- vector2D `liftM` genWeights + u <- vector2D `liftM` genU + foldM oneTrainingStep (BoltzmannData ws u b c d pats nr_h pats_classes) pats + where + genWeights = replicateM nr_h . replicateM nr_visible $ normal 0.0 0.01 + genU = replicateM nr_h . replicateM nr_classes $ normal 0.0 0.01 + b = V.fromList $ replicate nr_visible 0.0 + c = V.fromList $ replicate nr_h 0.0 + d = V.fromList $ replicate nr_classes 0.0 + nr_classes = length nub_pats + nub_pats = nub pats + pats_classes = zip nub_pats [0 .. ] + nr_visible = V.length $ head pats + + +-- | @matchPatternBoltzmann bm pat@ given the Boltzmann trained network @bm@ +-- recognizes @pat@, by classifying it to one of the patterns the network was +-- trained with. This is done by computing the free energy of @pat@ with +-- every possible classification, and choosing the classification with +-- lowest energy. +-- http://uai.sis.pitt.edu/papers/11/p463-louradour.pdf +matchPatternCBoltzmann :: BoltzmannData -> Pattern -> Int +matchPatternCBoltzmann bm v + | Just e <- validPattern Visible (weightsB bm) v = error e + | otherwise = fromJust $ maxPat `elemIndex` pats + where + pats_classes = pattern_to_class bm + pats = patternsB bm + patternsWithClassifications = [ (p, getClassificationVector pats_classes p) | p <- map fst pats_classes] + probability classification = exp $ - (getFreeEnergy bm v classification) + (maxPat, _) = maximumBy (compareBy $ probability . snd) patternsWithClassifications + + +-- | @getFreeEnergy bm visible classification_vector@ +-- Computes the free energy of @v@ with @classification_vector@, according +-- to the trained Boltzmann network @bm@. It is used for classifying a given +-- visible vector according to the classes used for training the network @bm@. +getFreeEnergy :: BoltzmannData -> Pattern -> Pattern -> Double +getFreeEnergy (BoltzmannData ws u _b c d _pats _nrH _pats_classes) v y + = - dotProduct d (toDouble y) - (V.sum $ V.map softplus hiddenSums) + where hiddenSums = getHiddenSums ws u c v y + + +-- | The activation function for the network (the logistic sigmoid). +-- http://en.wikipedia.org/wiki/Sigmoid_function +activation :: Double -> Double +activation x = 1.0 / (1.0 + exp (-x)) + +-- | The function used to compute the free energy +-- http://uai.sis.pitt.edu/papers/11/p463-louradour.pdf +softplus :: Double -> Double +softplus x = log (1.0 + exp x) + + +-- TODO move to tests +validClassificationVector :: Pattern -> Int -> Maybe String +validClassificationVector pat size + | V.length pat /= size = Just "classification vector does not match expected size" + | V.any (\x -> notElem x [0, 1]) pat = Just "Non binary element in classification pattern" + | V.sum pat /=1 = Just "Invalid classification vector" + | otherwise = Nothing + + +-- | @validPattern mode weights pattern@ +-- Returns an error string in a Just if the @pattern@ is not compatible +-- with @weights@ and Nothing otherwise. @mode@ gives the type of the pattern, +-- which is checked (Visible or Hidden). +validPattern :: Mode -> Weights -> Pattern -> Maybe String +validPattern mode ws pat + | getDimension mode ws /= V.length pat = Just $ "Size of pattern must match network size in " ++ show mode + | V.any (\x -> notElem x [0, 1]) pat = Just "Non binary element in Boltzmann pattern" + | otherwise = Nothing + +-- | @validWeights ws@ checks that a weight matrix is well formed. +validWeights :: Weights -> Maybe String +validWeights ws + | V.null ws = Just "The matrix of weights is empty" + | V.any (\x -> V.length x /= V.length (ws ! 0)) ws = Just "Weights matrix ill formed" + | otherwise = Nothing + diff --git a/src/Hopfield/Boltzmann/RestrictedBoltzmannMachine.hs b/src/Hopfield/Boltzmann/RestrictedBoltzmannMachine.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Boltzmann/RestrictedBoltzmannMachine.hs @@ -0,0 +1,228 @@ +{-# LANGUAGE PatternGuards, ScopedTypeVariables #-} + +-- | Base Restricted Boltzmann machine. +-- http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine +module Hopfield.Boltzmann.RestrictedBoltzmannMachine where + + +import Data.Maybe +import Control.Monad +import Control.Monad.Random +import Data.List +import Data.Vector ((!)) +import qualified Data.Vector as V +import qualified Numeric.Container as NC + +import Hopfield.Common +import Hopfield.Util + + +-- In the case of the Boltzmann Machine the weight matrix establishes the +-- weights between visible and hidden neurons +-- w i j - connection between visible neuron i and hidden neuron j + +-- | determines the rate in which the weights are changed in the training phase. +-- http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine#Training_algorithm +learningRate :: Double +learningRate = 0.1 + + +data Mode = Hidden | Visible + deriving(Eq, Show) + + +data Phase = Training | Matching + deriving(Eq, Show) + + +data BoltzmannData = BoltzmannData { + weightsB :: Weights -- ^ the weights of the network + , patternsB :: [Pattern] -- ^ the patterns which were used to train it + , nr_hiddenB :: Int -- ^ number of neurons in the hidden layer + , pattern_to_binaryB :: [(Pattern, [Int])] -- ^ the binary representation of the pattern index + -- the pattern_to_binary field will not replace the patternsB field as it does + -- not contain duplicated patterns, which might be required for statistical + -- analysis in clustering and super attractors +} + deriving(Show) + + +-- | Retrieves the dimension of the weights matrix corresponding to the given mode. +-- For hidden, it is the width of the matrix, and for visible it is the height. +getDimension :: Mode -> Weights -> Int +getDimension Hidden ws = V.length $ ws ! 0 +getDimension Visible ws = V.length $ ws + + +notMode :: Mode -> Mode +notMode Visible = Hidden +notMode Hidden = Visible + + +-- | @buildBoltzmannData patterns@ trains a boltzmann network with @patterns@. +-- The number of hidden neurons is set to the number of visible neurons. +buildBoltzmannData :: MonadRandom m => [Pattern] -> m BoltzmannData +buildBoltzmannData [] = error "Train patterns are empty" +buildBoltzmannData pats = + buildBoltzmannData' pats nr_visible + where nr_visible = fromIntegral $ V.length (head pats) + + +-- | @buildBoltzmannData' patterns nr_hidden@: Takes a list of patterns and +-- builds a Boltzmann network (by training) in which these patterns are +-- stable states. The result of this function can be used to run a pattern +-- against the network, by using 'matchPatternBoltzmann'. +buildBoltzmannData' :: MonadRandom m => [Pattern] -> Int -> m BoltzmannData +buildBoltzmannData' [] _ = error "Train patterns are empty" +buildBoltzmannData' pats nr_hidden + | first_len == 0 + = error "Cannot have empty patterns" + | any (\x -> V.length x /= first_len) pats + = error "All training patterns must have the same length" + | otherwise = do + (ws, pats_with_binary) :: (Weights, [(Pattern, [Int])]) <- trainBoltzmann pats nr_hidden + return $ BoltzmannData ws pats nr_hidden pats_with_binary + where + first_len = V.length (head pats) + + +-- Pure version of updateNeuron for testing +updateNeuron' :: Double -> Phase -> Mode -> Weights -> Pattern -> Int -> Int +updateNeuron' r phase mode ws pat index = if (r < a) then 1 else 0 + where a = getActivationProbability phase mode ws pat index + +-- +getActivationProbability :: Phase -> Mode -> Weights -> Pattern -> Int -> Double +getActivationProbability phase mode ws pat index = if a <=1 && a >=0 then a else error (show a) + where + a = activation . sum $ case mode of + Hidden -> [ (ws ! index ! i) *. (pat' ! i) | i <- [0 .. p-1] ] + Visible -> [ (ws ! i ! index) *. (pat' ! i) | i <- [0 .. p-1] ] + pat' = case phase of + Matching -> V.cons 1 pat + Training -> pat + p = V.length pat' + + +-- | @updateNeuron mode ws pat index@ , given a vector @pat@ of type @mode@ +-- updates the neuron with number @index@ in the layer with opposite type. +updateNeuron :: MonadRandom m => Phase -> Mode -> Weights -> Pattern -> Int -> m Int +updateNeuron phase mode ws pat index = do + r <- getRandomR (0.0, 1.0) + return $ updateNeuron' r phase mode ws pat index + + +-- | @getCounterPattern mode ws pat@, given a vector @pat@ of type @mode@ +-- computes the values of all the neurons in the layer of the opposite type. +getCounterPattern :: MonadRandom m => Phase -> Mode -> Weights -> Pattern -> m Pattern +getCounterPattern phase mode ws pat + | Just e <- validPattern phase mode ws pat = error e + | otherwise = V.fromList `liftM` mapM (updateNeuron phase mode ws pat) updatedIndices + where + updatedIndices = [0 .. getDimension (notMode mode) ws - diff] + diff = case phase of + Training -> 1 + Matching -> 2 + + +-- | One step which updates the weights in the CD-n training process. +-- The weights are changed according to one of the training patterns. +-- http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine#Training_algorithm +updateWeights :: MonadRandom m => Weights -> Pattern -> m Weights +updateWeights ws v = do + let biased_v = V.cons 1 v + h <- getCounterPattern Training Visible ws biased_v + v' <- getCounterPattern Training Hidden ws h + let f = fromDataVector . fmap fromIntegral + pos = NC.toLists $ (f biased_v) `NC.outer` (fromDataVector $ getSigmaH v) -- "positive gradient" + neg = NC.toLists $ (f v') `NC.outer` (fromDataVector $ getSigmaH v') -- "negative gradient" + d_ws = map (map (* learningRate)) $ combine (-) pos neg -- weights delta + new_weights = combine (+) (list2D ws) d_ws + nr_hidden = V.length $ ws ! 0 + getSigmaH y = V.fromList [getActivationProbability Training Visible ws y x | x <- [0.. nr_hidden - 1] ] + return $ vector2D new_weights + + +-- | The training function for the Boltzmann Machine. +-- We are using the contrastive divergence algorithm CD-1 +-- TODO see if making the vis +-- (we could extend to CD-n, but "In practice, CD-1 has been shown to work surprisingly well." +-- @trainBoltzmann pats nr_hidden@ where @pats@ are the training patterns +-- and @nr_hidden@ is the number of neurons to be created in the hidden layer. +-- http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine#Training_algorithm +trainBoltzmann :: MonadRandom m => [Pattern] -> Int -> m (Weights, [(Pattern, [Int])]) +trainBoltzmann pats nr_hidden = do + weights_without_bias <- genWeights + -- add biases as a dimension of the matrix, in order to include them in the + -- contrastive divergence algorithm + let ws = [0: x | x <- weights_without_bias] + ws_start = (replicate (nr_hidden + 1) 0) : ws + updated_ws <- foldM updateWeights (vector2D ws_start) pats' + return (updated_ws, paths_with_binary_indices) + where + genWeights = replicateM nr_visible . replicateM nr_hidden $ normal 0.0 0.01 + paths_with_binary_indices = getBinaryIndices pats + pats' = [(V.++) x $ encoding x | x <- pats] + encoding x = V.fromList . fromJust $ lookup x paths_with_binary_indices + nr_visible = V.length $ pats' !! 0 + + +-- | The activation function for the network (the logistic sigmoid). +-- http://en.wikipedia.org/wiki/Sigmoid_function +activation :: Double -> Double +activation x = 1.0 / (1.0 + exp (-x)) + + +-- | @validPattern mode weights pattern@ +-- Returns an error string in a Just if the @pattern@ is not compatible +-- with @weights@ and Nothing otherwise. @mode@ gives the type of the pattern, +-- which is checked (Visible or Hidden). +validPattern :: Phase -> Mode -> Weights -> Pattern -> Maybe String +validPattern phase mode ws pat + | checked_dim /= V.length pat = Just $ "Size of pattern must match network size in " ++ show phase ++ " " ++ show mode + | V.any (\x -> notElem x [0, 1]) pat = Just "Non binary element in Boltzmann pattern" + | otherwise = Nothing + where checked_dim = if phase == Training then actual_dim else actual_dim - 1 + actual_dim = getDimension mode ws + + +validWeights :: Weights -> Maybe String +validWeights ws + | V.null ws = Just "The matrix of weights is empty" + | V.any (\x -> V.length x /= V.length (ws ! 0)) ws = Just "weights matrix ill formed" + | otherwise = Nothing + + +-- | Updates a pattern using the Boltzmann machine +updateBoltzmann :: MonadRandom m => Weights -> Pattern -> m Pattern +updateBoltzmann ws pat = do + h <- getCounterPattern Matching Visible ws pat + getCounterPattern Matching Hidden ws h + + +-- see http://www.cs.toronto.edu/~hinton/absps/guideTR.pdf section 16.1 +-- And stack overflow discussion +-- http://stackoverflow.com/questions/9944568/the-free-energy-approximation-equation-in-restriction-boltzmann-machines +-- http://www.dmi.usherb.ca/~larocheh/publications/class_set_rbms_uai.pdf +getFreeEnergy :: Weights -> Pattern -> Double +getFreeEnergy ws pat + | Just e <- validWeights ws = error e + | Just e <- validPattern Matching Visible ws pat = error e + | otherwise = - biases - sum (map f xs) + where w i j = ((ws :: Weights) ! i ! j) :: Double + biases = sum [ w (i + 1) 0 *. (pat ! i) | i <- [0 .. p - 1] ] + xs = [ w 0 j + sum [ w (i + 1) j *. (pat ! i) | i <- [0 .. p - 1] ] | j <- [1 .. (V.length $ ws ! 0) - 1]] + f x = log (1 + exp x) + p = V.length pat + + +-- | Matches a pattern against the a given network +matchPatternBoltzmann :: MonadRandom m => BoltzmannData -> Pattern -> m Int +matchPatternBoltzmann (BoltzmannData ws pats _ pats_with_binary) pat = do + hot_pat <- updateBoltzmann ws ((V.++) pat (V.fromList $ snd $ head pats_with_binary)) + let h = V.take (V.length $ head pats) hot_pat + extendWithClass p = ((V.++) h (V.fromList . fromJust $ lookup p pats_with_binary) ) + getPatternProbability x = exp $ (- getFreeEnergy ws x) + fromPatToIndex p = fromJust $ p `elemIndex` pats + return $ fst $ maximumBy (compareBy snd) [(fromPatToIndex p, (getPatternProbability . extendWithClass) p) | p <- pats] + diff --git a/src/Hopfield/Clusters.hs b/src/Hopfield/Clusters.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Clusters.hs @@ -0,0 +1,157 @@ +{-# LANGUAGE PatternGuards #-} + +module Hopfield.Clusters where + + +-- Module which deals with pattern cluster generation and related functions. +-- Implements probabilistic rewiring using Hamming distance. + +import qualified Data.Vector as V +import Control.Monad.Random (MonadRandom) +import Control.Monad (liftM, replicateM) + +import Hopfield.Common +import Hopfield.Hopfield +import Hopfield.Measurement +import Hopfield.Util + + +-- |@getPatternInCluster pat p@ gets a pattern in a cluster given by @pat@ +-- by flipping each bit in the pattern with probability p. +getPatternInCluster :: MonadRandom m => Method -> Pattern -> Double -> m Pattern +getPatternInCluster method originPat p + = liftM V.fromList $ mapM transformBit (V.toList originPat) + where transformBit x = do + flip_bit <- gibbsSampling p + let bit = if (odd flip_bit) then (flipBit method x) else x + return bit + + +-- |@getPatternInCluster pat p@ gets a pattern in a cluster given by @pat@ +-- by flipping each bit in the pattern with probability p. +getCluster :: MonadRandom m => Method -> Pattern -> Int -> Double -> m [Pattern] +getCluster method originPat size p + = replicateM size (getPatternInCluster method originPat p) + + +-- Caller has to take care with setting the mean and stdDev such that +-- the sampled numbers tend to be in the interval [0 .. size -1] +-- Implements the T2 method described by Federico +-- Sample a Gaussian distribution with given mean and std dev +-- Round sampled numbers to integers +-- Use the integers to generate patters of the form 1 1 1... 1 -1 -1 -1 +-- which will have their Hamming distance normally distributed +getGaussianCluster :: MonadRandom m => Method -> Pattern -> Int -> Double -> Double -> m [Pattern] +getGaussianCluster method originPat size mean stdDev + | mean > fromIntegral patSize = error "the mean cannot be greater than the size of the pattern in getGaussianCluster" + | otherwise = do + normal_values <- replicateM size (normal mean stdDev) + return $ map encoding $ map round normal_values + where encoding x = V.fromList [ valueAtIndex y x | y <- [0 .. patSize - 1]] + patSize = V.length originPat + valueAtIndex y x = if (y <=x) then 1 else (smallerValue method) + smallerValue x = case x of + Hopfield -> -1 + _ -> 0 + +-- | @basinsGivenProbabilityT1 learning networkSize clusterSize p@ +-- Gets the average basin of attraction of a cluster of size @clusterSize@ +-- constructed using the T1 method given the flip probability @p@. +-- A hopfield network is trained (the type of training (Hebbian or Storkey) is +-- given by @learning@). +basinsGivenProbabilityT1 :: MonadRandom m => LearningType -> Int -> Int -> Double -> m Double +basinsGivenProbabilityT1 learning networkSize clusterSize p + = do + originPat <- randomSignVector networkSize + cluster <- getCluster Hopfield originPat clusterSize p + avgBasinsGivenPats learning cluster + + +-- | @experimentUsingT1 learning networkSize clusterSize@ +-- Gets the average basin of attraction obtained by iterating trough various +-- probabilities for flipping the bit when obtaining the cluster. +experimentUsingT1 :: MonadRandom m => LearningType -> Int -> Int -> m Double +experimentUsingT1 learning networkSize clusterSize + = do + basinAvgs <- mapM (basinsGivenProbabilityT1 learning networkSize clusterSize) probabilities + return $ average basinAvgs + where probabilities = [0.0, 0.1 .. 0.5] + +experimentUsingT1NoAvg :: MonadRandom m => LearningType -> Int -> Int -> m [(Double, Double)] +experimentUsingT1NoAvg learning networkSize clusterSize + = do + results <- mapM (basinsGivenProbabilityT1 learning networkSize clusterSize) probabilities + return $ zip probabilities results + where probabilities = [0.0, 0.1 .. 0.5] + + +------- + +basinsGivenProbabilityT1With2Clusters :: MonadRandom m => LearningType -> Int -> Int -> Double -> Double -> m (Double, Double) +basinsGivenProbabilityT1With2Clusters learning networkSize clusterSize p1 p2 = do + originPat1 <- randomSignVector networkSize + originPat2 <- randomSignVector networkSize + cluster1 <- getCluster Hopfield originPat1 clusterSize p1 + cluster2 <- getCluster Hopfield originPat2 clusterSize p2 + avg1 <- avgBasinsGivenPats learning cluster1 + avg2 <- avgBasinsGivenPats learning cluster2 + return $ (avg1, avg2) + + + +------- Experiments using Gaussian distributed patterns + +basinsGivenStdT2 :: MonadRandom m => LearningType -> Int -> Int -> Double -> Double -> m Double +basinsGivenStdT2 learning networkSize clusterSize mean std + = do + originPat <- randomSignVector networkSize + cluster <- getGaussianCluster Hopfield originPat clusterSize mean std + avgBasinsGivenPats learning cluster + + +experimentUsingT2 :: MonadRandom m => LearningType -> Int -> Int -> m Double +experimentUsingT2 learning networkSize clusterSize + = do + let mean = networkSize ./. (2 :: Int) + deviations = [0.0, 2.0 .. networkSize ./. (8 :: Int)] + basinAvgs <- mapM (basinsGivenStdT2 learning networkSize clusterSize mean) deviations + return $ average basinAvgs + +experimentUsingT2NoAvg :: MonadRandom m => LearningType -> Int -> Int -> m [(Double, Double)] +experimentUsingT2NoAvg learning networkSize clusterSize + = do + let mean = networkSize ./. (2 :: Int) + deviations = [0.0, 2.0 .. networkSize ./. (8 :: Int)] + basinAvgs <- mapM (basinsGivenStdT2 learning networkSize clusterSize mean) deviations + return $ zip deviations basinAvgs + + + +basinsGivenStdT2With2Clusters :: MonadRandom m => LearningType -> Int -> Int -> + Double -> Double -> Double -> Double -> m (Double, Double) +basinsGivenStdT2With2Clusters learning networkSize clusterSize mean1 mean2 std1 std2 = do + originPat1 <- randomSignVector networkSize + originPat2 <- randomSignVector networkSize + cluster1 <- getGaussianCluster Hopfield originPat1 clusterSize mean1 std1 + cluster2 <- getGaussianCluster Hopfield originPat2 clusterSize mean2 std2 + avg1 <- avgBasinsGivenPats learning cluster1 + avg2 <- avgBasinsGivenPats learning cluster2 + return $ (avg1, avg2) + + + + +--------------- General used functions, independent of method + +avgBasinsGivenPats :: MonadRandom m => LearningType -> [Pattern] -> m Double +avgBasinsGivenPats learning pats = do + basinSizes <- mapM (measurePatternBasin hopfield) pats + return $ average basinSizes + where hopfield = buildHopfieldData learning pats + + +-- Repeats an experiment for a single cluster, and averages the results obtained +-- in each of the experiments. +repeatExperiment :: MonadRandom m => (LearningType -> Int -> Int -> m Double) -> LearningType -> Int -> Int -> Int -> m Double +repeatExperiment experiment learning nrExperiments networkSize clusterSize + = liftM average $ replicateM nrExperiments (experiment learning networkSize clusterSize) diff --git a/src/Hopfield/Common.hs b/src/Hopfield/Common.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Common.hs @@ -0,0 +1,39 @@ +{-# LANGUAGE PatternGuards, ExistentialQuantification #-} + +module Hopfield.Common where + +-- This module contains data types and functions specific to the project +-- which are used for all different types of networks we support + +import Data.Vector (Vector) + +type Weights = Vector (Vector Double) +type Pattern = Vector Int +type Bias = Vector Double + +-- Data type used trought the project to choose a network to use +-- Boltzmann corresponds to the new method and CBoltzmann to the Classification +-- Boltzmann machine +data Method = Hopfield | Boltzmann | CBoltzmann + deriving (Eq, Enum, Ord, Show) + + +-- http://www.haskell.org/haskellwiki/Heterogenous_collections +data Showable = forall a . Show a => MkShowable a + +instance Show Showable + where showsPrec p (MkShowable a) = showsPrec p a + +pack :: Show a => a -> Showable +pack = MkShowable + + +packL :: Show a => [a] -> [Showable] +packL = map pack + + +-- flips a bit according to the method employed, as patterns +-- take different values if they are Hopfield or RBM. +flipBit :: Method -> Int -> Int +flipBit Hopfield x = - x +flipBit _ x = 1 - x diff --git a/src/Hopfield/Experiments/ClusterExperiments.hs b/src/Hopfield/Experiments/ClusterExperiments.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Experiments/ClusterExperiments.hs @@ -0,0 +1,85 @@ +{-# LANGUAGE ParallelListComp #-} + +module Hopfield.Experiments.ClusterExperiments where +-- Cluster Experiments which were performed by Federico + +import Control.Monad +import Control.Monad.Random +import Control.Parallel.Strategies + +import Hopfield.Clusters +import Hopfield.Hopfield +import Hopfield.Util + + +-- Data type which gives the type of the experiment +-- T1: bit flipping +-- T2: Gaussian distributed Hamming distance +data ExpType = T1 | T2 + deriving (Eq, Show, Read) + + +-- Runs one iteration of an experiment with 1 cluster +oneIteration1 :: ExpType -> LearningType -> Int -> Int -> Double -> Double -> Double -> Int-> [(Double, Double)] +oneIteration1 expType learnType networkSize clusterSize start stop p_step i + = zip cs values + where + f x = evalRand (evaluatedFunction x) (mkStdGen i) + unevaluated = map f values + cs = unevaluated `using` parList rdeepseq + values = [start, p_step .. stop] + evaluatedFunction = case expType of + T1 -> basinsGivenProbabilityT1 learnType networkSize clusterSize + T2 -> basinsGivenStdT2 learnType networkSize clusterSize (networkSize ./ 2.0) + + +-- Runs multiple iterations of an experiment with one cluster +-- Prints information to the user about the parameters of the experiment +performAndPrint1 :: ExpType -> LearningType -> Int -> Int -> Double -> Double -> Double -> Int -> IO () +performAndPrint1 expType learnType neurons clusterSize start stop step iterations = do + putStrLn $ "Experiment type" ++ show expType + putStrLn $ "Learning type " ++ show learnType + putStrLn $ "Only one clusters" + putStrLn $ "neurons " ++ show neurons ++ " cluster " ++ show clusterSize + putStrLn $ "performed for " ++ show iterations ++ " iterations" + mapM_ print $ map (oneIteration1 expType learnType neurons clusterSize start stop step) [0.. iterations] + + +-- Runs one iteration of an experiment with 2 clusters +oneIteration2 :: ExpType -> LearningType -> Int -> Int -> Double -> Double -> Double -> Double -> Int-> [(Double, (Double, Double))] +oneIteration2 expType learnType networkSize clusterSize val1 start2 stop2 p_step2 i + = zip values cs + where + f x = evalRand (evaluatedFunction x) (mkStdGen i) + unevaluated = map f values + cs = unevaluated `using` parList rdeepseq + values = [start2, start2 + p_step2 .. stop2] + evaluatedFunction = case expType of + T1 -> basinsGivenProbabilityT1With2Clusters learnType networkSize clusterSize val1 + T2 -> basinsGivenStdT2With2Clusters learnType networkSize clusterSize (networkSize ./ 2.0) (networkSize ./ 2.0) val1 + + +-- Runs multiple iterations of an experiment with 2 clusters +-- Prints information to the user about the parameters of the experiment +performAndPrint2 :: ExpType -> LearningType -> Int -> Int -> Double -> Double -> Double -> Double -> Int -> IO () +performAndPrint2 expType learnType neurons clusterSize val1 start2 stop2 step2 iterations = do + putStrLn $ "Experiment type " ++ show expType + putStrLn $ "Learning type " ++ show learnType + putStrLn $ "Two clusters" + putStrLn $ "neurons " ++ show neurons ++ " cluster " ++ show clusterSize + putStrLn $ "performed for " ++ show iterations ++ " iterations" + putStrLn $ "fixing the parameter(prob for T1 or std dev for T2) for the first cluster to be " ++ show val1 + putStrLn $ "varying the parameter for the second cluster between" ++ show start2 ++ "and " ++ show stop2 + seeds <- replicateM iterations $ getRandomR (0 :: Int, 1000 :: Int) + mapM_ print $ map (oneIteration2 expType learnType neurons clusterSize val1 start2 stop2 step2) seeds + + +-- Called from the main of apps/ExperimentsMain.hs when the first argument of +-- the executable is 'cluster' +run :: [String] -> IO () +run args = do + + case args of + ("1": t : l: n : c : start : stop : step: iterations: _)-> performAndPrint1 (read t) (read l) (read n) (read c) (read start) (read stop) (read step) (read iterations) + ("2": t : l: n : c : fixed : start : stop : step: iterations: _)-> performAndPrint2 (read t) (read l) (read n) (read c) (read fixed) (read start) (read stop) (read step) (read iterations) + _ -> error "invalid arguments" diff --git a/src/Hopfield/Experiments/Experiment.hs b/src/Hopfield/Experiments/Experiment.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Experiments/Experiment.hs @@ -0,0 +1,94 @@ +{-# LANGUAGE ParallelListComp #-} + +module Hopfield.Experiments.Experiment where + +import Control.Monad (replicateM) +import Control.Monad.Random +import Test.QuickCheck +import Test.QuickCheck.Gen (unGen) + +import Hopfield.Clusters +import Hopfield.Common +import Hopfield.Experiments.ExperimentUtil +import Hopfield.Hopfield +import Hopfield.Measurement +import Hopfield.SuperAttractors +import Hopfield.TestUtil (Type(H), patternGen) +import Hopfield.Util + + +genIO :: Gen a -> IO a +genIO g = do + rndInt <- randomIO + stdGen <- getStdGen + return $ unGen g stdGen rndInt + + +errorHeader :: String +errorHeader = "Degree\tExpected error" + + +basinHeader :: String +basinHeader = "Degree\tBasin size" + + +main :: IO () +main = do + + let n = 100 -- number of neurons + numRandoms = 8 -- number of random patterns to include + maxDegree = 32 -- maximum degree of super attractor + + + -- The super attractor - primary care giver + originPat <- genIO $ patternGen H n + + -- Sample random patterns with Hamming distance between 25-75% from origin + -- This is to ensure that this is a pure super attractor experiment + -- and not a cluster one! + let minHamming = round $ n .* (0.25 :: Double) + maxHamming = round $ n .* (0.75 :: Double) + dist = hammingDistribution n (minHamming, maxHamming) + + randomPats <- replicateM numRandoms $ sampleHammingRange originPat dist + + + let pats = originPat:randomPats + p = length pats + originIndex = 0 -- index of main pattern + degrees = powersOfTwo maxDegree + patCombiner = oneSuperAttr + + + putStrLn $ unwords [show n, "neurons.", "Super attractor plus", show numRandoms, "random patterns.\n"] + + + -- Compute probability of error + doErrorProb n p degrees + + + -- Compute hamming distance + doHamming originPat randomPats "origin" "random" + + + putStrLn "Building networks..." + let nets = buildNetworks pats degrees Hebbian patCombiner + + + --Check if pattern is fixed. + doCheckFixed (zip degrees nets) originIndex "degrees" + + + putStrLn "Measuring basins of attraction" + let results = measureMultiBasins measurePatternBasin nets originPat + + putStrLn basinHeader + printMList results [ \r -> attachLabel [pack d, pack r] | d <- degrees ] + + -- putStrLn "T1 experiment with 1 cluster" + -- putStrLn $ show $ evalRand (repeatExperiment experimentUsingT1 Storkey 1 50 8) (mkStdGen 1) + + putStrLn "T1 experiment with 1 cluster with no average but lists" + + let avgs = replicate 10 $ experimentUsingT1NoAvg Hebbian 100 10 + printMList avgs (replicate 10 show) diff --git a/src/Hopfield/Experiments/Experiment2SuperAttractors.hs b/src/Hopfield/Experiments/Experiment2SuperAttractors.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Experiments/Experiment2SuperAttractors.hs @@ -0,0 +1,94 @@ +{-# LANGUAGE ParallelListComp #-} + +-- | Performs experiments with two super attractors. +module Hopfield.Experiments.Experiment2SuperAttractors where + +import Control.Monad (replicateM) +import Control.Monad.Random +import Test.QuickCheck +import Test.QuickCheck.Gen (unGen) + +import Hopfield.Common +import Hopfield.Experiments.ExperimentUtil +import Hopfield.Hopfield (LearningType (..)) +import Hopfield.Measurement +import Hopfield.SuperAttractors +import Hopfield.TestUtil (Type(H), patternGen) +import Hopfield.Util + + +genIO :: Gen a -> IO a +genIO g = do + rndInt <- randomIO + stdGen <- getStdGen + return $ unGen g stdGen rndInt + + +basinHeader :: String +basinHeader = "Degree\tOrigin basin\tNew basin" + + +main :: IO () +main = do + + let n = 100 -- number of neurons + numRandoms = 8 -- number of random patterns to include + maxDegree = 32 -- maximum degree of second super attractor + fstDegree = 8 -- (fixed) degree of first super attractor + + + -- The first super attractor - primary care giver + originPat <- genIO $ patternGen H n + + -- Sample random patterns with Hamming distance between 25-75% from origin + -- This is to ensure that this is a pure super attractor experiment + -- and not a cluster one! + let minHamming = round $ n .* (0.25 :: Double) + maxHamming = round $ n .* (0.75 :: Double) + dist = hammingDistribution n (minHamming, maxHamming) + + randomPats <- replicateM numRandoms $ sampleHammingRange originPat dist + + -- The second super attractor - retraining + newPat <- sampleHammingRange originPat dist + + + + let pats = originPat:newPat:randomPats + originIndex = 0 -- index of main pattern + newIndex = fstDegree + 1 -- index of new pattern + degrees = powersOfTwo maxDegree + patCombiner = twoSuperAttrOneFixed fstDegree + + + putStrLn $ unwords [show n, "neurons.", "Two Super attractors plus", show numRandoms, "random patterns.\n"] + + + -- Check hamming distances + doHamming originPat randomPats "origin" "random" + doHamming newPat randomPats "new" "random" + + + putStrLn "Building networks...\n" + let nets = buildNetworks pats degrees Hebbian patCombiner + + + --Check if patterns are fixed. + putStrLn "Checking original pattern" + doCheckFixed (zip degrees nets) originIndex "degrees" + putStrLn "Checking new pattern" + doCheckFixed (zip degrees nets) newIndex "degrees" + + + putStrLn "Measuring basins of attraction of origin" + let resultsOrigin = measureMultiBasins measurePatternBasin nets originPat + let resultsNew = measureMultiBasins measurePatternBasin nets newPat + + + + let results = zipWith (\a b -> sequence [a, b]) resultsOrigin resultsNew + printResults d rs = attachLabel $ [pack d] ++ map pack rs + + putStrLn basinHeader + printMList results [ printResults d | d <- degrees ] + diff --git a/src/Hopfield/Experiments/ExperimentUtil.hs b/src/Hopfield/Experiments/ExperimentUtil.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Experiments/ExperimentUtil.hs @@ -0,0 +1,49 @@ +{-# LANGUAGE ParallelListComp #-} + +module Hopfield.Experiments.ExperimentUtil where + +import qualified Data.Vector as V + +import Hopfield.Analysis (computeErrorSuperAttractorNumbers) +import Hopfield.Common +import Hopfield.Hopfield +import Hopfield.Measurement +import Hopfield.SuperAttractors (Degree) +import Hopfield.Util + + +-- Measure hamming distance from p to each of ps +doHamming :: Pattern -> [Pattern] -> String -> String -> IO () +doHamming p ps pName psName = do + let msg = unwords ["Hamming distance between", pName, "pattern and", + psName, "patterns:"] + putStrLn msg + let n = V.length p + hammingDists = map (hammingDistance p) ps + hammingPct = map (./. n) hammingDists :: [Double] + putStrLn $ prettyList hammingDists + putStrLn $ toPercents hammingPct ++ "\n" + + +-- Check if pattern is a fixed fixed points +doCheckFixed :: Show a => [ (a, HopfieldData) ] -> Int -> String -> IO () +doCheckFixed pairs index labelsName = do + let patErrs = [ label | (label, net) <- pairs, not $ checkFixed net index] + msg = unwords ["WARNING: The following", labelsName, + "have produced networks where the pattern is NOT a fixed point:\n"] + + if not $ null patErrs + then putStrLn $ msg ++ prettyList patErrs ++ "\n" + else putStrLn "Pattern is always a fixed point\n" + + +-- Compute probabilities of error - i.e. pattern not fixed +doErrorProb :: Int -> Int -> [Degree] -> IO () +doErrorProb n p degrees = do + + putStrLn $ "Expected network errors: " + + let expErrs = [ computeErrorSuperAttractorNumbers d p n | d <- degrees ] + errorHeader = "Degree\tExpected error" + + putStrLn $ attachLabels errorHeader $ [packL degrees, packL expErrs] diff --git a/src/Hopfield/Experiments/SmallExperiments.hs b/src/Hopfield/Experiments/SmallExperiments.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Experiments/SmallExperiments.hs @@ -0,0 +1,30 @@ +{-# LANGUAGE ParallelListComp #-} + +module Hopfield.Experiments.SmallExperiments where + +-- Module use to perform small experiments that prove that Storkey learning +-- has a bigger basin size than Hebbian learning + +import Control.Applicative +import Control.Monad + +import Hopfield.Clusters +import Hopfield.Hopfield +import Hopfield.Util + +_REPETITIONS :: Int +_REPETITIONS = 10 + +-- Experiments are performed using the T2 method +-- (learning type, number of neurons, cluster size, mean of cluster, std dev) +runs :: [(LearningType, Int, Int, Double, Double)] +runs = [ (Hebbian, 50, 6, 25, 5) + , (Storkey, 50, 6, 25, 5) + , (Hebbian, 50, 4, 25, 10) + , (Storkey, 50, 4, 25, 10) + ] + +main :: IO () +main = do + forM_ runs $ \(method, a, b, c, d) -> + print =<< ((average <$> replicateM _REPETITIONS (basinsGivenStdT2 method a b c d)) :: IO Double) diff --git a/src/Hopfield/Hopfield.hs b/src/Hopfield/Hopfield.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Hopfield.hs @@ -0,0 +1,315 @@ +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} + +-- | Base Hopfield model, providing training and running. +module Hopfield.Hopfield ( + Pattern + , Weights + , LearningType (Hebbian, Storkey) + -- * Hopfield data structure + , HopfieldData () + , weights + , patterns + , buildHopfieldData + -- * Running + , update + , addPatterns + , repeatedUpdate + , updateChain + , matchPattern + , computeH + -- * Energy + , energy +) where + +import Control.Monad +import Control.Monad.Random (MonadRandom) +import Data.Maybe +import Data.Vector ((!)) +import qualified Data.Vector as V +import Data.Vector.Generic.Mutable (write) + +import Hopfield.Common +import Hopfield.Util + + + +data LearningType = Hebbian | Storkey deriving (Eq, Show, Read) + +--make Hopefield data implement show +-- | Encapsulates the network weights together with the patterns that generate +-- it with the patterns which generate it +data HopfieldData = HopfieldData { + weights :: Weights -- ^ the weights of the network + , patterns :: [Pattern] -- ^ the patterns which were used to train it +} deriving (Show) + + +-- | Checks if weights and pattern given to the function satisfy their constraints, +-- if yes, calling the function, otherwise erroring out. +-- Usage: `checkWsPat (functionTakingWeightsAndPattern)`. +checkWsPat :: (Weights -> Pattern -> a) -> Weights -> Pattern -> a +checkWsPat f ws pat + | Just e <- validWeights ws = error e + | Just e <- validPattern pat = error e + | Just e <- validWeightsPatternSize ws pat = error e + | otherwise = f ws pat + + +-- | @update weights pattern@: Applies the update rule on @pattern@ for the +-- first updatable neuron given the Hopfield network (represented by @weights@). +-- +-- Pre: @length weights == length pattern@ +update :: MonadRandom m => Weights -> Pattern -> m (Maybe Pattern) +update = checkWsPat update_ + + +-- | @repeatedUpdate weights pattern@: Performs repeated updates on the given +-- pattern until it reaches a stable state with respect to the Hopfield network +-- (represented by @weights@). +-- Pre: @length weights == length pattern@ +repeatedUpdate :: (MonadRandom m) => Weights -> Pattern -> m Pattern +repeatedUpdate = checkWsPat repeatedUpdate_ + + +-- | Computes the weighted sum of current neuron values, which will give us +-- the value of the neuron (by taking the sign) +computeH :: Weights -> Pattern -> Int -> Int +computeH ws pat i = checkWsPat (\w p -> computeH_ w p i) ws pat + + +-- | @energy weights pattern@: Computes the energy of a pattern given a Hopfield +-- network (represented by @weights@). +-- Pre: @length weights == length pattern@ +energy :: Weights -> Pattern -> Double +energy = checkWsPat energy_ + + + +-- | @buildHopfieldData patterns@: Takes a list of patterns and +-- builds a Hopfield network (by training) in which these patterns are +-- stable states. The result of this function can be used to run a pattern +-- against the network, by using 'matchPattern'. +buildHopfieldData :: LearningType -> [Pattern] -> HopfieldData +buildHopfieldData _ [] = error "Train patterns are empty" +buildHopfieldData learningType pats + | first_len == 0 + = error "Cannot have empty patterns" + | any (\x -> V.length x /= first_len) pats + = error "All training patterns must have the same length" + | otherwise + = HopfieldData (trainingFunction pats) pats + where + first_len = V.length (head pats) + trainingFunction = case learningType of + Hebbian -> train + Storkey -> trainStorkey + + +-- | @train patterns@: Trains and constructs network given a list of patterns +-- which are used to build the weight matrix. As a consequence, they will be +-- stable points in the network (by construction). +train :: [Pattern] -> Weights +train pats = vector2D ws + -- No need to check pats ws size, buildHopfieldData does it + where + ws = [ [ w i j ./. n | j <- [0 .. n-1] ] | i <- [0 .. n-1] ] + w i j + | i == j = 0 + | otherwise = sum [ (pat ! i) * (pat ! j) | pat <- pats ] + n = V.length (head pats) + + +-- | See `computeH`. +computeH_ :: Weights -> Pattern -> Int -> Int +computeH_ ws pat i = {-# SCC "computeHall" #-} if weighted >= 0 then 1 else -1 + where + weighted :: Double + wss = ws ! i + weighted = go 0 0.0 + go :: Int -> Double -> Double + go !j !s | j == p = s + | otherwise = let w = wss `V.unsafeIndex` j + x = if pat `V.unsafeIndex` j > 0 then w + else -w + in go (j+1) (s+x) + + p = {-# SCC "computeHvlength" #-} V.length pat + + + +-- | See `update`. +-- The update is done by finding a neuron that will change its value given the +-- current state. The search for this neuron is done in a random manner: +-- pick up a random neuron, check if it is updatable: if so, update the pattern +-- by updating this neuron. If not, continue until an updatable neuron is found. +-- (Note: Initially the update was performed by obtaining a list of all +-- updatable neurons and then picking a random one. The current method is 2 times +-- faster) +update_ :: MonadRandom m => Weights -> Pattern -> m (Maybe Pattern) +update_ ws pat = do + randomIndices <- shuffle . toArray $ [0 .. V.length pat - 1] + -- TODO avoid Array -> List -> Vector conversion + return $ case firstUpdatable (V.fromList randomIndices) of + Nothing -> Nothing + Just index -> Just $ flipAtIndex pat index + where + firstUpdatable indices = go 0 + where + go n + | n == V.length pat = Nothing + | pat ! i /= computeH_ ws pat i = Just i + | otherwise = go (n+1) + where i = indices ! n + + flipAtIndex vec index = let val = vec ! index -- seq only brings small saving here + in val `seq` V.modify (\v -> write v index (-val)) vec + + +-- | See `repeatedUpdate`. +repeatedUpdate_ :: (MonadRandom m) => Weights -> Pattern -> m Pattern +repeatedUpdate_ ws pat = repeatUntilNothing (update_ ws) pat + + +-- | @matchPatterns hopfieldData pattern@: +-- Computes the stable state of a pattern given a Hopfield network(represented +-- by @weights@) and tries to find a match in a list of patterns which are +-- stored in @hopfieldData@. +-- Returns: +-- +-- The index of the matching pattern in @patterns@, if a match exists +-- The converged pattern (the stable state), otherwise +-- +-- Pre: @length weights == length pattern@ +matchPattern :: MonadRandom m => HopfieldData -> Pattern -> m (Either Pattern Int) +matchPattern (HopfieldData ws pats) pat = do + converged_pattern <- repeatedUpdate_ ws pat + return $ findInList pats converged_pattern + + +-- | Like `repeatedUpdate`, but collecting all patterns until convergence. +-- The last pattern in the list is the converged pattern. +-- The argument pattern is NOT prepended to the result list. +-- +-- POST: The returned list is not empty. +updateChain :: (MonadRandom m) => HopfieldData -> Pattern -> m [Pattern] +updateChain (HopfieldData ws _pats) pat + | Just e <- validPattern pat = error e + | otherwise = (pat:) `liftM` unfoldrSelfM (update_ ws) pat + + +-- | Stores patterns in an already trained network. One has to ensure that this +-- function is not over used, as this will decrease the capacity of the network. +addPatterns :: LearningType -> HopfieldData -> [Pattern] -> HopfieldData +addPatterns learning (HopfieldData ws pats) addedPats + | any (isJust . validPattern) addedPats = error "invalid patterns in addMultiplePatterns" + | any (isJust . validWeightsPatternSize ws) addedPats = error "pattern does not match weights in addMultiplePatterns" + | otherwise = HopfieldData new_ws (pats ++ addedPats) + where new_ws = foldl (updateWeightsGivenNewPattern learning) ws addedPats + + +-- Updates the weight matrix when a new pattern is stored in the network +updateWeightsGivenNewPattern :: LearningType -> Weights -> Pattern -> Weights +updateWeightsGivenNewPattern Storkey ws pat = updateWeightsStorkey ws pat +updateWeightsGivenNewPattern Hebbian ws pat = vector2D updated_ws + where updated_ws = [ [ws ! i ! j + (pat ! i * pat ! j) ./. n | j <- neurons ] | i <- neurons] + n = V.length ws - 1 + neurons = [0 .. n] + + +-- | See `energy`. +energy_ :: Weights -> Pattern -> Double +energy_ ws pat = s / (-2.0) + where + p = V.length pat + w i j = ws ! i ! j + x i = pat ! i + s = sum [ w i j *. (x i * x j) | i <- [0 .. p-1], j <- [0 .. p-1] ] + + +-- | Checks if a pattern consists of only 1s and -1s. +-- Returns @Nothing@ on success, an error string on failure. +validPattern :: Pattern -> Maybe String +validPattern pat = case [ x | x <- V.toList pat, not (x == 1 || x == -1) ] of + [] -> Nothing + x:_ -> Just $ "Pattern contains invalid value " ++ show x + + +-- | @validWeightsPatternSize weights pattern@ +-- Returns an error string in a Just if the @pattern@ is not compatible +-- with @weights@ and Nothing otherwise. +validWeightsPatternSize :: Weights -> Pattern -> Maybe String +validWeightsPatternSize ws pat + | V.length ws /= V.length pat = Just "Pattern size must match network size" + | otherwise = Nothing + + +-- Checks the validity of a weight matrix by ensuring: +-- * It is non-empty +-- +-- * It is square +-- +-- * It is symmetric +-- +-- * All diagonal elements must be zero +-- These checks hold for both Hebbian and Storkey. +validWeights :: Weights -> Maybe String +validWeights ws + | n == 0 + = Just "Weight matrix must be non-empty" + | any (\x -> V.length x /= n) $ V.toList ws + = Just "Weight matrix has to be a square matrix" + | any (/= 0) [ ws ! i ! i | i <- [0..n-1] ] + = Just "Weight matrix first diagonal must be zero" + | not $ and [ abs( (ws ! i ! j) - (ws ! j ! i) ) < 0.0001 | i <- [0..n-1], j <- [0..n-1] ] + = Just "Weight matrix must be symmetric" + | null [ abs (ws ! i ! j) > 1 | i <- [0..n-1], j <- [0..n-1] ] + = Just "Weights should be between (-1, 1)" + | otherwise = Nothing + where + n = V.length ws + + +-- Storkey training provides advantages for the Hopfield network as +-- it gives it bigger capacity and higher basins of attraction. +-- For more details see: +-- http://homepages.inf.ed.ac.uk/amos/publications/Storkey1997IncreasingtheCapacityoftheHopfieldNetworkwithoutSacrificingFunctionality.pdf + + +-- | @storkeyHiddenSum ws pat i j@ computes the value at indices @i@ @j@ in the +-- hidden matrix which is used for updating in the weight matrix during trainig +-- given the training pattern @pat@. +storkeyHiddenSum :: Weights -> Pattern -> Int -> Int -> Double +storkeyHiddenSum ws pat i j + = sum [ ws ! i ! k *. (pat ! k) | k <- [0 .. n - 1] , k /= i , k /= j] + where n = V.length ws + +-- | @updateWeightsGivenIndicesStorkey ws pat i j@ computes the new value at +-- indices @i@ @j@ of the weights matrix for the training iteration of +-- pattern @pat@. +updateWeightsGivenIndicesStorkey :: Weights -> Pattern -> Int -> Int -> Double +updateWeightsGivenIndicesStorkey ws pat i j + | i == j = 0.0 + | otherwise = ws ! i ! j + (1 :: Int) ./. n * (fromIntegral (pat ! i * (pat ! j)) - h j i *. (pat ! i) - h i j *. (pat ! j)) + where n = V.length ws + h = storkeyHiddenSum ws pat + + +-- | @updateWeightsStorkey ws pat@ updates the weights matrix, given training +-- instance @pat@. +updateWeightsStorkey :: Weights -> Pattern -> Weights +updateWeightsStorkey ws pat + = vector2D [ [ updateWeightsGivenIndicesStorkey ws pat i j | j <- [0 ..n - 1] ] | i <- [0 ..n - 1] ] + where n = V.length ws + + +-- | @trainStorkey pats@ trains the Hopfield network by computing the weights +-- matrix by iterating trough all training instances (@pats@) and updating the +-- weights according to the Storkey learning rule. +trainStorkey :: [Pattern] -> Weights +-- No need to check pats ws size, buildHopfieldData does it +trainStorkey pats = foldl updateWeightsStorkey start_ws pats + where start_ws = vector2D $ replicate n $ replicate n 0 + n = V.length $ head pats + diff --git a/src/Hopfield/Images/ConvertImage.hsc b/src/Hopfield/Images/ConvertImage.hsc new file mode 100644 --- /dev/null +++ b/src/Hopfield/Images/ConvertImage.hsc @@ -0,0 +1,40 @@ +{-# LANGUAGE CPP, ForeignFunctionInterface #-} + +module Hopfield.Images.ConvertImage ( + loadPicture +, CBinaryPattern (..) +) where + +import Data.Word +import Foreign.C +import Foreign.Ptr +import Foreign.Storable +import Foreign.Marshal.Array + +#include "Images/convertImage.h" + +-- From: http://www.haskell.org/haskellwiki/FFI_cook_book +#let alignment t = "%lu", (unsigned long) offsetof(struct { char x__; t (y__); }, y__) + +data CBinaryPattern = CBinaryPattern { + cPatternSize :: Word32 + , cPattern :: [Word32] +} deriving (Eq, Show) + + +foreign import ccall "convertImage.h load_picture" load_picture :: CString -> CInt -> CInt -> Ptr CBinaryPattern + +instance Storable CBinaryPattern where + alignment _ = #{alignment binary_pattern_t} + sizeOf _ = #{size binary_pattern_t} + peek ptr = do s <- #{peek binary_pattern_t, size} ptr + pattern_ptr <- #{peek binary_pattern_t, pattern} ptr + pattern_01s <- peekArray (fromIntegral s) pattern_ptr + return $ CBinaryPattern s pattern_01s + poke _ptr (CBinaryPattern _s _p) = error "Storable CBinaryPattern: poke not implemented" + + +loadPicture :: String -> Int -> Int-> IO CBinaryPattern +loadPicture path w h = do + cpath <- newCString path + peek (load_picture cpath (fromIntegral w) (fromIntegral h)) diff --git a/src/Hopfield/Images/convertImage.c b/src/Hopfield/Images/convertImage.c new file mode 100644 --- /dev/null +++ b/src/Hopfield/Images/convertImage.c @@ -0,0 +1,93 @@ +#include "convertImage.h" +#include +#include + +void ThrowWandException(MagickWand *wand) + { + char *description; + ExceptionType severity; + description=MagickGetException(wand,&severity); + (void) fprintf(stderr,"%s %s %lu %s\n",GetMagickModule(),description); + description=(char *) MagickRelinquishMemory(description); + exit(-1); + } + +/* converts a list of doubles to binary + flattens the matrix to a vector */ +binary_pattern_t * +mapToBinary(double** pattern, int width, int height){ + binary_pattern_t* binaryPattern = (binary_pattern_t *) malloc(sizeof(*binaryPattern)); + const int size = width * height; + binaryPattern->size = size; + binaryPattern->pattern = (uint32_t *) malloc(sizeof(*binaryPattern->pattern) * size); + + int i=0; + + for(int h = 0; h < height; h++) + for(int w = 0; w < width; w++) + { + binaryPattern->pattern[i] = pattern[w][h] < 0.5 ? 0 : 1; + i++; + } + + return binaryPattern; +} + +/* loads a picture from the inputImg filepath, scales it to the specified + width and height, and then converts it to a binary pattern, usable by + the Hopfield Network */ +binary_pattern_t * +load_picture(char* inputImg, size_t width, size_t height) +{ + /* load a picture in the "wand" */ + MagickWand *mw = NULL; + + MagickWandGenesis(); + + /* Create a wand */ + mw = NewMagickWand(); + + /* Read the input image */ + MagickBooleanType retVal = MagickReadImage(mw, inputImg); + + if (retVal == MagickFalse) + ThrowWandException(mw); + + PixelWand** pixels; + /* rescale the image */ + int resizeSuccess = MagickResizeImage(mw, height, width, LanczosFilter, 0); + + if (!resizeSuccess) + { + printf("resize failed\n"); + exit(1); + } + + PixelIterator* pixelIt = NewPixelIterator(mw); + double** outputPattern = (double**) malloc (sizeof(double*) * width); + for(size_t i=0; i < width; i++) + { + outputPattern[i] = (double*) malloc(sizeof(double) * height); + } + + long y; + /* get pixel grayscale values */ + for (y=0; y < (long) height; y++) + { + size_t iter_width; + pixels=PixelGetNextIteratorRow(pixelIt,&iter_width); + assert (iter_width == width); + for (long x=0; x < (long) iter_width; x++) { + outputPattern[x][y] = (PixelGetRed(pixels[x]) + + PixelGetGreen(pixels[x]) + PixelGetBlue(pixels[x]))/3; + } + } + + /* Tidy up */ + if(mw) mw = DestroyMagickWand(mw); + + MagickWandTerminus(); + + /* since outputPattern is a list of doubles, convert it to a list of + binary values */ + return mapToBinary(outputPattern, width, height); +} diff --git a/src/Hopfield/Measurement.hs b/src/Hopfield/Measurement.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Measurement.hs @@ -0,0 +1,116 @@ +-- | Functions to measure various properties of a network +module Hopfield.Measurement ( + -- * Basin of attraction + BasinMeasure + , hammingDistribution + , sampleHammingRange + , sampleHammingDistance + , samplePatternRing + , samplePatternBasin + , measurePatternBasin + -- * Fixed point errors + , checkFixed + , measureError +) where + +import Control.Monad (liftM, replicateM) +import Control.Monad.Random (MonadRandom) +import Data.List +import Data.Maybe +import qualified Data.Vector as V +import Math.Combinatorics.Exact.Binomial (choose) +import Numeric.Probability.Distribution (Spread, relative) +import Numeric.Probability.Random (T, pick) + +import Hopfield.Hopfield +import Hopfield.Util ((./.), toArray, shuffle, runT) + + +-- A function computing some measure of a pattern's basin in the given network +type BasinMeasure m a = HopfieldData -> Pattern -> m a + + +-- ----------------------------------------------------------------------------- +-- Functions relating to measuring a pattern's basin of attraction + + +-- Create a probability distribution for Hamming distances in the given range +hammingDistribution :: Int -> (Int, Int) -> T Int +hammingDistribution n (mini, maxi) = pick $ dist rs + where + dist = relative probs :: Spread Double Int + probs = [ fromIntegral $ n `choose` r | r <- rs] + rs = [mini..maxi] + + +-- Sample a pattern in the Hamming distance range specified by dist +sampleHammingRange :: MonadRandom m => Pattern -> T Int -> m Pattern +sampleHammingRange pat dist = do + r <- runT dist + (sample:_) <- sampleHammingDistance pat r 1 + return sample + + +-- Samples patterns of hamming distance r of the given pattern +sampleHammingDistance :: MonadRandom m => Pattern -> Int -> Int -> m [Pattern] +sampleHammingDistance pat r numSamples + = liftM (map (V.fromList . multByPat)) coeffSamples + where + n = V.length pat + basePerm = toArray $ replicate r (-1) ++ replicate (n-r) 1 + coeffSamples = replicateM numSamples $ shuffle basePerm + multByPat coeffs = zipWith (*) coeffs (V.toList pat) + + +-- Percentage of sampled patterns in the ring of 'pat' which converge to 'pat' +-- pre: pattern of same size as network + +-- A pattern ring of radius 'r' around 'pat' is the set of states with hamming +-- distance 'r' from 'pat'. +samplePatternRing :: MonadRandom m => HopfieldData -> Pattern -> Int -> m Double +samplePatternRing hs pat r = do + samples <- sampleHammingDistance pat r 100 + convergedPatterns <- mapM (repeatedUpdate $ weights hs) samples + let numConverging = length $ filter (==pat) convergedPatterns + + return $ numConverging ./. (length samples) + + +-- Percentage convergence for each ring of 'pat' (excluding the trivial ring 0) +-- pre: pattern of same size as network +samplePatternBasin :: (MonadRandom m) => BasinMeasure m [Double] +samplePatternBasin hs pat = mapM (samplePatternRing hs pat) [1..n] + where + n = V.length pat + + +-- Measures pattern's basin of attraction using the Storkey-Valabregue method +-- pre: pattern of same size as network +measurePatternBasin :: (MonadRandom m) => BasinMeasure m Int +measurePatternBasin hs pat = do + t_mus <- samplePatternBasin hs pat + return $ fromMaybe n $ findIndex (<0.9) t_mus + where + n = V.length pat + + +-- ----------------------------------------------------------------------------- +-- Functions relating to measuring errors in a network + +compTerm :: HopfieldData -> Int -> Int -> Int +compTerm hs index n = - (pat V.! n) * (computeH (weights hs) pat n - pat V.! n) + where pat = (patterns hs) !! index + + +checkFixed :: HopfieldData -> Int -> Bool +checkFixed hs index = all (\x -> compTerm hs index x <= 1) [0.. V.length ((patterns hs) !! index) - 1] + + +-- | @measureError hopfield@: Measures the percentage of patterns in the network +-- which are NOT fixed points. That is, it measures the *actual* error +measureError :: HopfieldData -> Double +measureError hs = num_errors ./. num_pats + where + fixed_points = map (checkFixed hs) [0..num_pats-1] + num_errors = length $ filter not fixed_points + num_pats = length $ patterns hs diff --git a/src/Hopfield/SuperAttractors.hs b/src/Hopfield/SuperAttractors.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/SuperAttractors.hs @@ -0,0 +1,133 @@ +{-# LANGUAGE PatternGuards #-} + +-- Provides functions to construct and measure networks with super attractors +-- with varying pattern degrees +module Hopfield.SuperAttractors where + +import Control.Monad.Random (MonadRandom) +import qualified Data.Vector as V + +import Hopfield.Hopfield +import Hopfield.Measurement + + +-- Degree of a pattern is the number of instances it has in a network +type Degree = Int + + +-- A function combining some input and given degree into patterns for a network +type PatternCombiner a = a -> Degree -> [Pattern] + + +-- Produces all powers of two <= ceil +powersOfTwo :: Degree -> [Degree] +powersOfTwo ceil = takeWhile (<= ceil) $ iterate (*2) 1 + + +-- For each degree in 'ds', builds a network combining the degree and the list +-- of patterns (or some variant) 'ps' using the given function 'combine' +buildNetworks :: a -> [Degree] -> LearningType-> PatternCombiner a -> [HopfieldData] +buildNetworks ps ds learnType combine + = [ buildHopfieldData learnType $ combine ps d | d <- ds ] + + +-- ----------------------------------------------------------------------------- +-- Combine functions. 'buildNetworks' uses these to build super attractors + +-- Replicates the first pattern k times. +oneSuperAttr :: PatternCombiner [Pattern] +oneSuperAttr [] _ = [] +oneSuperAttr (p: ps) k = replicate k p ++ ps + + +-- Replicates the first pattern j times, and the second pattern k times +twoSuperAttrOneFixed :: Degree -> PatternCombiner [Pattern] +twoSuperAttrOneFixed j (pa:pb: ps) k = replicate j pa ++ replicate k pb ++ ps +twoSuperAttrOneFixed _ _ _ = [] + + +-- Replicates each pattern k times. +allSuperAttr :: PatternCombiner [Pattern] +allSuperAttr ps k = concatMap (replicate k) ps + + +-- Aggregate list of combiner functions of input [Pattern] into a single +-- combiner function of input [[Pattern]] +aggregateCombiners :: [PatternCombiner [Pattern]] -> PatternCombiner [[Pattern]] +aggregateCombiners combiners patList degree + | length combiners /= length patList + = error "Number of [Pattern] in list must match number of functions" + | otherwise + = concat $ zipWith ($) funcs patList + where + funcs = map (($ degree) . flip) combiners + +-- ----------------------------------------------------------------------------- +-- Experiments to measure super attractors + +-- Training (pre) patterns +p1, p2 :: Pattern +p1 = V.fromList [1,1,1,-1,-1,1,1,-1,1,-1] +p2 = V.fromList [-1,-1,1,1,-1,-1,1,-1,-1,1] + + +-- Retraining (post) patterns +q1 :: Pattern +q1 = V.fromList [1,-1,-1,-1,1,-1,-1,1,1,1] + + +-- Networks with first pattern as a super attractor +oneSuperNets :: LearningType -> [HopfieldData] +oneSuperNets learnType = buildNetworks ps degrees learnType oneSuperAttr + where + ps = [p1,p2] + degrees = powersOfTwo $ V.length $ head ps + + +-- Networks with all patterns as (equal) super attractors +allSuperNets :: LearningType -> [HopfieldData] +allSuperNets learnType = buildNetworks ps degrees learnType allSuperAttr + where + ps = [p1,p2] + degrees = powersOfTwo $ V.length $ head ps + + + +-- Convenience function for building networks with multiple training phases +buildMultiPhaseNetwork :: LearningType -> [PatternCombiner [Pattern]] -> [HopfieldData] +buildMultiPhaseNetwork learnType combFuncs = buildNetworks patList degrees learnType aggComb + where + patList = [ [p1,p2], [q1] ] + degrees = powersOfTwo $ (V.length . head . head) patList + aggComb = aggregateCombiners combFuncs + + +retrainNormalWithOneSuper :: LearningType -> [HopfieldData] +retrainOneSuperWithNormal :: LearningType -> [HopfieldData] +retrainOneSuperWithOneSuper :: LearningType -> [HopfieldData] +retrainAllSuperWithNormal :: LearningType -> [HopfieldData] +retrainAllSuperWithOneSuper :: LearningType -> [HopfieldData] + +-- A normal network (i.e. no super attractor) retrained with one super attractor +retrainNormalWithOneSuper l = buildMultiPhaseNetwork l [const, oneSuperAttr] + +-- A network with one super attractor retrained with a normal pattern (i.e. a +-- non-super attractor) +retrainOneSuperWithNormal l = buildMultiPhaseNetwork l [oneSuperAttr, const] + +-- A network with one super attractor retrained with another super attractor +retrainOneSuperWithOneSuper l = buildMultiPhaseNetwork l [oneSuperAttr, oneSuperAttr] + +-- A network with all super attractors retrained with a normal pattern (i.e. a +-- non-super attractor) +retrainAllSuperWithNormal l = buildMultiPhaseNetwork l [allSuperAttr, const] + +-- A network with all super attractors retrained with another super attractor +retrainAllSuperWithOneSuper l = buildMultiPhaseNetwork l [allSuperAttr, oneSuperAttr] + + + +-- Measure basin of multiple networks +measureMultiBasins :: MonadRandom m => BasinMeasure m a -> [HopfieldData] -> Pattern -> [m a] +measureMultiBasins measureBasin hs p = map (\h -> measureBasin h p) hs + diff --git a/src/Hopfield/TestUtil.hs b/src/Hopfield/TestUtil.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/TestUtil.hs @@ -0,0 +1,256 @@ +{-# LANGUAGE ParallelListComp #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module Hopfield.TestUtil where + +-- Util functions used for test +-- Mostly properties and generators + +import Control.Applicative +import Control.Monad +import Control.Monad.Random +import Data.Vector ((!)) +import qualified Data.Vector as V +import Test.QuickCheck + +import Hopfield.Hopfield +import Hopfield.Measurement +import Hopfield.Boltzmann.RestrictedBoltzmannMachine +import Hopfield.Util + + +data Type = H | BM + +-- Warning: orphan instance. For more details see: +-- http://stackoverflow.com/questions/3079537/orphaned-instances-in-haskell +-- The way to avoid it is to make a warpper around Vector. +instance (Arbitrary a) => Arbitrary (V.Vector a) where + arbitrary = fmap V.fromList arbitrary + +nonempty :: forall a. Gen [a] -> Gen [a] +nonempty = (`suchThat` (not . null)) + +mapMonad :: Monad m => (a -> b) -> m [a] -> m [b] +mapMonad f m_xs = do + xs <- m_xs + return $ map f xs + + +-- | Convert a list generator to a vector generator +toGenVector :: Gen [a] -> Gen (V.Vector a) +toGenVector listGen = fmap V.fromList listGen + + +-- | Generate a random sign (+/- 1) +signGen :: Gen Int +signGen = do + n <- choose (0,1) + return $ n*2 - 1 + +binaryGen :: Gen Int +binaryGen = do + n <- choose (0,1) + return n + +-- | @patternGen n@: Generates patterns of size n +patternGen :: Type -> Int -> Gen Pattern +patternGen H n = toGenVector $ vectorOf n signGen +patternGen BM n = toGenVector $ vectorOf n binaryGen + + +-- | @patternRangeGen (min, max)@ Generates patterns of size ranging +-- between min and max +patternRangeGen :: Type -> (Int, Int) -> Gen Pattern +patternRangeGen t bounds = choose bounds >>= patternGen t + + +-- | @boundedListGen g n@: Generates lists (max length n) of the given Gen +boundedListGen :: Gen a -> Int -> Gen [a] +boundedListGen g n = do + len <- choose (0, n) + vectorOf len g + + +patListGen :: Type -> Int -> Int -> Gen [Pattern] +patListGen t maxPatSize maxPatListSize = do + i <- choose (1, maxPatSize) + nonempty $ boundedListGen (patternGen t i) maxPatListSize + + +-- | @patternsTupleGen g m1 m2@Generates a tuple of lists, as follows: +-- Uses patListGen to generate 1 list of patterns with length less than m2. +-- The list itself has to have length less than m1. +-- The second element of a tuple is a list of patterns which have the same size +-- as the patterns of the first list. +patternsTupleGen :: Type -> Int -> Int -> Gen ([Pattern], [Pattern]) +patternsTupleGen t m1 m2 = do + fst_list <- patListGen t m1 m2 + i <- choose (0, m2) + snd_list <- vectorOf i (patternGen t $ V.length $ head fst_list) + return $ (fst_list, snd_list) + + +-- Generate lists containing only 'x' +sameElemList :: a -> Gen [a] +sameElemList x = do + len <- arbitrary + return $ replicate len x + + +-- | Generate vectors containing the same element replicated +sameElemVector :: a -> Gen (V.Vector a) +sameElemVector = toGenVector . sameElemList + + +-- | Produces a matrix with 0's along the diagonal and 1's otherwise +allWeightsSame :: Int -> [[Double]] +allWeightsSame n + = [ [ if i==j then 0 else w | i <- [0..n-1] ] | j <- [0..n-1] ] + where w = (1 :: Int) ./. n + + +-- | @boundedReplicateGen n g@ Generates lists containing 'g' replicated. +-- The list is bounded in size by n. +boundedReplicateGen :: Int -> Gen a -> Gen [a] +boundedReplicateGen n g = liftM2 replicate (choose (0, n)) g + + +-- | Replaces the nth element in the list with 'r' +replaceAtN :: Int -> a -> [a] -> [a] +replaceAtN _ _ [] = error "index greater than list size" +replaceAtN 0 r (_:xs) = (r:xs) +replaceAtN n r (x:xs) + | n > 0 = (x:(replaceAtN (n-1) r xs)) + | otherwise = error "negative index" + + +-- | Compute crosstalk term for a pattern and a given neuron +-- @crosstalk hopfield index neuron +-- todo think if it is better to actually pass in the hopfield data +-- strucutre +-- the pattern on which we do this has to be one of the traninig patterns +-- todo error checks +-- note that this is a very basic check +-- one should try and implement the probability error thing as +-- that would give as a good idea of how to +-- scale +crosstalk :: HopfieldData -> Int -> Int -> Int +-- the cross talk term is h(xi k ) - xi k +crosstalk hs index n = computeH (weights hs) pat n - pat ! n + where pat = (patterns hs) !! index + + +-- | Used as a property to check that patterns which +-- are used to create the network are stable in respect to update +trainingPatsAreFixedPoints :: LearningType -> [Pattern] -> Gen Bool +trainingPatsAreFixedPoints method pats = + and <$> mapM checkFixedPoint [0.. length pats - 1] + where + hs = buildHopfieldData method pats + ws = weights hs + checkFixedPoint index = do + i <- arbitrary + return $ evalRand (update ws (pats !! index)) (mkStdGen i) == Nothing || (not $ checkFixed hs index) + + +-- | Trains a network using @training_pats@ and then updates each +-- pattern in pats according to the weights of that network. +-- The aim is to check that the energy decreases after each update. +energyDecreasesAfterUpdate :: LearningType -> ([Pattern], [Pattern]) -> Gen Bool +energyDecreasesAfterUpdate method (training_pats, pats) + = and <$> (forM pats $ \pat -> do + i <- arbitrary + return $ evalRand (energyDecreases pat) (mkStdGen i) + ) + where + ws = weights $ buildHopfieldData method training_pats + check pat afterPat = energy ws pat >= energy ws afterPat || energy ws afterPat - energy ws pat <= 0.00000001 + energyDecreases :: (MonadRandom m) => Pattern -> m Bool + energyDecreases pat = do + maybe_pat <- update ws pat + case maybe_pat of + Nothing -> return True + Just updatedPattern -> return $ check pat updatedPattern + + +-- TODO mihaela unused? +repeatedUpdateCheck :: LearningType -> ([Pattern], [Pattern]) -> Gen Bool +repeatedUpdateCheck method (training_pats, pats) + = and <$> mapM s pats + where + ws = weights $ buildHopfieldData method training_pats + stopped pat = do + p <- converged_pattern + maybe_new_p <- update ws p + return $ maybe_new_p == Nothing + where + converged_pattern = repeatedUpdate ws pat + s pat = do + i <- arbitrary + return $ evalRand (stopped pat) (mkStdGen i) + + +boltzmannBuildGen :: Int -> Int -> Int -> Gen ([Pattern], Int) +boltzmannBuildGen m1 m2 max_hidden = do + pats <- patListGen BM m1 m2 + i <- choose (1, max_hidden) + return $ (pats, i) + + +buildBoltzmannCheck :: ([Pattern], Int) -> Gen Bool +buildBoltzmannCheck (pats, nr_h) = do + i <- arbitrary + let bd = evalRand (buildBoltzmannData' pats nr_h) (mkStdGen i) + return $ patternsB bd == pats && nr_hiddenB bd == nr_h + + +boltzmannAndPatGen :: Int -> Int -> Int -> Gen ([Pattern], Int, Pattern) +boltzmannAndPatGen m1 m2 max_hidden = do + pats_train <- patListGen BM m1 m2 + i <- choose (1, max_hidden) + pats_check <- patternGen BM (V.length $ pats_train !! 0) + return $ (pats_train, i, pats_check) + + +probabilityCheck :: ([Pattern], Int, Pattern) -> Gen Bool +probabilityCheck (pats, nr_h, pat) = do + seed <- arbitrary + let bd = evalRand (buildBoltzmannData' pats nr_h) (mkStdGen seed) + ws = weightsB bd + return $ all (\x -> c $ getActivationProbability Matching Visible ws pat x) [0 .. nr_h - 1] + where c x = x <= 1 && x >=0 + + +-- -- r should only be 0 or 1 for this test +updateNeuronCheck :: Int -> ([Pattern], Int, Pattern) -> Gen Bool +updateNeuronCheck r (pats, nr_h, pat) + | not (r == 0 || r == 1) = error "r has to be 0 or 1 for updateNeuronCheck" + | otherwise = do + i <- choose (0, nr_h -1) + seed <- arbitrary + let bd = evalRand (buildBoltzmannData' pats nr_h) (mkStdGen seed) + return $ updateNeuron' (fromIntegral r) Matching Visible (weightsB bd) pat i == (1 - r) + + +-- TODO write comment and change the name to show the restrictions +buildIntTuple :: Gen (Int, Int) +buildIntTuple = do + i <- choose (1, 100) + let min_size = ceiling $ log2 $ fromIntegral i + j <- choose (min_size + 1, min_size + 2) + return (i, j) + + +binaryCheck :: (Int, Int) -> Bool +binaryCheck (x, y) = x == refold + where + refold = sum [ b * 2^pos | b <- reverse bits | pos <- [(0:: Int)..] ] + bits = toBinary x y + + +-- Runs expressions requiring random numbers (e.g. RandomMonad) in the Gen monad +evalRandGen :: Rand StdGen a -> Gen a +evalRandGen e = do + rndInt <- arbitrary + return $ evalRand e (mkStdGen rndInt) diff --git a/src/Hopfield/Util.hs b/src/Hopfield/Util.hs new file mode 100644 --- /dev/null +++ b/src/Hopfield/Util.hs @@ -0,0 +1,335 @@ +{-# LANGUAGE ParallelListComp, ScopedTypeVariables #-} + +-- | This module uses general purpose functions which are use trought the project. +-- Should not contain any project defined data types. Needs to be kept +-- as general as possible. + +module Hopfield.Util ( + average + , (*.) + , (.*) + , (./) + , (./.) + , (/.) + , attachLabel + , attachLabels + , columnVector + , combine + , combineVectors + , compareBy + , dotProduct + , findInList + , fromDataVector + , getBinaryIndices + , getElemOccurrences + , gibbsSampling + , hammingDistance + , list2D + , log2 + , normal + , numDiffs + , prettyList + , printMList + , randomBinaryVector + , randomElem + , randomSignVector + , repeatUntilEqual + , repeatUntilEqualOrLimitExceeded + , repeatUntilNothing + , runT + , shuffle + , toArray + , toBinary + , toDouble + , toPercents + , vector2D + , unfoldrSelfM + , patternToAsciiArt +) where + + +import Control.Monad +import Control.Monad.Random (MonadRandom) +import qualified Control.Monad.Random as Random +import Data.Array.ST +import Data.List.Split (chunksOf) +import qualified Data.Random as DR +import Data.List +import qualified Data.Vector as V +import Data.Word (Word32) +import Foreign.Storable +import GHC.Arr as Arr +import qualified Numeric.Container as NC +import Numeric.Probability.Random (T, runSeed) +import System.Random (mkStdGen) + +import Hopfield.Common + + +(./.) :: (Fractional a, Integral a1, Integral a2) => a1 -> a2 -> a +x ./. y = fromIntegral x / fromIntegral y + +(./) :: (Fractional a, Integral a1) => a1 -> a -> a +x ./ y = fromIntegral x / y + +(/.) :: (Fractional a, Integral a2) => a -> a2 -> a +x /. y = x / fromIntegral y + +(*.) :: (Integral a1, Num a) => a -> a1 -> a +x *. y = x * fromIntegral y + +(.*) :: (Fractional a, Integral a1) => a1 -> a -> a +x .* y = fromIntegral x * y + + +toDouble :: (Integral a, Num b) => V.Vector a -> V.Vector b +toDouble = fmap fromIntegral + + +compareBy :: Ord b => (a -> b) -> a -> a -> Ordering +compareBy f x1 x2 = compare (f x1) (f x2) + + +getElemOccurrences :: Ord a => [a] -> [(a, Int)] +getElemOccurrences = map (\xs@(x:_) -> (x, length xs)) . group . sort + + +log2 :: Double -> Double +log2 = logBase 2.0 + + + +-- | Generates a number sampled from a random distribution, given the mean and +-- standard deviation. +normal :: forall m . MonadRandom m => Double -> Double -> m Double +normal m std = do + r <- DR.runRVar (DR.normal m std) (Random.getRandom :: MonadRandom m => m Word32) + return r + + +-- | @gibbsSampling a@ Gives the binary value of a neuron (0 or 1) from the +-- activation sum +gibbsSampling :: MonadRandom m => Double -> m Int +gibbsSampling a + | (a < 0.0 || a > 1.0) = error "argument of gibbsSampling is not a probability" + | otherwise = do + r <- Random.getRandomR (0.0, 1.0) + return $ if (r < a) then 1 else 0 + + +randomElem :: MonadRandom m => [a] -> m a +randomElem [] = error "randomElem: empty list" +randomElem xs = Random.fromList (zip xs (repeat 1)) + + +repeatUntilEqual :: (MonadRandom m, Eq a) => (a -> m a) -> a -> m a +repeatUntilEqual f a = + do + new_a <- f a + if a == new_a then return a else repeatUntilEqual f new_a + +repeatUntilNothing :: (MonadRandom m) => (a -> m (Maybe a)) -> a -> m a +repeatUntilNothing f x = + do + new_x <- f x + case new_x of + Nothing -> return x + Just y -> repeatUntilNothing f y + + +repeatUntilEqualOrLimitExceeded :: (MonadRandom m, Eq a) => Int -> (a -> m a) -> a -> m a +repeatUntilEqualOrLimitExceeded limit f a + | limit < 0 = error "negative limit in repeatUntilEqualOrLimitExceeded" + | otherwise = repeatUntilEqualOrLimitExceeded' 0 limit f a + + +repeatUntilEqualOrLimitExceeded' :: (MonadRandom m, Eq a) => Int -> Int -> (a -> m a) -> a -> m a +repeatUntilEqualOrLimitExceeded' current limit f a + | current == limit = return a + | otherwise = do + new_a <- f a + if a == new_a then return a else repeatUntilEqualOrLimitExceeded' (current + 1) limit f new_a + + +-- | Converts a list of lists to a 2D vector +vector2D :: [[a]] -> V.Vector (V.Vector a) +vector2D ll = V.fromList $ map V.fromList ll + + +-- | Converts a 2D vector into a list of lists +list2D :: V.Vector (V.Vector a) -> [[a]] +list2D vv = map V.toList $ V.toList vv + +-- Returns the coumn vector of a matrix +-- Caller needs to ensure that the matrix is well formed +columnVector :: V.Vector (V.Vector a) -> Int -> V.Vector a +columnVector m i = V.map (V.! i) m + + +-- from Data.Vector to Numeric.Container.Vector +fromDataVector:: (Foreign.Storable.Storable a) => V.Vector a -> NC.Vector a +fromDataVector v = NC.fromList $ V.toList v + +-- the caller has to ensure that the dimensions are the same +combine :: (a-> b -> c) -> [[a]] -> [[b]] -> [[c]] +combine f xs ys + | length xs /= length ys = error "list sizes do not match in Utils.combine" + | otherwise = zipWith (zipWith f) xs ys + + +combineVectors :: (a -> b -> c) -> V.Vector a -> V.Vector b -> V.Vector c +combineVectors f v_a v_b + | V.length v_a /= V.length v_b = error "vector sizes do not match in dot product" + | otherwise = V.fromList (zipWith f (V.toList v_a) (V.toList v_b) ) + + +-- assertion same size and move to Util +dotProduct :: Num a => V.Vector a -> V.Vector a -> a +dotProduct xs ys + | V.length xs /= V.length ys = error "vector sizes do not match in dot product" + | otherwise = sum [ xs V.! i * (ys V.! i ) | i <- [0.. V.length xs - 1]] + + +-- Tries to find a element in a list. In case of success, returns the index +-- of the element (the first one, in case of multiple occurences). In case of +-- failure, returns the search element itself. +findInList :: Eq a => [a] -> a -> Either a Int +findInList xs x = + case m_index of + Nothing -> Left x + Just i -> Right i + where m_index = x `elemIndex` xs + + +-- @toBinary n size@. Returns the binary representation of n in size bits. +-- The caller has to ensure that n fits in size bits, or an error will be raised. +toBinary :: Int -> Int -> [Int] +toBinary n size + | n < 0 = error "toBinary requires positive arguments" + | n > 2 ^ size - 1 = error "cannot fit binary representation into given size" + | otherwise = [ (n `div` 2 ^ i) `mod` 2 | i <- [size - 1, size - 2 .. 0] ] + +-- returns the binary represenation of the indices of the elements in a list +-- after the duplicates have been removed +getBinaryIndices :: Eq a => [a] -> [(a, [Int])] +getBinaryIndices xs = [ (x, toBinary i bitsNeeded) | i <- [0 ..] | x <- nub_xs] + where + nub_xs = nub xs + bitsNeeded = 1 + (floor $ log2 $ fromIntegral (length nub_xs)) :: Int + + +-- Counts the number of pairwise differences in two lists +numDiffs :: (Eq a) => [a] -> [a] -> Int +numDiffs xs ys = length $ filter id $ zipWith (/=) xs ys + + +-- Convert list to Array +toArray :: [a] -> Array Int a +toArray xs = listArray (0, l-1) xs + where l = length xs + + +-- Efficient O(n) random shuffle of an array +-- Modified from http://www.haskell.org/haskellwiki/Random_shuffle +shuffle :: MonadRandom m => Array Int a -> m [a] +shuffle xs = do + let len = Arr.numElements xs + rands <- take len `liftM` Random.getRandomRs (0, len-1) + let shuffledArray = runSTArray $ do + ar <- Arr.thawSTArray xs + forM_ (zip [0..(len-1)] rands) $ \(i, j) -> do + vi <- Arr.readSTArray ar i + vj <- Arr.readSTArray ar j + Arr.writeSTArray ar j vi + Arr.writeSTArray ar i vj + return ar + return (elems shuffledArray) + + +-- Run a random generator T (Numeric.Probability.Random) in MonadRandom +runT :: forall m a . MonadRandom m => T a -> m a +runT dist = do + rndInt <- Random.getRandom + return $ runSeed (mkStdGen rndInt) dist + + +randomBinaryVector :: MonadRandom m => Int -> m (V.Vector Int) +randomBinaryVector size = liftM V.fromList $ replicateM size $ Random.getRandomR (0, 1) + + +randomSignVector :: MonadRandom m => Int -> m (V.Vector Int) +randomSignVector size = do + binaryVec <- randomBinaryVector size + return $ V.map (\x -> 2 * x - 1) binaryVec + + +-- Returns the average of the elements in a list +average :: (Real a, Fractional b) => [a] -> b +average xs = realToFrac (sum xs) / genericLength xs + + +hammingDistance :: V.Vector Int -> V.Vector Int -> Int +hammingDistance p1 p2 = length $ filter (== -1) $ zipWith (*) l1 l2 + where [l1, l2] = map V.toList [p1, p2] + + +-- Convert lists of double to a pretty string of (rounded) percentages +-- e.g. toPercents [0.123, 0.999] = "12% 99%" +toPercents :: [Double] -> String +toPercents ns = unwords [ show (round $ n * 100.0 :: Int) ++ "%" | n <- ns] + + +-- Prints given elements separated by a tab +attachLabel :: [Showable] -> String +attachLabel xs = concat $ intersperse "\t" $ map show xs + + +-- Tabulates the two given lists as columns +attachLabels :: String -> [[Showable]] -> String +attachLabels header is + = header ++ "\n" ++ concat list + where list = [ attachLabel i ++ "\n" | i <- is ] + + +-- Format list for output +prettyList :: Show a => [a] -> String +prettyList xs = unwords $ map show xs + + +-- Prints a list of IO actions, applying a corresponding function to it +-- e.g. printMList [IO a1, IO a2] [f1, f2] +-- Outputs the equivalent of: +-- show f1 a1 ++ show f2 a2 +printMList :: (Show a) => [IO a] -> [a -> String] -> IO () +printMList [] _ = return () +printMList _ [] = error "Function list shorter than IO action list" +printMList (x:xs) (f:fs) = do + value <- x + putStrLn $ f value + printMList xs fs + + +-- | Executes the monadic action returning a maybe until 'Nothing' is returned, +-- collecting the results in a list. +-- +-- Like `unfoldr`, the initial value is not part of the result list. +unfoldrSelfM :: Monad m => (a -> m (Maybe a)) -> a -> m [a] +-- Could be the following with `unfoldrM` from monad-loops: +-- unfoldrSelfM f seed = unfoldrM (\x -> ((\z -> (z,z)) <$>) `liftM` f x) seed +-- but monad-loops < 0.4.2 has a bug: +-- https://github.com/mokus0/monad-loops/commit/7ede550ecd2df61d12f5148b86bd5f3daaf6eb24 +unfoldrSelfM f seed = go seed + where + go a = do + mx <- f a + case mx of + Nothing -> return [] + Just x -> do xs <- go x + return $ x : xs + + +patternToAsciiArt :: Int -> Pattern -> String +patternToAsciiArt width = unlines . chunksOf width . V.toList . fmap toChar + where + toChar i | i > 0 = '1' + | otherwise = ' ' diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -0,0 +1,21 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +module Main where + +import Test.Hspec + +import TestBoltzmann +import TestHopfield +import TestBinary +import TestMeasurement + +main :: IO () +main = hspec $ do + + testBoltzmannMachine + + testHopfield + + testBinary + + testMeasurement