/*! Simulation copyright (c) 2009 Reed A. Cartwright, PhD <reed@scit.us>
 *   Contains  a binomial RNG copyright (C) 1996-2003 James Theiler, Brian Gough
 *   http://www.gnu.org/licenses/gpl-2.0.html
 *   Ported from gsl/randist/binomial_tpe.c by Dan Earle <danearle@orielsystems.com>
 */

(function($) {

$.widget('ui.redlynx', {
	histcolors: [
		"#FF0000", "#FF8000", "#FFFF00", "#80FF00", "#00FF00", "#00FF80",
		"#00FFFF", "#0080FF", "#0000FF", "#8000FF", "#FF00FF", "#FF0080",
	],
	histpos: 0,
	ylabs: [
		'0.0', '', '0.1', '', '0.2', '',
		'0.0', '', '0.1', '', '0.2', '',
		'0.0', '', '0.1', '', '0.2', '',
		'0.0', '', '0.1', '', '0.2', '',
		'0.0', '', '0.1', '', '0.2', '',
		'0.0', '', '0.1', '', '0.2', '',
		'0.0', '', '0.1', '', '0.2', '',
	],
	_setID: function(a) {
		var id = this.options.idPrefix+$.data(a);
		a.attr('id', id);
		return id;
	},
	_mkButton: function(o) {
		return $('<a href="#" class="ui-redlynx-'+o.name+' ui-state-default ui-corner-all">'
			+'<span class="ui-icon ui-icon-refresh"></span>'+o.text+'</a>').hover(
			function(event, ui){ $(this).addClass('ui-state-hover'); }, 
			function(event, ui){ $(this).removeClass('ui-state-hover'); }
		).mousedown(
			function(event, ui){ $(this).addClass('ui-state-active'); }
		).mouseup(
			function(event, ui){ $(this).removeClass('ui-state-active'); }		
		).click(o.click);
	},	
	_mkOption: function(o) {
		var $p = $('<p class="ui-helper-clearfix"/>');
		var $input = $('<input type="text" value="'+ o.value.value +'" size="6"/>');
		var $slide = $('<div/>').slider({min: 0, max:1024, step: 1, range: "min",
			orientation: "horizontal", value: Math.round(1024*o.funcX(o.value.value)),
			slide: function(event, ui) {
				var $v = o.funcY(ui.value/1024);
				o.value.value = parseFloat($v);
				$input.val($v);
			}
		});
		$input.change(function() {
			var $v = $(this).val();
			o.value.value = parseFloat($v);
			$slide.slider('value', Math.round(1024*o.funcX($v)));
		});
		$p.append('<label for="'+this._setID($input)+'" title="'+o.title+'">'
			+o.label+'</label>')
		.append($input).append($slide);
		return $p;;
	},
	_init: function() {
		var self = this;
		this.$canvus = $('<div class="canvus"/>');
		this.$warn = $(
			'<div class="ui-widget"><div class="ui-state-error ui-corner-all">'
			+  '<p><span class="ui-icon ui-icon-alert"></span>'
			+  '<strong>Error:</strong> Your browser does not support SVG.  Try Firefox.</p>'
			+ '</div></div>'
		);
		this.element.addClass('ui-redlynx').append(this.$canvus).append(this.$warn);
		this.$canvus.svg({initPath: this.options.initPath, onLoad: function(zvg) {
			self.$svg = zvg;
			self.$warn.hide();
			zvg.graph.noDraw().title('Frequency of Allele A');
			zvg.graph.yAxis.title('Percentage in Population').scale(-5, 105)
				.ticks(10, 5);
			zvg.graph.xAxis.title('Generation').scale(0, 2000).ticks(2000/4,2000/16);
			zvg.graph.legend.show(0);
			zvg.graph.noDraw().chartArea([0.1, 0.1, 0.95, 0.85])
				.chartType("line").redraw();
		}});
		this.$control = $('<div class="control ui-widget"/>');
		var $p, $s, $t;
		
		// Button Bar
		$p = $('<p class="ui-helper-clearfix ui-redlynx-buttons"/>');
		
		// Run Button
		$p.append(this._mkButton({name: 'run', text: 'Run Simulation', click: function(event, ui){
			event.preventDefault();
			self.element.redlynx('simulate');
		}}));

		// Clear Button
		$p.append(this._mkButton({name: 'clear', text: 'Clear Graph', click: function(event, ui){
			event.preventDefault();
			self.element.redlynx('clear_graph');
		}}));
		
		this.$control.append($p);			
				
		// Population Size
		this.$control.size = {value: 800};
		this.$control.append(this._mkOption({
			value: this.$control.size,
			label: 'Population Size:',
			title: 'Number of diploid individuals in the population',
			funcY: function(x) { return Math.round(Math.pow(10, 5*x+1));},
			funcX: function(y) { return (Math.log(y)/Math.LN10-1)/5; }
		}));

		// Initial Allele Percent
		this.$control.start = {value: 50};
		this.$control.append(this._mkOption({
			value: this.$control.start,
			label: 'Initial % of Allele A:',
			title: 'Initial percentage of allele A in the population',
			funcY: function(x) { return Math.round(x*10000)/100; },
			funcX: function(y) { return y/100; }
		}));
				
		// Selection
		this.$control.sel = {value: 1};
		this.$control.append(this._mkOption({
			value: this.$control.sel,
			label: 'Fitness of AA:',
			title: 'Fitness of genotype AA relative to the fitness of genotype aa',
			funcY: function(x) {
					var f = Math.pow((x-0.5)*2,3);
					return Math.round(Math.pow(10,f)*10000)/10000
			},
			funcX: function(y) {
				var f = Math.log(y)/Math.LN10;
				return Math.pow(f,1/3)*0.5+0.5;
			}
		}));
		
		// Dominance
		this.$control.dom = {value: 0.5};
		this.$control.append(this._mkOption({
			value: this.$control.dom,
			label: 'Dominance:',
			title: "Determines how Aa's fitness relates to AA's and aa's",
			funcY: function(x) { return Math.round((3*x-1)*10000)/10000; },
			funcX: function(y) { return (y+1)/3; }
		}));

		// Forward Mutation
		this.$control.fmut = {value: 0};
		this.$control.append(this._mkOption({
			value: this.$control.fmut,
			label: 'Forward Mutation:',
			title: 'Probability of a->A mutation',
			funcY: function(x) {
				if(x == 0) {
					return 0;
				} else if(x >= 0.8) {
					var f = 10*x-10;
					return Math.round(Math.pow(10,f)*10000)/10000;
				} else {
					var f = 10*x-10;
					var n = Math.floor(f);
					f = Math.round(Math.pow(10,f-n)*100)/100;
					return (f+'e'+n);
				}
			},
			funcX: function(y) {
				if(y == 0)
					return 0;
				return (Math.log(y)/Math.LN10+10)/10;
			}
		}));
		
		// Backward Mutation
		this.$control.bmut = {value: 0};
		this.$control.append(this._mkOption({
			value: this.$control.bmut,
			label: 'Backward Mutation:',
			title: 'Probability of A->a mutation',
			funcY: function(x) {
				if(x == 0) {
					return 0;
				} else if(x >= 0.8) {
					var f = 10*x-10;
					return Math.round(Math.pow(10,f)*10000)/10000;
				} else {
					var f = 10*x-10;
					var n = Math.floor(f);
					f = Math.round(Math.pow(10,f-n)*100)/100;
					return (f+'e'+n);
				}
			},
			funcX: function(y) {
				if(y == 0)
					return 0;
				return (Math.log(y)/Math.LN10+10)/10;
			}
		}));
		
		// Generations
		this.$control.gens = {value: 2000};
		this.$control.append(this._mkOption({
			value: this.$control.gens,
			label: 'Generations:',
			title: 'Length of a single simulation',
			funcY: function(x) { return Math.round(10000*x); },
			funcX: function(y) { return y/10000; }
		}));
				
		this.element.append(this.$control);		
	},
	simulate: function() {
		var g = Math.round(this.$control.gens.value);
		var n = Math.round(this.$control.size.value);
		var p = this.$control.start.value/100;
		var u = this.$control.fmut.value;
		var v = this.$control.bmut.value;
		var w = this.$control.sel.value;
		var h = this.$control.dom.value;
		
		var np = Math.round(p*2*n);
		
		var y = [];
		y.push(100.0*np/(2*n));
		for(var i =1;i<=g;i+=1) {
			var p = (np*(1-v)+((2*n)-np)*u)/(2*n);
			var q = 1.0-p;
			var w12 = Math.max(1-h+h*w,0);
			var ma = p*w+q*w12;
			var mb = p*w12+q;
			p = p*ma/(p*ma+q*mb);
			
			np = this.binom.get(p, 2*n);
			y.push(100.0*np/(2*n));
		}
		this.addHistory(y);
	},
	addHistory: function(hist) {
		if(this.histpos == this.histcolors.length) {
			this.histpos = 0;
			this.$svg.graph._series = [];
		}
		var col = this.histcolors[this.histpos++];
		
		var s = this.$svg.graph.noDraw().xAxis.scale();
		if(s.max != hist.length) {
			var olen = hist.length-1;
			this.$svg.graph.xAxis.scale(0,olen).ticks(olen/4,olen/16);
		}
		
		this.$svg.graph.addSeries('',hist,col,col,2).redraw();
	},
	clear_graph: function() {
		this.histpos = 0;
		this.$svg.graph._series = [];	
		this.$svg.graph.redraw();
	},
	binom : {
		SMALL_MEAN: 14,
		BINV_CUTOFF: 110,
		FAR_FROM_MEAN: 20,
		Stirling: function (y1) {
			var y2 = y1 * y1;
			var s = (13860.0 - (462.0 - (132.0 - (99.0 - 140.0
				/ y2) / y2) / y2) / y2) / y1 / 166320.0;
			return s;
		},
		get: function(p, n) {
			var ix;                       /* return value */
			var flipped = 0;
			var q, s, np;

			if (n == 0)
				return 0;

			if (p > 0.5) {
			  p = 1.0 - p;              /* work with small p */
			  flipped = 1;
			}

			q = 1 - p;
			s = p / q;
			np = n * p;

			/* Inverse cdf logic for small mean (BINV in K+S) */

			if (np < this.SMALL_MEAN) {
				var f0 = Math.pow (q, n);   /* f(x), starting with x=0 */

				while (1) {
				/* This while(1) loop will almost certainly only loop once; but
				 * if u=1 to within a few epsilons of machine precision, then it
				 * is possible for roundoff to prevent the main loop over ix to
				 * achieve its proper value.  following the ranlib implementation,
				 * we introduce a check for that situation, and when it occurs,
				 * we just try again.
				*/
					var f = f0;
					var u = Math.random();

					for (ix = 0; ix <= this.BINV_CUTOFF; ++ix) {
						if (u < f)
							break;
						u -= f;
						/* Use recursion f(x+1) = f(x)*[(n-x)/(x+1)]*[p/(1-p)] */
						f *= s * (n - ix) / (ix + 1);
					}
					if (u < f)		// If reached satisfactory end condition, break out to return value
						break;
				}
			} else {
				/* For n >= SMALL_MEAN, we invoke the BTPE algorithm */
				var k;

				var ffm = np + p;      /* ffm = n*p+p             */
				var m = Math.floor(ffm);        /* m = int floor[n*p+p]    */
				var fm = m;            /* fm = double m;          */
				var xm = fm + 0.5;     /* xm = half integer mean (tip of triangle)  */
				var npq = np * q;      /* npq = n*p*q            */

				var p1 = Math.floor (2.195 * Math.sqrt (npq) - 4.6 * q) + 0.5;

				/* xl, xr: left and right edges of triangle */
				var xl = xm - p1;
				var xr = xm + p1;

				/* Parameter of exponential tails */

				var c = 0.134 + 20.5 / (15.3 + fm);
				var p2 = p1 * (1.0 + c + c);

				var al = (ffm - xl) / (ffm - xl * p);
				var lambda_l = al * (1.0 + 0.5 * al);
				var ar = (xr - ffm) / (xr * q);
				var lambda_r = ar * (1.0 + 0.5 * ar);
				var p3 = p2 + c / lambda_l;
				var p4 = p3 + c / lambda_r;

				var zvar, accept;
				var u, v;              /* random variates */

				while (1) { // Keep trying again.. 
					/* generate random variates, u specifies which region: Tri, Par, Tail */
					u = Math.random() * p4;
					v = Math.random();

					if (u <= p1) {
						/* Triangular region */
						ix = Math.floor (xm - p1 * v + u);
						break;
					} else if (u <= p2) {
						/* Parallelogram region */
						var x = xl + (u - p1) / c;
						v = v * c + 1.0 - Math.abs(x - xm) / p1;
						if (v > 1.0 || v <= 0.0)
							continue;
						ix = Math.floor(x);
					} else if (u <= p3) {
						/* Left tail */
						ix = Math.floor(xl + Math.log (v) / lambda_l);
						if (ix < 0)
							continue;
						v *= ((u - p2) * lambda_l);
					} else {
						/* Right tail */
						ix = Math.floor(xr - Math.log (v) / lambda_r);
						if (ix > n)
							continue;
						v *= ((u - p3) * lambda_r);
					}
					/* More efficient determination of whether v < f(x)/f(M) */

					k = Math.abs (ix - m);
					if (k <= this.FAR_FROM_MEAN) {
						/*
						* If ix near m (ie, |ix-m|<FAR_FROM_MEAN), then do
						* explicit evaluation using recursion relation for f(x)
						*/
						var g = (n + 1) * s;
						var f = 1.0;

						zvar = v;

						if (m < ix) {
							var i;
							for (i = m + 1; i <= ix; i++) {
								f *= (g / i - s);
							}
						} else if (m > ix) {
							var i;
							for (i = ix + 1; i <= m; i++) {
								f /= (g / i - s);
							}
						}
						accept = f;
					} else {
						/* If ix is far from the mean m: k=ABS(ix-m) large */

						zvar = Math.log (v);

						if (k < npq / 2 - 1) {
							/* "Squeeze" using upper and lower bounds on
							* log(f(x)) The squeeze condition was derived
							* under the condition k < npq/2-1 */
							var amaxp =
								k / npq * ((k * (k / 3.0 + 0.625) + (1.0 / 6.0)) / npq + 0.5);
							var ynorm = -(k * k / (2.0 * npq));
							if (zvar < ynorm - amaxp)
								break;

							if (zvar > ynorm + amaxp)
								continue;
						}

						/* Now, again: do the test log(v) vs. log f(x)/f(M) */

						{
						var x1 = ix + 1.0;
						var w1 = n - ix + 1.0;
						var f1 = fm + 1.0;
						var z1 = n + 1.0 - fm;

						accept = xm * Math.log (f1 / x1) + (n - m + 0.5) * Math.log (z1 / w1)
							+ (ix - m) * Math.log (w1 * p / (x1 * q))
							+ this.Stirling (f1) + this.Stirling (z1)
							- this.Stirling (x1) - this.Stirling (w1);
						}
					}

					if (zvar <= accept) {
						break;
					}
				}
			}
			if (flipped) return Math.floor (n - ix); else return Math.floor (ix);
		}
	}
});

$.ui.redlynx.defaults = {
	idPrefix: 'ui-redlynx-',
	initPath: ''
};

})(jQuery);