Discussion :

[vasilpapa]

Bonjour,
Passionné depuis des années 70 par le graphisme et les fractales, j'ai voulu recemant m'initier à Haskell. A l'aide des consieil obtenus ici j'ai reussi de créer des petits codes pour une vingaine des fractales les plus répendues. Les codes el les images sont visible sur http://haskellaskvasil.wordpress.com.
Les problèmes sont apparus quant j'ai voulu m'attaquer aux complexes. Pour m'instruire j'ai consulté tous ce que j'ai pu trouver sur le net. Hélas, à part deux codes pour Mandel les "index hors limite", "not in scope" et modul manquant son la monnée courante. Le point commun : ils utilisent tous le module Complex.
Voici deux petits codes VB6 dessinant Mandelbrot et Julia en 1026 x 750 points et ceci en une vingtaine de secondes et sans complexe.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
Public Function Simple_Mandel()
minX = -2.4: maxX = 2.4
minY = -1.5: maxY = 1.5
w = Me.ScaleWidth: h = Me.ScaleHeight
        DoEvents
        If Touche <> 0 Then Exit Function
For x = 0 To w - 1
  rc = minX + (maxX - minX) / w * x
  For y = 0 To h - 1
    ic = minY + (maxY - minY) / h * y
    rz = 1
    iz = 0
    For a = 1 To 29: Rem choix de 32 itérations
      DoEvents
      r = rz
      i = iz
      rz = r * r - i * i + rc
      iz = 2 * r * i + ic
      module = rz * rz + iz * iz
      If module > 4 Then Exit For
    Next a
    PSet (x, y) , RGB(64, 10 * a Mod 255, 20 * a Mod 255)
Next y
Next x
While Touche = 0
    DoEvents
    Wend
End Function
 
Public Function JuliaSet()
    Dim cRe As Double, cIm  As Double 
    Dim newRe As Double, newIm As Double, oldRe As Double, oldIm As Double 
    Dim zoom As Double: zoom = 1
    moveX = 0: moveY = 0 
    Dim color As Double 
    Dim maxIterations As Integer: maxIterations = 300  
    w = Me.ScaleWidth: h = Me.ScaleHeight
    cRe = -0.7
    cIm = 0.27015
    Dim x As Integer, y As Integer
    For x = 0 To w
    For y = 0 To h
        newRe = 1.5 * (x - w / 2) / (0.5 * zoom * w) + moveX
        newIm = (y - h / 2) / (0.5 * zoom * h) + moveY
        Dim i As Integer
        For i = 0 To maxIterations
            oldRe = newRe
            oldIm = newIm
            newRe = oldRe * oldRe - oldIm * oldIm + cRe
            newIm = 2 * oldRe * oldIm + cIm
            'if the point is outside the circle with radius 2: stop
            If ((newRe * newRe + newIm * newIm) > 4) Then Exit For
        Next i
        color = RGB(255 Mod (i * 30 + 1), 255 Mod (i + 1), 20 * i) '< maxIterations))
        PSet (x, y), color
    Next y, x
    Me.Refresh
End Function

Malgré mes efforts, je n'arrive pas à faire les codes correspondants en Haskell, que j'ai commencé à apprendre depuis peu.
Alors possible ou pas possible ?

[stendhal666]

Possible sans aucun doute!
Mais il faudrait comprendre quel est ce problème avec le module Complex, j'ai du mal à croire qu'il contienne des bugs qui t'empêche de programmer tes fractales.

En tout cas je plussoie ta question, je suis allé sur ton blog qui est très sympa!

[vasilpapa]

Non, ce n'est pas le module Complexe qui pose problème.
Je voudrai tout simplement m'en passer. Hélas toutes les variantes que j'ai faites , soit tournent sans fin , soit n'affichent pas l'image voulu, mais n'importe quoi. Alors je continue à chercher.

[gorgonite]

• Pourquoi t'en passer ? Y a-t-il une incompatibilité avec un autre module que tu utilises ? juste un problème de conversion vers l'entrée de ton module d'affichage ?
• Qu'as-tu essayé en pratique ?

[vasilpapa]

Il n'y a pas de problème !
Je désire faire autrement puisque c'est possible, à quoi ça sert de faire ce que les autres ont déjà fait.

[gorgonite]

Il n'y a pas de problème !
Je désire faire autrement puisque c'est possible, à quoi ça sert de faire ce que les autres ont déjà fait.

Ah ok...


Aurais-tu un exemple de ta création qui pose problème ? histoire qu'on puisse te conseiller quelque chose...

[vasilpapa]

Bah, j'ai plusieurs tentatives de différentes manières, voici un code maladroit qui au bout de 13 minutes dessine juste une ligne pour y = 0 !
Je vois déjà que nr et ni se mordent la queue avec or et oi et il manque aussi la sortie de la boucle ( ou la changement de l'index n) dans le cas où le module est > 4.
J'ai honte.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
 
import Data.Monoid
import System.IO
import System.Random
import Graphics.Gloss
 
--(-0.7)
--0.27015
w = 300
h = 300
 
main = do
        let x = 0
        let y = 0
	let re = [1..300]
	let im = [1..300]
        let reim = [(x,y) | x<- re, y <- im]
	affiche reim x y
	main
 
affiche reim x y = display  
	 (InWindow
	       "Julia Set" 	 -- window title
		(300, 300) 	 -- window size
		(10, 10)) 	 -- window position 20
         black                   -- background color
         (picture reim x y)     -- picture to display
 
picture reim x y = pictures (listOfPoints n reim x y nr ni or oi c m)
   where
     n = 2699999
     nr = 1.5*(x-w/2)/(0.5*w)
     ni = (y-h)/(0.5*h)
     or = 0
     oi = 0
     c = 1
     m = 0
 
listOfPoints 0 _ _ _ _ _ _ _ _ _ = []
listOfPoints n reim x y nr ni or oi c m = point : listOfPoints (n-1) reim' x' y' nr' ni' or' oi' c' m' 
    where
	reim' = reim
	or' = nr
	oi' = ni
	c' = comptec  c n
	x' = (fst(reim' !! (c'-1)))
	y' = (snd(reim' !! (c'-1)))
	nr' = calcnr n x' w or' oi'
	ni' = calcni n y' h or' oi'
	m' = (n`mod`16)
	point = choixp nr' ni' x' y' m'
 
comptec c n
	|n`mod`30 <=0 = c+1
	|n`mod`30 > 0 = c 
 
calcnr n x' w or' oi'
	|n`mod`9000 == 0 = 1.5*(x'-w/2)/(0.5*w)  
	|otherwise = or'*or' - oi'*oi' -0.7
 
calcni n y' h or' oi'
	|n`mod`30 == 0 = (y'-h)/(0.5*h) 
	|otherwise = 2*or'*oi' + 0.27015
 
choixp nr' ni' x' y' m'
	|ceiling(nr'*nr'+ni'*ni') > 4 = translate (x'-150) (y'-150) (color (myColor(m')) (circle (1.0)))
	|ceiling(nr'*nr'+ni'*ni' ) <= 4 = translate (0) (0) (color black (circle (1.0)))
 
myColor m' = makeColor8 (255`mod`(m'*30+1)) (255`mod`(m'+1)) (20*m') 255

[vasilpapa]

Après pas mal de réflexions je suis en train d'aboutir.
Évidemment c'est très, très lent. Reste aussi à améliorer le code pour la couleur, afin de faire apparaître toutes les spires de Julia. Voici le résultat actuel :

le code sera visible sur http://haskellaskvasil.wordpress.com
A bientôt.

[stendhal666]

Hey Vasilpapa, en regardant les papiers de recherche importants rédigés avec Haskell, je suis tombé sur ça:
http://web.cecs.pdx.edu/~mpj/pubs/co...g-fractals.pdf
qui cible en particulier Julia et Mandelbrot, avec de bonnes idées sur la composition modulaire du programme. Du coup j'ai pensé à toi!

[vasilpapa]

LOL!
Tu pense bien que j'ai épluché le web . Hélas, rien pour Haskell a part Julia et Mandelbrot avec module Complex.
J'ai tout de même reussi, mais que c'est lent!. Un programmeur chevronné en Haskell pourrai certainement faire un code meilleur et plus rapide. Le mandelbrot de Yan Esposito et autres ne sont pa telement rapides non plus en comparaison avec VB6.
Merci de suivre mon poste.

[Jedai]

Vu que le papier proposé par Stendhal666 ne fait pas allusion aux complexes (qui sont simplement une façon de représenter les points du plan et de décrire simplement les calculs qu'on effectue dessus, il n'y a donc aucune raison de refuser de regarder des solutions juste parce qu'ils utilisent des complexes...), il serait bon que tu lises entièrement les réponses qui te sont faites.
Par ailleurs je noterais que l'une des raisons pour lesquelles ton code peut être lent est parce qu'il crée un objet OpenGL par point (puisque tu dessines des cercles pour tes points) au lieu de décider de la couleur d'un pixel d'une image.

De plus tu présentes ton code VB6 comme "sans complexe" mais il contient ce morceau :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
        For i = 0 To maxIterations
            oldRe = newRe
            oldIm = newIm
            newRe = oldRe * oldRe - oldIm * oldIm + cRe
            newIm = 2 * oldRe * oldIm + cIm
            'if the point is outside the circle with radius 2: stop
            If ((newRe * newRe + newIm * newIm) > 4) Then Exit For
        Next i

qui montre bien que l'unique raison pour laquelle il n'y a pas de complexe dans ton code est parce que VB6 (ou toi) ne sait pas les représenter. Ces calculs sont exactement calqués sur ce que tu ferais bien plus simplement avec des complexes, ici les *Re contiennent la partie réelle et *Im la partie imaginaire d'un complexe. Je ne comprend vraiment pas ton obsession de te passer des complexes.
--
Jedaï

[vasilpapa]

Je ne veux pas me passer des complexes, mais du module Complexe. D'autre part, les codes disponible sur le net qui utilisent ce module et que j'ai testé sont aussi lents, même si il sont plus rapides ( et de loin) que mon code.
La manière de dessiner les points plus rapidement m'intéresse.
Merci.

[Jedai]

Simplement au lieu de dessiner les points, tu peux générer une image pixel par pixel et l'afficher, cela sera probablement plus rapide :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
{-# LANGUAGE FlexibleContexts #-}
module Main (main) where
 
import Codec.Picture (generateImage, PixelRGBA8(..))
import Graphics.Gloss.Juicy
import Graphics.Gloss hiding (Point)
import Data.Complex
import Data.Array.Unboxed
 
main = do
  let xsteps = 500
      ysteps = 500
      start@(xs,ys) = (-2,-1.5)
      precision = 0.006
      end = (xs+ fromIntegral xsteps*precision, ys+ fromIntegral ysteps*precision)
      limit = 10
      n = 50
      gen x y = chooseColor n $ mandel ! (x,y)
      mandel = grid xsteps ysteps start end (mandelbrot limit n)
  display
    (InWindow "Fractale"
     (xsteps, ysteps) (10, 10))
    (colorFromPixel . head $ palette n)
    (fromImageRGBA8 $ generateImage gen xsteps ysteps)
 
type Point = (Double, Double)
 
mandelbrot :: Double -> Int -> (Point -> Int)
mandelbrot limit n (x,y) = length . take n . takeWhile fairlyClose . iterate (jc (x:+y)) $ 0
  where
    jc c z = sqr z + c
    fairlyClose z = magSqr z < limit
 
chooseColor :: Int -> Int -> PixelRGBA8
chooseColor n = \iter -> colors !! iter
  where
    colors = palette n
 
palette n = map (convert . (255-) . round) $ for (n+1) 0 255
  where
    convert g = PixelRGBA8 g g g 255
 
----------------------------------------------
-- Fonctions utilitaires
----------------------------------------------
 
colorFromPixel (PixelRGBA8 r g b a) = makeColorI (fI r) (fI g) (fI b) (fI a)
  where
    fI = fromIntegral
 
grid :: (IArray UArray a) => Int -> Int -> Point -> Point -> (Point -> a) -> UArray (Int,Int) a
grid r c (xi,yi) (xf,yf) f = listArray ((0,0),(r-1,c-1)) [f (x,y)
                                                         | x <- for r xi xf
                                                         , y <- for c yi yf
                                                         ]
 
for :: Int -> Double -> Double -> [Double]
for n i f = [i, i+delta .. f]
  where
    delta = (f - i) / fromIntegral (n-1)
 
magSqr :: Complex Double -> Double
magSqr (x:+y) =  scaleFloat k (sqr (scaleFloat mk x) + sqr (scaleFloat mk y))
  where
    k  = max (exponent x) (exponent y)
    mk = - k
 
sqr :: (Num a) => a -> a
sqr z = z * z

Mon usage des complexes ici est des plus minime, tu peux facilement t'en passer en modifiant "jc" et "magSqr" si tu le souhaites.
J'espère que ce sera plus rapide que ton code, pour moi avec les paramètres choisis ici et compilé avec -O2, l'affichage est pratiquement instantané.

La palette n'est pas très intéressante pour l'instant mais je suis sûr que tu sauras y remédier.
--
Jedaï

[vasilpapa]

Ah, merci, je n'avais pas vue cette réponse.
Voila le simple Mandel, hélas aussi lent.

Je vais tout revoir.

[vasilpapa]

Effectivement on peut colorier :
Je voie que seules sont définies deux régions et la frontière. On doit plutôt obtenir ça :

je vais essayer.

[Jedai]

Je voie que seules sont définies deux régions et la frontière.

Non, tu as une gradation, tu peux régler "n" et "limit" pour affecter de différente façon la précision de cette gradation, le choix de la palette peut également changer le résultat visuellement du tout au tout.

Voici une version interactive plus avancée, nettement plus rapide d'origine de façon à permettre la régénération rapide des résultats lorsque tu changes de point de vue (en attrapant avec la souris) ou que tu zoomes (roulette), tu peux aussi revenir au point de vue de départ avec "r". J'ai désactivé la rotation pour gagner un peu de temps et parce qu'elle n'apportait pas grand chose.

Cette version n'est clairement pas terriblement optimisée mais sur mon ordinateur avec les paramètres choisis et en compilant avec :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

-O2 -fexcess-precision -optc-ffast-math -optc-O3 -optc-march=x86-64 -optc-mfpmath=sse -threaded -Odph -fno-liberate-case -funfolding-use-threshold1000 -funfolding-keeness-factor1000 -with-rtsopts="-N"

J'obtiens un rendu fluide, largement utilisable pour explorer la fractale (l’auto-similarité est vraiment une propriété étonnante) jusqu'aux limites de ce que peuvent faire 100 itérations max par pixel et des Doubles pour représenter les coordonnées...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
{-# LANGUAGE FlexibleContexts, BangPatterns #-}
module Main (main) where
 
import Graphics.Gloss
import Graphics.Gloss.Raster.Field
import Graphics.Gloss.Data.ViewPort
import Graphics.Gloss.Data.ViewState
import Graphics.Gloss.Interface.Pure.Game
import Data.Complex
import GHC.Float
 
main = play displayMode background refreshs initWorld render handleEvent handleRefresh
  where
    sizeX = 500
    sizeY = 350
    limit = 4
    n = 100
    colors = palette n
    initWorld   = World (sizeX, sizeY) viewStateInit
    gen w = (colors !!) . mandelbrot limit n . rectifyView w
    displayMode = InWindow "Ensemble de Mandelbrot" (sizeX, sizeY) (10, 10)
    background = head $ palette n
    refreshs = 2
    render w@(World (sizeX,sizeY) _) = makePicture sizeX sizeY 1 1 (gen w)
    handleEvent evt w = case evt of
      EventResize size -> w {winSize = size}
      _                -> w {vs = updateViewStateWithEvent evt (vs w)}
    handleRefresh t w = w
 
data World = World {winSize :: (Int,Int), vs :: ViewState}
 
mandelbrot :: Double -> Int -> ((Double,Double) -> Int)
mandelbrot limit n (x,y) = go 0 0
  where
    go !step z
      | fairlyClose z && step < n = go (step + 1) (z*z + c)
      | otherwise                 = step
    c = x :+ y
    fairlyClose (x:+y) = x*x+y*y < limit
 
palette n = map (convert . (255-) . round) $ for (n+1) 0 255
  where
    convert g = rgb8w g g g
 
----------------------------------------------
-- Fonctions utilitaires
----------------------------------------------
 
rectifyView :: World -> Point -> (Double,Double)
rectifyView (World (sizeX, sizeY) vs) = transAndScale (viewStateViewPort vs) . rectifyAspect
  where
    rectifyAspect (x,y) = (fromIntegral sizeX / fromIntegral sizeY * float2Double x, float2Double y)
    transAndScale ViewPort {viewPortScale = zoom, viewPortTranslate = (vtx,vty)} (x,y) =
      let z = 1/float2Double zoom in
      (x*z - float2Double vtx/fromIntegral sizeX, y*z - float2Double vty/fromIntegral sizeY)
 
for :: Int -> Double -> Double -> [Double]
for n i f = [i, i+delta .. f]
  where
    delta = (f - i) / fromIntegral (n-1)

Toujours une palette monochrome qui pourrait être égayée par quelques couleurs.

--
Jedaï

[vasilpapa]

Oui, j'ai déjà joué sur limit et n mais ce n'est pas concluant. Parfois même ca fait apparaître des formes étonnantes. Finalement les meilleurs resultat sont obtenu avec le code de Yan Esposito en ajoutant un "display" et en jouant sur la définition des couleurs. Ca a donné l'image ci-dessus.

[Jedai]

Essaye donc avec ma dernière version, apparemment j'avais essayé de faire trop compliqué avec mon "magSqr", les scaleFloat produisaient des aberrations, une comparaison directe est plus efficace. Par ailleurs la limite de 4 (soit 2 pour le module) est la plus correcte mathématiquement et donne de bons résultats, 100 itérations est très suffisant et avec une bonne palette tu peux obtenir des résultats saisissants, tu peux par exemple utiliser "mixColor" :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
module Main (main) where
import Graphics.Gloss.Raster.Field
import Graphics.Gloss.Data.ViewPort
import Graphics.Gloss.Data.ViewState
import Graphics.Gloss.Interface.Pure.Game
import Data.Complex
import GHC.Float
 
data World = World {winSize :: (Int,Int), vs :: ViewState}
 
main = play displayMode background refreshs initWorld render handleEvent handleRefresh
  where
    n = 100
    palette = makePalette n
    initWorld   = World (700, 600) viewStateInit
    gen w = (palette !!) . mandelbrot n . rectifyView w
    displayMode = InWindow "Ensemble de Mandelbrot" (winSize initWorld) (10, 10)
    background = head $ palette
    refreshs = 2
    render w@(World (sizeX,sizeY) _) = makePicture sizeX sizeY 1 1 (gen w)
    handleEvent evt w = case evt of
      EventResize size -> w {winSize = size}
      _                -> w {vs = updateViewStateWithEvent evt (vs w)}
    handleRefresh t w = w
 
makePalette :: Int -> [Color]
makePalette n = white : colorsFromTo n1 white blue ++ colorsFromTo n2 blue orange ++ colorsFromTo n3 orange green ++ [black]
  where
    n' = n - 1
    (n1,n2,n3) = (n' `div` 3, n' `div` 3, n' - n1 - n2)
    colorsFromTo steps cfrom cto = map mix [1..steps]
      where mix step = let p = fromIntegral step / fromIntegral steps
                       in mixColors (1-p) p cfrom cto
 
mandelbrot :: Int -> ((Double,Double) -> Int)
mandelbrot n (x,y) = go 0 0
  where
    go step z
      | fairlyClose z && step < n = go (step + 1) (z*z + c)
      | otherwise                 = step
    c = x :+ y
    fairlyClose (x:+y) = x*x+y*y < 4
 
----------------------------------------------
-- Fonctions utilitaires
----------------------------------------------
 
rectifyView :: World -> Point -> (Double,Double)
rectifyView (World (sizeX, sizeY) vs) = transAndScale (viewStateViewPort vs) . rectifyAspect
  where
    rectifyAspect (x,y) = (fromIntegral sizeX / fromIntegral sizeY * float2Double x, float2Double y)
    transAndScale ViewPort {viewPortScale = zoom, viewPortTranslate = (vtx,vty)} (x,y) =
      let z = 1/float2Double zoom in
      (x*z - float2Double vtx/fromIntegral sizeX, y*z - float2Double vty/fromIntegral sizeY)

Cette fractale est vraiment plus intéressante lorsque tu peux l'explorer.

--
Jedaï

[vasilpapa]

Oups :
Could not find module `Graphics.Gloss.Raster.Field'
Use -v to see a list of the files searched for.

[Jedai]

Il faut installer gloss-raster : ce package utilise les tableaux repa pour définir un makePicture qui calcule automatiquement en parallèle la couleur des pixels à afficher. C'est idéal pour ce type de fractale et plus rapide. Juste un soucis c'est que makePicture et les autres fonctions de ce module demande des points de (-1,-1) à (1,1), il faut donc les modifier avant de les passer à la fonction "Point -> Color" de façon à obtenir un résultat correspondant au ViewPort actuel, c'est le rôle de ma fonction rectifyView.

Il y a aussi une petite faiblesse dans mon code actuel : je m'appuie sur ViewState pour gérer les commandes permettant de se déplacer dans la fractale et ViewState utilise (Float, Float) pour stocker la position du centre de l'écran alors que mon code calcule les coordonnées de ses points avec une finesse de (Double, Double)... Il arrive donc un moment où tu ne peux plus te déplacer parce que tu as atteint la limite des Float mais le dessin est encore net ce qui donne une impression étrange... Tu peux continuer à zoomer par contre. Pour rectifier cela, il faudrait que je gère moi-même les évènements claviers/souris (il n'y a pas de limitation fondamentale de ce côté) mais j'ai préféré la simplicité de ViewState vu qu'on peut déjà zoomer d'un facteur 1 million (?) avant que le problème apparaisse.

--
Jedaï