-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathex35_mandelbrot_julia.clj
182 lines (162 loc) · 6.01 KB
/
ex35_mandelbrot_julia.clj
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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
;; Manderlbrot and Julia explorer
;;
;; Complex number can be described by 2x2 matrix with values set to:
;;
;; ```
;; (a+bi) = [ a -b
;; b a]
;; ```
;;
;; With normal matrix multiplication and addition as complex operations, absolute value is square root of determinant.
;;
;; Here I test various matrix representations (not only complex) to draw mandelbrots and julias
;;
;; Press:
;;
;; * SPACE to select random matrix representation
;; * 'm' for normal Mandelbrot
;;
;; Move mouse to see Julia related to mouse position
(ns ex35-mandelbrot-julia
(:require [clojure2d.core :as c2d]
[fastmath.core :as m]
[fastmath.random :as r]
[fastmath.vector :as v])
(:import [fastmath.vector Vec4]))
(set! *warn-on-reflection* true)
(set! *unchecked-math* :warn-on-boxed)
(m/use-primitive-operators)
(def ^:const w 800)
(def ^:const hw (/ w 2))
;; First define basic matrix operations and apply them to Vec4 type
(defprotocol Matrix
(det [m])
(inverse [m])
(mulm [m1 m2]))
(extend Vec4
Matrix
{:det (fn [^Vec4 m]
(- (* (.x m) (.w m)) (* (.y m) (.z m))))
:mulm (fn [^Vec4 m1 ^Vec4 m2]
(Vec4. (+ (* (.x m1) (.x m2))
(* (.y m1) (.z m2)))
(+ (* (.x m1) (.y m2))
(* (.y m1) (.w m2)))
(+ (* (.z m1) (.x m2))
(* (.w m1) (.z m2)))
(+ (* (.z m1) (.y m2))
(* (.w m1) (.w m2)))))
:inverse (fn [^Vec4 m]
(v/mult (Vec4. (.w m)
(- (.y m))
(- (.z m))
(.x m)) (/ ^double (det m))))})
(defn absm
"Return absolute value of determinant."
^double [m]
(m/abs (det m)))
(defn draw-mandelbrot
"Draw Mandelbrot type"
[canvas make-matrix]
(dotimes [x w]
(dotimes [y w]
(let [xx (m/mnorm x 0.0 w -3.0 3.0)
yy (m/mnorm y 0.0 w -3.0 3.0)
sz (make-matrix 0.0 0.0)
c (make-matrix xx yy)
^int idx (loop [iter (int 0)
z sz]
(if (and (< iter 512)
(< (absm z) 4.0))
(recur (inc iter) (v/add (mulm z z) c))
(dec iter)))
col (* 255.0 (m/pow (/ idx 512.0) 0.4))]
(c2d/set-color canvas col col col)
(c2d/rect canvas x y 1 1))))
make-matrix)
(defn draw-julia
"Draw Julia type"
[canvas make-matrix c]
(dotimes [x hw]
(dotimes [y hw]
(let [xx (m/mnorm x 0.0 hw -3.0 3.0)
yy (m/mnorm y 0.0 hw -3.0 3.0)
sz (make-matrix xx yy)
^int idx (loop [iter (int 0)
z sz]
(if (and (< iter 64)
(< (absm z) 4.0))
(recur (inc iter) (v/add (mulm z z) c))
(dec iter)))
col (* 255.0 (m/pow (/ idx 64.0) 0.4))]
(c2d/set-color canvas col col col)
(c2d/rect canvas x y 1 1))))
make-matrix)
(defn check-matrix-def
"Returns true if matrix definition is valid (ie. reversible)"
[matrix-creator]
(let [generator (map #(vector %1 %2)
(repeatedly #(r/drand -100.0 100.0))
(repeatedly #(r/drand -100.0 100.0)))
not-zerof (filter #(not (v/is-zero? %)))
matrix-m (map matrix-creator)
det-m (map det)
checker (comp not-zerof matrix-m det-m (take 100000))]
(not-every? clojure.core/zero? (transduce checker conj generator))))
(defn make-matrix-maker
"Create matrix maker. Values are result of applying some random functions or their combinations"
[]
(let [funs [(fn ^double [^double a _] a)
(fn ^double [_ ^double b] b)
m/fast+
m/fast*
clojure.core/- #(m/atan2 %1 %2) #(m/hypot %1 %2)
(fn ^double [a _] (m/sin a))
(fn ^double [_ b] (m/sin b))
(fn ^double [a b] (r/noise a b))
(fn ^double [^double a ^double b] (if (zero? b) (/ a m/EPSILON) (/ a b)))
(fn ^double [^double a _] (- a))
(fn ^double [_ ^double b] (- b))]
rand-fn (repeatedly #(rand-nth funs))
newfuns (concat funs (take 5 (map #(let [nf (first rand-fn)]
(fn [a b]
(nf (%1 a b) (%2 a b)))) rand-fn rand-fn)))
[f1 f2 f3 f4] (repeatedly 4 #(rand-nth newfuns))
make-matrix-fn (fn make-matrix
([a b]
(Vec4. (f1 a b) (f2 a b) (f3 a b) (f4 a b)))
([[a b]]
(make-matrix a b)))]
(if (check-matrix-def make-matrix-fn) ;; recur if given matrix is invalid.
make-matrix-fn
(recur))))
(defn complex-matrix
"True Complex representation"
([a ^double b]
(Vec4. a (- b) b a))
([[a b]]
(complex-matrix a b)))
;; window, context and events
(def mcanvas (c2d/canvas w w :low))
(def mwindow (c2d/show-window {:canvas mcanvas
:window-name "Mandelbrot"
:state complex-matrix}))
(def jcanvas (c2d/canvas hw hw :low))
(def jwindow (c2d/show-window jcanvas "Julia"))
(c2d/with-canvas-> mcanvas
(draw-mandelbrot (c2d/get-state mwindow)))
(defmethod c2d/key-pressed ["Mandelbrot" \space] [_ _]
(c2d/with-canvas-> mcanvas
(draw-mandelbrot (make-matrix-maker))))
(defmethod c2d/key-pressed ["Mandelbrot" \m] [_ _]
(c2d/with-canvas-> mcanvas
(draw-mandelbrot complex-matrix)))
(defmethod c2d/mouse-event ["Mandelbrot" :mouse-moved] [e matrix-maker]
(let [nx (m/mnorm (c2d/mouse-y e) 0.0 w -3.0 3.0)
ny (m/mnorm (c2d/mouse-x e) 0.0 w -3.0 3.0)]
(c2d/with-canvas-> jcanvas
(draw-julia matrix-maker (matrix-maker nx ny)))))
(defmethod c2d/key-pressed ["Mandelbrot" \s] [_ s]
(c2d/save mcanvas (c2d/next-filename "results/ex35/mandelbrot" ".jpg"))
(c2d/save jcanvas (c2d/next-filename "results/ex35/julia" ".jpg"))
s)