model.c

// A GSList list (Glib) of interacting oscillators.
// The Kuramoto model is used for interaction between oscillators
// The Runge–Kutta method is implemented as an iterator on a GSList of oscillators



#define N 128 			// number of Kutamoto ocillators (each a cell agent here). N
#define Avg 1.99			// und
#define StDev 0.25 		// parameters for the distribution of oscillator frequencies

#define T0 0  			// initial time. t0
#define TN 16 			// final time. tn
#define H 0.075			// integration step (time differential dt). h
#define K 2.5			// interaction strength of the Kuramoto model. K

#define N_POINTS 256	// time-buffer length to store coherence valuers. r_t. 

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <complex.h>
#include <glib.h>

// pipe to gnuplot for plotting
#include "gnuplot_i.h"
// pseudo random number generation
#include "mt64.h"

// coherence and other records
double coherence[(int) N_POINTS]; 	// array to store coherence dynamics
double complex Phase[N];			// a complex number
double X[N],Y[N];

// cycle of a given oscillator
double two_pi = 2 * M_PI; // 2 * pi


// pseudo random normal disreibution from unifeorm mt64
// using the Box–Muller transform and a met64.h uniform prng.

double rand_normal()
        {
        double u1 = genrand64_real3();
        double u2 = genrand64_real3();
        return sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2);
        }


// ODE
// differential equation RH side. Declaration. See definition below
float f(float t,float theta, float omega, float psi, float r); 

// Runge-Kutta integrator for Kuramoto model. Declaration. See definition below
// it integrates the differential equation defined in f(t, theta) for frequency omega
float integrateKuramoto(float t0, float theta0, /*float tn,*/ float h, float omega, double psi, double r);

// Environ/Interactome structure
typedef struct{
		double complex z; // complex order parameter
		double r;		  // coherence
		double psi;		  // average phase
} Patch;

// Single cell/agent structure
typedef struct {
        char * name;	// id
        float  phase;  // should run from 0 to 2Pi. or 0 to 1.
		float frequency; // intrinsic frequency value
} Cell;


// constructor for a cell (oscillator)
Cell *  make_cell(char *name, float phase, float frequency)
	{
	Cell * cell = g_new(Cell, 1);
	cell->name = name;
	cell->phase = phase;
	cell->frequency = frequency;
	return cell;
	}

// constructor for a patch (environment)
Patch * make_patch(GSList *population, GSList *iterator)
	{
	Patch * habitat = g_new(Patch, 1);
	double complex coherence_number = 0;

	// get statistics over the collection (coherence)
	for (iterator = population; iterator ; iterator = iterator->next)
			coherence_number += (cexp(I*((Cell *)iterator->data)->phase))/N;
	// print environment
	printf("\t# The order parameter, a complex number, is: z = ");
	printf("%lf + %lfi\n", creal(coherence_number), cimag(coherence_number));
	// set up value state
	habitat->z = coherence_number;
	habitat->r = cabs(coherence_number);
	habitat->psi = carg(coherence_number);
	return habitat;
	};


Patch * update_patch(Patch * habitat, GSList *population, GSList *iterator){
	double complex coherence_number = 0;
	// get sum of complex exponentials over the collection of phase values
	for (iterator = population; iterator ; iterator = iterator->next)
			coherence_number += (cexp(I*((Cell *)iterator->data)->phase))/N;
	// complex order parameter captured
	printf("\t# The order parameter, a complex number, is: z = ");
	printf("%lf + %lfi\n", creal(coherence_number), cimag(coherence_number));
	// set up value state
	habitat->z = coherence_number;
	habitat->r = cabs(coherence_number);
	habitat->psi = carg(coherence_number);	
	return habitat;
	}

// print ocillator values
void print_iterator(Cell *item){
	printf("\t# %s: theta: %f  [x: %f  y: %f] omega: %f\n",
	(char *)item->name, (float)item->phase, cos((float)item->phase), sin((float)item->phase), (float)item->frequency/M_PI);
	}



// integrate ODE
void evolve_iterator(Cell *item, Patch *habitat, float t){
	item->phase = integrateKuramoto (t, item->phase, /*(float)TN,*/ (float)H, item->frequency, habitat->psi, habitat->r);
	}

// plot the phases in the population of oscillators at a given time
void plot_population(gnuplot_ctrl *h, GSList *population)
	{
	// plot around here
	for(guint i = 0;i<N;i++)
		{
			Cell * este = g_slist_nth_data (population, i);
			Phase[i]= cexp(I*(este->phase));
			X[i]= creal(Phase[i]);
			Y[i]= cimag(Phase[i]);
		}
	gnuplot_cmd(h, "set pointsize 2.5");
	gnuplot_cmd(h, "set size square");
	gnuplot_cmd(h, "set title 'Kuramoto model'");
	gnuplot_cmd(h, "set xrange [-1.5:1.5]");
	gnuplot_cmd(h, "set yrange [-1.5:1.5]");
	gnuplot_set_axislabel(h, "x", "Real");
	gnuplot_set_axislabel(h, "y", "Imaginary");
	gnuplot_cmd(h, "set object 1 circle front at 0.0,0.0 size 1.0  lw 0.25");
	gnuplot_plot_coordinates(h, X , Y, N, "Oscillator Phase");
}

// plot the dynamics of the oscillator's coherence
void plot_coherence_dynamics(gnuplot_ctrl *h)
	{
	gnuplot_set_axislabel(h, "x", "time (t)");
	gnuplot_set_axislabel(h, "y", "coherence (r)");
	gnuplot_cmd(h,"set title 'Coherence dynamics'");
	gnuplot_plot_coordinates(h, coherence, NULL, N_POINTS, "Oscillators coherence");
}

// Main
int main(int argc, char** argv)
	{
	//initialize the mt algorithm for random number genration
	unsigned long long seed = (unsigned int)time(NULL);
	init_genrand64(seed);

	// Get a handle to control a GnuPlot session (via unix pipe)
	gnuplot_ctrl *h1, *h2;
	// Initialize a gnuplot session on the handle
	h1 = gnuplot_init(); h2 = gnuplot_init();

	// Set an Initial value (condition)
	printf("\n\n Setting up initial conditions\n\n");

	// Make a population of N agents/ocillators
	GSList *population = NULL, *iterator = NULL;
	int number;  char name[20];
	for (number = 0; number < N; number ++){
	sprintf(name, "oscillator %d", number);
	population = g_slist_append(population, make_cell(g_strdup(name), 
	(float) two_pi*genrand64_real1(), Avg + StDev * rand_normal()));
	// (-1)*two_pi/8 + genrand64_real1()/10));
	} // done.

	// print the list of oscillators
	printf("\n\t# The population size is %d oscillators\n", g_slist_length(population));
	g_slist_foreach(population, (GFunc)print_iterator, NULL);

	// plot both, the oscillator phases and coherence dynamics
	plot_population(h1,population);
	plot_coherence_dynamics(h2);

	// evaluate (patch environment) --- complex order parameter
	Patch * oscillators_habitat = make_patch(population, iterator);
	printf("\n\t# The environment induced by these oscillators is:\n");

	// print the env
	printf("\t# The argument 'psi' is: %f \n", oscillators_habitat->psi);
	printf("\t# The absolute value 'r' is: %f \n", oscillators_habitat->r);

	// Start simulation
	printf("\n\n Integrating ODEs for time evolution\n\n");
	int i = 0;
	float t, t0, tn, h;
	t0 = (float) T0; h = (float) H; tn= (float)TN;
	t = t0 + h;
	// evolve the oscillators (ODEs)
	while(t<tn+h)
	{	printf("\nGeneration %f\n", t);
		// move the oscillator phase
		g_slist_foreach(population, (GFunc)evolve_iterator, (gpointer *) oscillators_habitat);
		// print oscillator
		printf("\n\n\tPopulation:\n");
		g_slist_foreach(population, (GFunc)print_iterator, NULL);	
		// update (environ) order parameter
		printf("\n\n\tEnvironment:\n");
		oscillators_habitat = update_patch(oscillators_habitat, population, iterator);
		t=t+h;
		// print environ
		printf("\t# The argument 'psi' is: %f \n", oscillators_habitat->psi);
		printf("\t# The absolute value 'r' is: %f \n", oscillators_habitat->r);
		// capture coherence vaue at a given instant
		coherence[i]= oscillators_habitat->r; i++;
		// plot
		gnuplot_resetplot(h1);plot_population(h1,population);  
		gnuplot_resetplot(h2); plot_coherence_dynamics(h2);
		// pause 1 sec
		sleep(1);
	}

	// say goodbye	
	printf("\n\n\tEnd of Simulation!  In 15 secs, both (gnu)plots are out\n");
	sleep(15); 
	printf("\n\n\tEnd of Program. Thanks for running me!\n");

	// clean up and close
	g_slist_foreach(population, (GFunc)g_free, NULL);
	g_slist_free(population);
	g_free(oscillators_habitat);
	gnuplot_close(h1);gnuplot_close(h2);
	return 0;
}



// f: RHS function
// ODE to be integrated
float f(float t, float theta, float omega, float psi, float r)
	{
	// Here we implement the Kuramoto
	float m;
	m = omega + (K*r)*sin(psi-theta);
	// omega is the intrinsic frequency of the ocillator
	// theta is the current phase
	// r and psi the order parameters
	// K the interaction strength
	return m;
	}


// Runge–Kutta method to integrate f above
float integrateKuramoto(float t0, float theta0, /*float tn,*/ float h, float omega, double psi, double r)
        {
        float m1, m2, m3, m4, m, theta;
        // Set initial condition
        theta = theta0;
        // interate the simulation one single step of length h
                // Set up Runge-Kutta parameters
                m1=f(t0,theta0, omega, psi, r);
                m2=f((t0+h/2.0),(theta0+m1*h/2.0), omega, psi, r);
                m3=f((t0+h/2.0),(theta0+m2*h/2.0), omega, psi, r);
                m4=f((t0+h),(theta0+m3*h), omega, psi, r);
                // Calculate the slope/flow
                m=((m1+2*m2+2*m3+m4)/6);
            // iterate the system by updating theta
            theta = theta+m*h;
            theta = fmod(theta + two_pi, two_pi);
        	return theta;
        }