Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Activity
  • Sep 18 17:37
    lehins closed #106
  • Sep 18 17:37
    lehins commented #106
  • Aug 01 00:11
    lehins closed #97
  • Aug 01 00:11
    lehins closed #114
  • Aug 01 00:11
    lehins closed #115
  • Jul 31 23:29
    lehins synchronize #115
  • Jul 31 23:26
    lehins synchronize #115
  • Jul 31 22:25
    lehins synchronize #115
  • Jul 31 20:15
    lehins synchronize #115
  • Jul 31 16:10
    lehins synchronize #115
  • Jul 30 13:26
    lehins synchronize #115
  • Jul 30 12:03
    lehins synchronize #115
  • Jul 30 00:44
    lehins synchronize #115
  • Jul 29 23:36
    lehins synchronize #115
  • Jul 29 21:48
    lehins synchronize #115
  • Jul 29 20:48
    lehins synchronize #115
  • Jul 29 15:17
    lehins synchronize #115
  • Jul 29 10:42
    lehins synchronize #115
  • Jul 29 00:09
    lehins synchronize #115
  • Jul 28 13:32
    lehins synchronize #115
idontgetoutmuch
@idontgetoutmuch

I think we should try and switch random-fu to use random directly instead. I might play around with this idea and see where I can get with it.

Exactly what I was thinking - I think we could get rid of random-source

idontgetoutmuch
@idontgetoutmuch
This seems to work but is there a better way?
idontgetoutmuch
@idontgetoutmuch
simulatedDataCreate :: forall r m . (Mutable r Ix2 Double, PrimMonad m, MonadRandom m) =>
                       Double -> Double -> Int -> m (Array r Ix2 Double)
simulatedDataCreate g deltaT n = createArrayS_ (Sz2 n 2) l
  where
    l marr = do
      write marr (Ix2 0 0) 0.01
      write marr (Ix2 0 1) 0.00
      Prelude.mapM_ (f marr) [0 .. n - 1]
    f :: MArray (PrimState m) r Ix2 Double -> Int -> m Bool
    f marr i = do
      x1Prev <- read marr (Ix2 i 0)
      case x1Prev of
        Nothing -> error "Whatever"
        Just y1Prev -> do
          x2Prev <- read marr (Ix2 i 1)
          case x2Prev of
            Nothing -> error "Whatever"
            Just y2Prev -> do
              let x1New = y1Prev + y2Prev * deltaT
                  x2New = y2Prev - g * sin y1Prev * deltaT
              eta <- sample (R.Normal 0.0 0.001)
              write marr (Ix2 (i + 1) 0) (x1New + eta)
              write marr (Ix2 (i + 1) 1) x2New
Alexey Kuleshevich
@lehins
@idontgetoutmuch I thing the same way as I wrote before, just adopted to mutable array:
simulatedDataPrim :: PrimMonad m => Double -> Double -> Ix1 -> m (Array P Ix2 Double)
simulatedDataPrim g deltaT bigT =
  createArrayS_ (Sz (bigT :. 2)) $ \ma -> foldlM_ (f ma) (0.01, 0.00) (0 ..: bigT)
  where
    f ma a@(x1Prev, x2Prev) i = do
      let x1New = x1Prev + x2Prev * deltaT
          x2New = x2Prev - g * sin x1Prev * deltaT
      write_ ma (i :. 0) x1New
      write_ ma (i :. 1) x2New
      show a `trace` pure (x1New, x2New)
FYI, I am half way done of migrating from random-sourse to random for random-fu :wink: I might have a draft PR either later on tonight or tomorrow night
idontgetoutmuch
@idontgetoutmuch
Let me know if I can help
idontgetoutmuch
@idontgetoutmuch
random-fu depends on rvar which depends on random-source also
Also I think the dependency of rvar on MonadPrompt needs looking at
You are far braver than I - I have been thinking about doing what you are doing for some time
MonadPrompt doesn't even have a github repo
Alexey Kuleshevich
@lehins
@idontgetoutmuch lo end behold haskell-numerics/random-fu#67
Alexey Kuleshevich
@lehins

Let me know if I can help

If you can finish up the rough edges, that will be great:

  • there are some numeric specific things that I left out on purpose (tagged with FIXME/TODO)
  • Investigate the performance. I think a few regressions are either missing some inline pragmas or are due to the fact that floating point generation is better in random but slower than in the replaced approach of random-fu
idontgetoutmuch
@idontgetoutmuch
Running the performance tests now...
idontgetoutmuch
@idontgetoutmuch
@lehins I am not seeing any performance regressions so far - can you post what you found? thanks
idontgetoutmuch
@idontgetoutmuch
Compare1.1Vs1.2E.png
Compare1.1Vs1.2D.png
Compare1.1Vs1.2B.png
Compare1.1Vs1.2A.png
Compare1.1Vs1.2B.png
idontgetoutmuch
@idontgetoutmuch

stratified resampling is where the weights are added, divided into equal pieces and a uniform is sampled on each piece, which means that points with large weights are sampled at least once and those with small weights at most once - see here for more info: https://xianblog.wordpress.com/2012/03/13/resampling-and-gpu-parallelism/

Here's a julia implementation (which itself is based off a matlab implementation)

function resample_stratified( weights )

    N = length(weights)
    # make N subdivisions, and chose a random position within each one
    positions =  (rand(N) + collect(range(0, N - 1, length = N))) / N

    indexes = zeros(Int64, N)
    cumulative_sum = cumsum(weights)
    i, j = 1, 1
    while i <= N
        if positions[i] < cumulative_sum[j]
            indexes[i] = j
            i += 1
        else
            j += 1
        end
    end
    return indexes
end

and here is my massiv translation

resample_stratified :: forall g m r . (StatefulGen g m, PrimMonad m, Manifest r Ix1 Double) =>
                       Array r Ix1 Double -> g -> m (Array P Int Int)
resample_stratified weights stateGen = indices stateGen

  where

    bigN = elemsCount weights

    cumulative_sum :: m (Array P Int Double)
    cumulative_sum = createArrayS_ (Sz bigN) $
      (\ma -> foldM_ (f ma) 0.0 (0 ..: bigN))
      where
        f ma s i = do
          let v = weights!i
              t = s + v
          _ <- write ma i t
          return t

    -- Make N subdivisions, and chose a random position within each one
    positions :: g -> m (Array P Int Double)
    positions stateGen = createArrayS_ (Sz bigN) $
      \ma -> foldlM_ (f ma) 0.0 (0 ..: bigN)
      where
        f ma _ i = do
          epsilon <- sampleFrom stateGen (Uniform 0.0 1.0)
          let t = (epsilon + fromIntegral i) / (fromIntegral bigN)
          _ <- write ma i t
          return t

    indices stateGen = do
      ps <- positions stateGen
      cs <- cumulative_sum
      let f ma s i = do
            let go j =
                  if (ps!i) < (cs!j)
                  then do
                    _ <- write ma i j
                    return j
                  else go (j + 1)
            go s
      createArrayS_ (Sz bigN) $ \ma -> foldlM_ (f ma) 0 (0 ..: bigN)

Both in Julia and massiv, it is not clear what is going on so probably I should re-write it to be as clear as possible in Haskell without massiv before trying to optimise things.

Thoughts?

Actually this picture is quite helpful in explaining what the algorithm does:
image.png
idontgetoutmuch
@idontgetoutmuch
(b) is stratified resampling
idontgetoutmuch
@idontgetoutmuch
RS.runSTGen_ (mkStdGen 44) (resample_stratified ((fromList Seq ([0.5] ++ P.replicate 5 0.1)) :: Array B Ix1 Double)) :: Array P Int Int
I will at some point create a usable github repo for this
idontgetoutmuch
@idontgetoutmuch

This works fine

pfOneStep' :: (PrimMonad m, MonadThrow m, StatefulGen g m, MonadReader g m) =>
             Array D Ix2 Double -> Array P Int Double -> m (Array DL Ix2 Double)
pfOneStep' x_pf wn = do
  let bigN = elemsCount wn
  is <- resample_stratified' wn

  let y_pf = backpermute' (Sz (2 :. bigN)) (\(i :. j) -> i :. (is!j)) x_pf
      y1Prev = y_pf !> 0
      y2Prev = y_pf !> 1
      y1New  = y1Prev !+! y2Prev .* deltaT
      y2New  = y2Prev !-! M.map (\x -> g * sin x * deltaT) y1Prev
  stackSlicesM (Dim 2) [y1New, y2New]

but now I want to add some gaussian noise

pfOneStep :: (PrimMonad m, MonadThrow m, StatefulGen g m, MonadReader g m) =>
             Array D Ix2 Double -> Array P Int Double -> m (Array DL Ix2 Double)
pfOneStep x_pf wn = do
  let bigN = elemsCount wn
  is <- resample_stratified' wn
  eta :: Array DS Ix1 (LA.Vector Double) <- M.sreplicateM (Sz bigN) $ sample (Normal (LA.vector [0.0, 0.0]) bigQH)

  let eta1 = M.map (LA.! 0) eta

  let y_pf = backpermute' (Sz (2 :. bigN)) (\(i :. j) -> i :. (is!j)) x_pf
      y1Prev = y_pf !> 0
      y2Prev = y_pf !> 1
      y1New  = y1Prev !+! y2Prev .* deltaT !+! eta1
      y2New  = y2Prev !-! M.map (\x -> g * sin x * deltaT) y1Prev
  stackSlicesM (Dim 2) [y1New, y2New]

But I can't do this Could not deduce (Source DS Ix1 (LA.Vector Double)) arising from a use of ‘M.map’

What's my best approach here?
I can create e.g.
sampleNormal :: forall g m . (PrimMonad m, StatefulGen g m, MonadReader g m) =>
                Array P Ix1 Double -> Array P Ix2 Double -> m (Array P Ix1 Double)
sampleNormal mu sigma = ...
But I am guessing I can't add a P array to a DS array
idontgetoutmuch
@idontgetoutmuch
This is not going to be efficient but I can re-write multivariate normal to only use massiv
sampleNormal :: forall g m . (PrimMonad m, StatefulGen g m, MonadReader g m) =>
                Array P Ix1 Double -> Array P Ix2 Double -> m (Array P Ix1 Double)
sampleNormal mu sigma = do
  xs <- sample (Normal mu' sigma')
  return $ fromList Seq $ LA.toList xs
  where
    mu' = LA.fromList $ toList mu
    sigma' = LA.sym $ LA.fromLists $ toLists sigma
idontgetoutmuch
@idontgetoutmuch
Ignore my messages above - I now have something working which I will make available for comment soonish
idontgetoutmuch
@idontgetoutmuch
In Julia I can say x_pfs[:, :, :, k] = particles to set a slice of an array to that of another (lower rank) array - what's the best approach in massiv?
Alexey Kuleshevich
@lehins
@idontgetoutmuch this particular one would be: particles = x_pfs <! k. I haven't designed some cool syntax for arbitrary slice picking yet, but it can all be achieved by some combination of !>, <!>, <! and extract'
Alexey Kuleshevich
@lehins
Oh, I guess it is the opposite. x_pfs = particles <! k. However it is not exactly the same, as you will definitely notice. x_pfs[:, :, :, k] is mutation of an existing array, which I haven't implemented anything like that yet. However a writing a manual loop that iterates over the desired indices and assignes values to the mutable array at those indices is of course possible.
idontgetoutmuch
@idontgetoutmuch
Prelude.mapM_ (\(i, j, k) -> write ma (i :> j :> k :. l) (particles ! (i :> j :. k))) [(i, j, k) | i <- [1 .. 2], j <- [0 .. bigN - 1], k <- [0 .. bigT - 1]]
Alexey Kuleshevich
@lehins
@idontgetoutmuch I'd avoid lists. It is never a guarantee that they will fuse away. Delayed massiv arrays are always converted into a tight loop, so something like that: mapM_ (\ix@(i :> j :. k) -> write ma (i :> j :> k :. l) (particles ! ix)) ((1 :> 0 :. 0) ..: (3 :> bigN :. bigT))
λ> bigN = 5
λ> bigT = 6
λ> ((1 :> 0 :. 0) ..: (3 :> bigN :. bigT))
Array D Seq (Sz (2 :> 5 :. 6))
  [ [ [ 1 :> 0 :. 0, 1 :> 0 :. 1, 1 :> 0 :. 2, 1 :> 0 :. 3, 1 :> 0 :. 4, 1 :> 0 :. 5 ]
    , [ 1 :> 1 :. 0, 1 :> 1 :. 1, 1 :> 1 :. 2, 1 :> 1 :. 3, 1 :> 1 :. 4, 1 :> 1 :. 5 ]
    , [ 1 :> 2 :. 0, 1 :> 2 :. 1, 1 :> 2 :. 2, 1 :> 2 :. 3, 1 :> 2 :. 4, 1 :> 2 :. 5 ]
    , [ 1 :> 3 :. 0, 1 :> 3 :. 1, 1 :> 3 :. 2, 1 :> 3 :. 3, 1 :> 3 :. 4, 1 :> 3 :. 5 ]
    , [ 1 :> 4 :. 0, 1 :> 4 :. 1, 1 :> 4 :. 2, 1 :> 4 :. 3, 1 :> 4 :. 4, 1 :> 4 :. 5 ]
    ]
  , [ [ 2 :> 0 :. 0, 2 :> 0 :. 1, 2 :> 0 :. 2, 2 :> 0 :. 3, 2 :> 0 :. 4, 2 :> 0 :. 5 ]
    , [ 2 :> 1 :. 0, 2 :> 1 :. 1, 2 :> 1 :. 2, 2 :> 1 :. 3, 2 :> 1 :. 4, 2 :> 1 :. 5 ]
    , [ 2 :> 2 :. 0, 2 :> 2 :. 1, 2 :> 2 :. 2, 2 :> 2 :. 3, 2 :> 2 :. 4, 2 :> 2 :. 5 ]
    , [ 2 :> 3 :. 0, 2 :> 3 :. 1, 2 :> 3 :. 2, 2 :> 3 :. 3, 2 :> 3 :. 4, 2 :> 3 :. 5 ]
    , [ 2 :> 4 :. 0, 2 :> 4 :. 1, 2 :> 4 :. 2, 2 :> 4 :. 3, 2 :> 4 :. 4, 2 :> 4 :. 5 ]
    ]
  ]
λ>
Alexey Kuleshevich
@lehins
Or this slightly shorter version: mapM_ (\ix -> write ma (snocDim ix l) (particles ! ix)) (1 :> 0 :. 0 ..: 3 :> bigN :. bigT)
Alexey Kuleshevich
@lehins
FYI, using ... instead would look closer to what you had: mapM_ (\(i :> j :. k) -> write ma (i :> j :> k :. l) (particles ! (i :> j :. k))) (1 :> 0 :. 0 ... 2 :> (bigN - 1) :. (bigT - 1))
One way or another your approach with list comprehension should work just as well :+1:
idontgetoutmuch
@idontgetoutmuch
Cool - I didn't know about ...
or ..: for that matter
Callan McGill
@Boarders
Hey, I am trying to do a stencil computation which starts out with a matrix with the first row and column filled in and then propagates information through the rest of the matrix by reading the cells above, to the left and above and to the left. Is that something the stencil stuff in massiv is set up to be able to do?
image.png
visually the computation look like this
Alexey Kuleshevich
@lehins
@Boarders I've had plans to implement this functionality since the beginning of massiv. Here is a ticket that provides a simple solution to this problem, which you can use with copy+paste and a bit adjustment: lehins/massiv#88 It is a bit complex to implement it safely as a general stencil solution, because it relies on mutation, but I am confident that it can be done.
Callan McGill
@Boarders
Thanks @lehins, that is perfect. I am building out some bioinformatics tooling built on top of massiv. I'll let you know how it goes ;)
Alexey Kuleshevich
@lehins
Nice! Looking forward to hearing from you about it!