1 /* -*- mode: c; c-basic-offset: 4; -*-
2  *
3  * bifurcation-diagram.c - A HistogramImager that renders a bifurcation diagram
4  *                         in which one axis interpolates across de Jong parameters
5  *                         and the other axis shows a 1 dimensional projection of the
6  *                         image at those parameters.
7  *
8  * Fyre - rendering and interactive exploration of chaotic functions
9  * Copyright (C) 2004-2006 David Trowbridge and Micah Dowty
10  *
11  * This program is free software; you can redistribute it and/or
12  * modify it under the terms of the GNU General Public License
13  * as published by the Free Software Foundation; either version 2
14  * of the License, or (at your option) any later version.
15  *
16  * This program is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with this program; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
24  *
25  */
26 
27 #include <string.h>
28 #include <math.h>
29 #include "bifurcation-diagram.h"
30 #include "math-util.h"
31 
32 static void               bifurcation_diagram_class_init        (BifurcationDiagramClass *klass);
33 static void               bifurcation_diagram_init              (BifurcationDiagram      *self);
34 static void               bifurcation_diagram_dispose           (GObject                 *gobject);
35 
36 static void               bifurcation_diagram_init_columns      (BifurcationDiagram      *self);
37 static BifurcationColumn* bifurcation_diagram_next_column       (BifurcationDiagram      *self);
38 static void               bifurcation_diagram_get_column_params (BifurcationDiagram      *self,
39 								 BifurcationColumn       *column,
40 								 DeJongParams            *param);
41 
42 static gpointer parent_class = NULL;
43 
44 
45 /************************************************************************************/
46 /**************************************************** Initialization / Finalization */
47 /************************************************************************************/
48 
bifurcation_diagram_get_type(void)49 GType bifurcation_diagram_get_type(void) {
50     static GType dj_type = 0;
51 
52     if (!dj_type) {
53 	static const GTypeInfo dj_info = {
54 	    sizeof(BifurcationDiagramClass),
55 	    NULL, /* base_init */
56 	    NULL, /* base_finalize */
57 	    (GClassInitFunc) bifurcation_diagram_class_init,
58 	    NULL, /* class_finalize */
59 	    NULL, /* class_data */
60 	    sizeof(BifurcationDiagram),
61 	    0,
62 	    (GInstanceInitFunc) bifurcation_diagram_init,
63 	};
64 
65 	dj_type = g_type_register_static(HISTOGRAM_IMAGER_TYPE, "BifurcationDiagram", &dj_info, 0);
66     }
67 
68     return dj_type;
69 }
70 
bifurcation_diagram_class_init(BifurcationDiagramClass * klass)71 static void bifurcation_diagram_class_init(BifurcationDiagramClass *klass) {
72     GObjectClass *object_class;
73     parent_class = g_type_class_ref(G_TYPE_OBJECT);
74     object_class = (GObjectClass*) klass;
75 
76     object_class->dispose      = bifurcation_diagram_dispose;
77 }
78 
bifurcation_diagram_init(BifurcationDiagram * self)79 static void bifurcation_diagram_init(BifurcationDiagram *self) {
80 }
81 
bifurcation_diagram_dispose(GObject * gobject)82 static void bifurcation_diagram_dispose(GObject *gobject) {
83     BifurcationDiagram *self = BIFURCATION_DIAGRAM(gobject);
84 
85     if (self->columns) {
86 	g_free(self->columns);
87 	self->columns = NULL;
88     }
89 
90     if (self->interpolant) {
91 	g_object_unref(self->interpolant);
92 	self->interpolant = NULL;
93     }
94 
95     if (self->interp_data && self->interp_data_free)
96 	self->interp_data_free(self->interp_data);
97     self->interp_data = NULL;
98     self->interp_data_free = NULL;
99 
100     G_OBJECT_CLASS(parent_class)->dispose(gobject);
101 }
102 
bifurcation_diagram_new()103 BifurcationDiagram* bifurcation_diagram_new() {
104     return BIFURCATION_DIAGRAM(g_object_new(bifurcation_diagram_get_type(), NULL));
105 }
106 
107 
108 /************************************************************************************/
109 /************************************************************************* Settings */
110 /************************************************************************************/
111 
bifurcation_diagram_set_interpolator(BifurcationDiagram * self,ParameterInterpolator * interp,gpointer interp_data,GFreeFunc interp_data_free)112 void bifurcation_diagram_set_interpolator (BifurcationDiagram    *self,
113 					   ParameterInterpolator *interp,
114 					   gpointer               interp_data,
115 					   GFreeFunc              interp_data_free) {
116     /* Free the old interpolator data */
117     if (self->interp_data && self->interp_data_free)
118 	self->interp_data_free(self->interp_data);
119 
120     self->interp = interp;
121     self->interp_data = interp_data;
122     self->interp_data_free = interp_data_free;
123 
124     self->calc_dirty_flag = TRUE;
125 }
126 
bifurcation_diagram_set_linear_endpoints(BifurcationDiagram * self,DeJong * first,DeJong * second)127 void bifurcation_diagram_set_linear_endpoints (BifurcationDiagram     *self,
128 					       DeJong                 *first,
129 					       DeJong                 *second) {
130     ParameterHolderPair *pair, *oldpair;
131 
132     /* If the old interpolator was also linear, see if we can avoid this update... */
133     if (self->interp == PARAMETER_INTERPOLATOR(parameter_holder_interpolate_linear)) {
134 	oldpair = (ParameterHolderPair*) self->interp_data;
135 
136 	if ((!memcmp(&DE_JONG(oldpair->a)->param, &first->param, sizeof(first->param))) &&
137 	    (!memcmp(&DE_JONG(oldpair->b)->param, &second->param, sizeof(second->param))))
138 	    return;
139     }
140 
141     pair = g_new(ParameterHolderPair, 1);
142     pair->a = g_object_ref(first);
143     pair->b = g_object_ref(second);
144 
145     bifurcation_diagram_set_interpolator(self, PARAMETER_INTERPOLATOR(parameter_holder_interpolate_linear),
146 					 pair, (GFreeFunc) parameter_holder_pair_free);
147 }
148 
149 
150 /************************************************************************************/
151 /************************************************************************** Columns */
152 /************************************************************************************/
153 
bifurcation_diagram_init_columns(BifurcationDiagram * self)154 static void bifurcation_diagram_init_columns (BifurcationDiagram *self) {
155     int hist_width;
156     int i, tmp, j;
157 
158     /* Do we need to resize the column array? */
159     histogram_imager_get_hist_size(HISTOGRAM_IMAGER(self), &hist_width, NULL);
160     if (hist_width != self->num_columns) {
161 
162 	if (self->columns)
163 	    g_free(self->columns);
164 
165 	/* Create columns and number them */
166 	self->num_columns = hist_width;
167 	self->current_column = 0;
168 	self->columns = g_new0(BifurcationColumn, self->num_columns);
169 	for (i=0; i<self->num_columns; i++)
170 	    self->columns[i].ix = i;
171 
172 	/* Shuffle them, so we render in a seemingly-random order */
173 	for (i=self->num_columns-1; i>=0; i--) {
174 	    j = int_variate(0, i+1);
175 	    tmp = self->columns[i].ix;
176 	    self->columns[i].ix = self->columns[j].ix;
177 	    self->columns[j].ix = tmp;
178 	}
179 
180 	self->calc_dirty_flag = TRUE;
181     }
182 
183     /* Should we be resetting the columns, either due to
184      * a change in interpolants or a column array resize?
185      */
186     if (self->calc_dirty_flag || HISTOGRAM_IMAGER(self)->histogram_clear_flag) {
187 	/* Clear the histogram if it isn't already */
188 	if (!HISTOGRAM_IMAGER(self)->histogram_clear_flag)
189 	    histogram_imager_clear(HISTOGRAM_IMAGER(self));
190 
191 	for (i=0; i<self->num_columns; i++) {
192 	    /* Invalidate each column and each interpolated parameter set */
193 	    self->columns[i].point.valid = FALSE;
194 	    for (j=0; j<(sizeof(self->columns[0].interpolated)/
195 			 sizeof(self->columns[0].interpolated[0])); j++) {
196 		self->columns[i].interpolated[j].valid = FALSE;
197 	    }
198 	}
199 
200 	HISTOGRAM_IMAGER(self)->histogram_clear_flag = FALSE;
201 	self->calc_dirty_flag = FALSE;
202     }
203 }
204 
bifurcation_diagram_next_column(BifurcationDiagram * self)205 static BifurcationColumn* bifurcation_diagram_next_column (BifurcationDiagram *self) {
206     /* Get the next column, wrapping around when we hit the end */
207     BifurcationColumn *column = &self->columns[self->current_column];
208     if (++self->current_column >= self->num_columns)
209 	self->current_column = 0;
210 
211     /* Initialize this column's point if it isn't yet */
212     if (!column->point.valid) {
213 	column->point.x = uniform_variate();
214 	column->point.y = uniform_variate();
215 	column->point.valid = TRUE;
216     }
217 
218     return column;
219 }
220 
bifurcation_diagram_get_column_params(BifurcationDiagram * self,BifurcationColumn * column,DeJongParams * param)221 static void bifurcation_diagram_get_column_params (BifurcationDiagram *self,
222 						   BifurcationColumn  *column,
223 						   DeJongParams       *param) {
224     /* Get a random parameter set from the given column, creating it if necessary */
225     int interpIndex = int_variate(0, sizeof(self->columns[0].interpolated)/
226 				  sizeof(self->columns[0].interpolated[0]));
227 
228     if (!column->interpolated[interpIndex].valid) {
229 
230 	/* Create an interpolant if we don't have one yet */
231 	if (!self->interpolant)
232 	    self->interpolant = de_jong_new();
233 
234 	/* Pick a random place within the column to perform the interpolation */
235 	self->interp(PARAMETER_HOLDER(self->interpolant),
236 		     (column->ix + uniform_variate()) / (self->num_columns - 1),
237 		     self->interp_data);
238 	column->interpolated[interpIndex].param = self->interpolant->param;
239 
240 	column->interpolated[interpIndex].valid = TRUE;
241     }
242 
243     *param = column->interpolated[interpIndex].param;
244 }
245 
246 
247 /************************************************************************************/
248 /********************************************************************** Calculation */
249 /************************************************************************************/
250 
bifurcation_diagram_calculate(BifurcationDiagram * self,guint iterations_total,guint iterations_per_column)251 void bifurcation_diagram_calculate (BifurcationDiagram *self,
252 				    guint               iterations_total,
253 				    guint               iterations_per_column) {
254     DeJongParams param;
255     BifurcationColumn *column;
256     HistogramPlot plot;
257     int hist_width, hist_height;
258     double x, y, point_x, point_y;
259     int i, col_i, ix, iy;
260 
261     const float y_min = -3;
262     const float y_max = 3;
263 
264     bifurcation_diagram_init_columns(self);
265     histogram_imager_prepare_plots(HISTOGRAM_IMAGER(self), &plot);
266     histogram_imager_get_hist_size(HISTOGRAM_IMAGER(self), &hist_width, &hist_height);
267 
268     for (i=iterations_total; i;) {
269 
270 	column = bifurcation_diagram_next_column(self);
271 	bifurcation_diagram_get_column_params(self, column, &param);
272 	ix = column->ix;
273 	point_x = column->point.x;
274 	point_y = column->point.y;
275 
276 	for(col_i=iterations_per_column; i && col_i; --i, --col_i) {
277 	    /* These are the actual Peter de Jong map equations. The new point value
278 	     * gets stored into 'point', then we go on and mess with x and y before plotting.
279 	     */
280 	    x = sin(param.a * point_y) - cos(param.b * point_x);
281 	    y = sin(param.c * point_x) - cos(param.d * point_y);
282 	    point_x = x;
283 	    point_y = y;
284 
285 	    if (y >= y_min && y < y_max) {
286 		iy = (int)( (y - y_min) / (y_max - y_min) * hist_height );
287 
288 		HISTOGRAM_IMAGER_PLOT(plot, ix, iy);
289 	    }
290 	}
291 
292 	column->point.x = point_x;
293 	column->point.y = point_y;
294     }
295 
296     histogram_imager_finish_plots(HISTOGRAM_IMAGER(self), &plot);
297 }
298 
299 /* The End */
300