hopfield (empty) → 0.1.0.0

raw patch · 23 files changed

+2856/−0 lines, 23 filesdep +HUnitdep +JuicyPixelsdep +MonadRandomsetup-changed

Dependencies added: HUnit, JuicyPixels, MonadRandom, QuickCheck, array, base, deepseq, directory, erf, exact-combinatorics, hmatrix, hopfield, hspec, monad-loops, optparse-applicative, parallel, probability, random, random-fu, rvar, split, vector

Files

• Setup.hs +2/−0
• apps/ExperimentMain.hs +34/−0
• apps/Recognize.hs +181/−0
• hopfield.cabal +154/−0
• src/Hopfield/Analysis.hs +76/−0
• src/Hopfield/Benchmark.hs +12/−0
• src/Hopfield/Boltzmann/ClassificationBoltzmannMachine.hs +312/−0
• src/Hopfield/Boltzmann/RestrictedBoltzmannMachine.hs +228/−0
• src/Hopfield/Clusters.hs +157/−0
• src/Hopfield/Common.hs +39/−0
• src/Hopfield/Experiments/ClusterExperiments.hs +85/−0
• src/Hopfield/Experiments/Experiment.hs +94/−0
• src/Hopfield/Experiments/Experiment2SuperAttractors.hs +94/−0
• src/Hopfield/Experiments/ExperimentUtil.hs +49/−0
• src/Hopfield/Experiments/SmallExperiments.hs +30/−0
• src/Hopfield/Hopfield.hs +315/−0
• src/Hopfield/Images/ConvertImage.hsc +40/−0
• src/Hopfield/Images/convertImage.c +93/−0
• src/Hopfield/Measurement.hs +116/−0
• src/Hopfield/SuperAttractors.hs +133/−0
• src/Hopfield/TestUtil.hs +256/−0
• src/Hopfield/Util.hs +335/−0
• test/Main.hs +21/−0

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -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"

@@ -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"

@@ -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 <mail@nh2.me>+category:      AI, Machine Learning+stability:     experimental+synopsis:      Hopfield Networks, Boltzmann Machines and Clusters+description:   Attractor Neural Networks for Modelling Associative Memory+               .+               Report: <https://github.com/imperialhopfield/hopfield/raw/master/report/report.pdf>+               .+               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

@@ -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)

@@ -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

@@ -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+

@@ -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]+

@@ -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)

@@ -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

@@ -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"

@@ -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)

@@ -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 ]+

@@ -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]

@@ -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)

@@ -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+

@@ -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))

@@ -0,0 +1,93 @@+#include "convertImage.h"+#include <assert.h>+#include <wand/magick_wand.h>++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);+}

@@ -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

@@ -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+

@@ -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)

@@ -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 = ' '

@@ -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