Mandelbrot et Julia sans le module Complex

Version imprimable

Voir 40 message(s) de cette discussion en une page

2 pièce(s) jointe(s)

Mandelbrot et Julia sans le module Complex

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:

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
Pièce jointe 173072
Pièce jointe 173073

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 ?

30/03/2015, 19h13
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!
30/03/2015, 20h19
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.
31/03/2015, 09h37
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 ?
31/03/2015, 12h28
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.
01/04/2015, 10h24
gorgonite

Citation:

Envoyé par 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.

Ah ok... :)

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

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:

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

13/04/2015, 17h55
vasilpapa

1 pièce(s) jointe(s)

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 :

Pièce jointe 174698
le code sera visible sur http://haskellaskvasil.wordpress.com
A bientôt.
14/04/2015, 13h55
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!
14/04/2015, 14h21
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.

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:

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ï

14/04/2015, 20h56
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.

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:

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ï

18/04/2015, 20h21
vasilpapa

1 pièce(s) jointe(s)

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

Pièce jointe 175260
Je vais tout revoir.
18/04/2015, 21h26
vasilpapa

2 pièce(s) jointe(s)

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

Pièce jointe 175264
je vais essayer.

Citation:

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:

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

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ï

19/04/2015, 06h45
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.

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:

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ï

19/04/2015, 13h38
vasilpapa

Oups :
Could not find module `Graphics.Gloss.Raster.Field'
Use -v to see a list of the files searched for.
19/04/2015, 14h51
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ï

Voir 40 message(s) de cette discussion en une page