package org.singinst.uf.math;

import junit.framework.Assert;
import junit.framework.TestCase;

import org.apache.commons.math.analysis.RombergIntegrator;
import org.apache.commons.math.analysis.UnivariateRealFunction;

public class MathUtilTest extends TestCase {

	static MathUtil.EventDiscreteDistributionSchedule heavisideSchedule(double t) {
		MathUtil.EventDiscreteDistributionSchedule o = new MathUtil.EventDiscreteDistributionSchedule();
		double[] o_dot_time = {0, t*(1-MathUtil.DOUBLE_EPS), t, 1e100};
		o.time = o_dot_time;
		double[] o_dot_clogProbEvent = {0,0,10,10};
		o.clogProbEvent = o_dot_clogProbEvent;
		o.logitSubevents = new double[0][];
		return o;
	}
	
	static MathUtil.ScienceSpeedModelParameters nullScienceModel() {
		MathUtil.ScienceSpeedModelParameters o = new MathUtil.ScienceSpeedModelParameters();
	 	o.sequencedata_base_year = 0;
	 	o.sequencedata_mean_log_year = 0;
	 	o.sequencedata_stddev_log_year = 0;
	 	o.sequencedata_mean_init_log_factor = 0;
	 	o.sequencedata_stddev_init_log_factor = 0;
	 	o.sequencedata_mean_slope_log_factor = 0;
	 	o.sequencedata_stddev_slope_log_factor = 0;
	 	o.population_base_year = 0;
	 	o.population_mean_log_year = 0;
	 	o.population_stddev_log_year = 0;
	 	o.population_mean_init_log_rate = 0;
	 	o.population_stddev_init_log_rate = 0;
	 	o.population_mean_slope_log_rate = 0;
	 	o.population_stddev_slope_log_rate = 0;
	 	o.institutional_base_year = 0;
	 	o.institutional_mean_slope_log_factor = 0;
	 	o.institutional_stddev_slope_log_factor = 0;
	 	return o;
	}

	static MathUtil.ScienceSpeedZValues nullScienceSpeedZValues () {
		MathUtil.ScienceSpeedZValues o = new MathUtil.ScienceSpeedZValues();
		return o;
	}
	
	static MathUtil.ScienceSpeedModelParameters pointMassScienceSpeedModel(MathUtil.ScienceSpeedModelParameters m, MathUtil.ScienceSpeedScenario s) {
		MathUtil.ScienceSpeedModelParameters o = nullScienceModel();
		o.institutional_mean_slope_log_factor = s.institutional_slope_log_factor;
		if (s.population_year > m.population_base_year) {
			o.population_mean_log_year = Math.log(s.population_year - m.population_base_year);
			o.population_base_year = m.population_base_year;
		} else {
			o.population_mean_log_year = Double.NEGATIVE_INFINITY;
			o.population_base_year = s.population_year;
		}
		o.population_mean_init_log_rate = s.population_init_log_rate;
		o.population_mean_slope_log_rate = s.population_slope_log_rate;
		if (s.sequencedata_year > m.sequencedata_base_year) {
			o.sequencedata_mean_log_year = Math.log(s.sequencedata_year - m.sequencedata_base_year);
			o.sequencedata_base_year = m.sequencedata_base_year;
		} else {
			o.sequencedata_mean_log_year = Double.NEGATIVE_INFINITY;
			o.sequencedata_base_year = s.sequencedata_year;
		}
		o.sequencedata_mean_init_log_factor = s.sequencedata_init_log_factor;
		o.sequencedata_mean_slope_log_factor = s.sequencedata_slope_log_factor;
		return o;
	}

	static MathUtil.ScienceSpeedScenario applyScienceSpeedZValues (MathUtil.ScienceSpeedModelParameters m, MathUtil.ScienceSpeedZValues z) {
		MathUtil.ScienceSpeedScenario o = new MathUtil.ScienceSpeedScenario();
		o.institutional_slope_log_factor = 
			  m.institutional_mean_slope_log_factor
			+ z.institutional_rate_Z*m.institutional_stddev_slope_log_factor;
		o.population_year =
			  m.population_base_year
			+ Math.exp(
				  m.population_mean_log_year
				+ z.population_year_Z*m.population_stddev_log_year);
		o.sequencedata_year =
			  m.sequencedata_base_year
			+ Math.exp(
				  m.sequencedata_mean_log_year
				+ z.sequencedata_year_Z*m.sequencedata_stddev_log_year);
		o.population_init_log_rate = 
			  m.population_mean_init_log_rate
			+ z.population_rate_Z*m.population_stddev_init_log_rate;
		o.population_slope_log_rate = 
			  m.population_mean_slope_log_rate
			+ z.population_rate_Z*m.population_stddev_slope_log_rate;
		o.sequencedata_init_log_factor = 
			+ m.sequencedata_mean_init_log_factor
			+ z.sequencedata_factor_Z*m.sequencedata_stddev_init_log_factor;
		o.sequencedata_slope_log_factor =
			+ m.sequencedata_mean_slope_log_factor
			+ z.sequencedata_factor_Z*m.sequencedata_stddev_slope_log_factor;
		return o;
	}

	static MathUtil.ScienceSpeedScenarioReduced applyScienceSpeedZValuesReduce (MathUtil.ScienceSpeedModelParameters m, MathUtil.ScienceSpeedZValues z) {
		MathUtil.ScienceSpeedScenarioReduced o = new MathUtil.ScienceSpeedScenarioReduced();
		o.institutional_slope_log_factor = 
			  m.institutional_mean_slope_log_factor
			+ z.institutional_rate_Z*m.institutional_stddev_slope_log_factor;
		o.sciencetalent_rel_year =
			  m.population_base_year
			+ Math.exp(
				  m.population_mean_log_year
				+ z.population_year_Z*m.population_stddev_log_year)
			- m.institutional_base_year;
		double sequencedata_rel_year =
			  m.sequencedata_base_year
			+ Math.exp(
				  m.sequencedata_mean_log_year
				+ z.sequencedata_year_Z*m.sequencedata_stddev_log_year)
			- m.institutional_base_year;
		o.sciencetalent_init_log_rate = 
			  m.population_mean_init_log_rate
			+ z.population_rate_Z*m.population_stddev_init_log_rate
			+ m.sequencedata_mean_init_log_factor
			+ (o.sciencetalent_rel_year - sequencedata_rel_year)
			  * (m.sequencedata_mean_slope_log_factor)
			+ z.sequencedata_factor_Z * (
			    m.sequencedata_stddev_init_log_factor
			  + (o.sciencetalent_rel_year - sequencedata_rel_year)
			    * (m.sequencedata_stddev_slope_log_factor)
			  );
		o.sciencetalent_slope_log_rate =
			  m.population_mean_slope_log_rate
			+ z.population_rate_Z*m.population_stddev_slope_log_rate
			+ m.sequencedata_mean_slope_log_factor
			+ z.sequencedata_factor_Z*m.sequencedata_stddev_slope_log_factor;
		return o;
	}
	
	static void logDoubleList(double[] l) {
		StringBuilder b = new StringBuilder();
		for (int i=0; i<l.length; ++i) {
			b.append(String.format("%20.16g ", l[i]));
		}		
		org.singinst.uf.common.LogUtil.info(b.toString());
	}

	static void logDoubleList(String s, double[] l) {
		StringBuilder b = new StringBuilder(s);
		for (int i=0; i<l.length; ++i) {
			b.append(String.format("%20.16g ", l[i]));
		}		
		org.singinst.uf.common.LogUtil.info(b.toString());
	}
	
	static double[] projectScienceYearsNumerical(MathUtil.ScienceSpeedScenario s, double[] t) throws org.apache.commons.math.FunctionEvaluationException,org.apache.commons.math.ConvergenceException {
		double[] o = new double[t.length];
		final double t0_seq, m_seq, b_seq;
		final double t0_pop, m_pop, b_pop;
		final double t0_inst, m_inst;
		
		t0_seq = s.sequencedata_year;
		m_seq = s.sequencedata_slope_log_factor;
		b_seq = s.sequencedata_init_log_factor;
		t0_pop = s.population_year;
		m_pop = s.population_slope_log_rate;
		b_pop = s.population_init_log_rate;
		t0_inst = s.institutional_base_year;
		m_inst = s.institutional_slope_log_factor;
		
		UnivariateRealFunction science_rate_after_pop = new UnivariateRealFunction() {
			public double value(double t) {
				return (t-t0_pop)*MathUtil.exprel((m_seq+m_pop)*(t-t0_pop))*Math.exp(m_inst*(t-t0_inst)+b_pop+m_seq*(t0_pop-t0_seq)+b_seq);
				//return (t-t0_pop)*MathUtil.exprel(m_pop*(t-t0_pop))*Math.exp(b_pop+m_inst*(t-t0_inst)+b_seq+m_seq*(t-t0_seq));
			}
		};
		double baseyear = s.institutional_base_year;
		double prevyear = t0_pop; double science_from_pop = 0;
		double year; double science;
		org.apache.commons.math.analysis.UnivariateRealIntegrator science_from_pop_integrator = new RombergIntegrator(science_rate_after_pop);
		science_from_pop_integrator.setRelativeAccuracy(32*MathUtil.DOUBLE_EPS);
		for (int i=0; i<t.length; ++i) {
			year = t[i];
			science = (year-baseyear)*MathUtil.exprel(m_inst*(year-baseyear));
			if (year >= t0_pop) {
				if (year > prevyear) {
					science_from_pop = science_from_pop + science_from_pop_integrator.integrate(prevyear, year);
					prevyear = year;
				} else if (year < prevyear) {
					science_from_pop = science_from_pop - science_from_pop_integrator.integrate(year, prevyear);
					prevyear = year;					
				}
				science = science + science_from_pop;
			}
			o[i] = science + baseyear;
		}
		return o;
	}
	
	static double[] projectScienceYearsAnalytic(MathUtil.ScienceSpeedScenario s, double[] rescheduled_years) {
		double[] o = new double[rescheduled_years.length];
		double dt0_seq, m_seq, b_seq;
		double dt0_pop, m_pop, b_pop;
		double t0_inst, m_inst;
		
		double dt, science_years;
		
		t0_inst = s.institutional_base_year;
		m_inst = s.institutional_slope_log_factor;
		dt0_seq = s.sequencedata_year - t0_inst;
		m_seq = s.sequencedata_slope_log_factor;
		b_seq = s.sequencedata_init_log_factor;
		dt0_pop = s.population_year - t0_inst;
		m_pop = s.population_slope_log_rate;
		b_pop = s.population_init_log_rate;
		
		for (int j=0; j<rescheduled_years.length; ++j) {
			dt = rescheduled_years[j] - t0_inst;
			science_years = MathUtil.exprel( m_inst * dt ) * dt;
			if (dt > dt0_pop)
				science_years += (dt-dt0_pop)*(dt-dt0_pop) * Math.exp( m_inst*((1*dt0_pop+2*dt)/3/* - 0 */) + m_seq*((2*dt0_pop+1*dt)/3-dt0_seq)+b_seq + m_pop*((2*dt0_pop+1*dt)/3-dt0_pop)+b_pop ) * MathUtil.rpn_b_expm1_a_t_a_expm1_b_t_m_a_b_p_3_d_exp_d_a_d_b_d_a_b_m_d( -m_inst*(dt-dt0_pop), (m_pop+m_seq)*(dt-dt0_pop) ) ;
			o[j] = science_years + t0_inst;
		}
		return o;
	}

	static double[] projectScienceYearsAnalytic(MathUtil.ScienceSpeedScenarioReduced s, double[] reschedule_rel_years, double base_year) {
		double[] o = new double[reschedule_rel_years.length];
		double dt0_tal, m_tal, b_tal;
		double m_inst;
		
		double dt, science_years;
		
		m_inst = s.institutional_slope_log_factor;
		dt0_tal = s.sciencetalent_rel_year;
		m_tal = s.sciencetalent_slope_log_rate;
		b_tal = s.sciencetalent_init_log_rate;
		
		for (int j=0; j<reschedule_rel_years.length; ++j) {
			dt = reschedule_rel_years[j] - base_year;
			science_years = MathUtil.exprel( m_inst * dt ) * dt;
			if (dt > dt0_tal)
				science_years += (dt-dt0_tal)*(dt-dt0_tal) * Math.exp( m_inst*((1*dt0_tal+2*dt)/3/* - 0 */) + m_tal*((2*dt0_tal+1*dt)/3-dt0_tal)+b_tal ) * MathUtil.rpn_b_expm1_a_t_a_expm1_b_t_m_a_b_p_3_d_exp_d_a_d_b_d_a_b_m_d( -m_inst*(dt-dt0_tal), m_tal*(dt-dt0_tal) ) ;
			o[j] = science_years + base_year;
		}
		return o;
	}

	static public void testShowMeWhatClogIntervalLooksLike() throws Exception {
		double d = MathUtil.clogIntervalAverage(0, 4);
		org.singinst.uf.common.LogUtil.info(new Double(d).toString());
	}
	
	void notclearonthis() throws Exception {
		for (int i=0; i<2; ++i) {
		final double d = MathUtil.RANDOM.nextDouble();
		UnivariateRealFunction f = new UnivariateRealFunction() {
			public double value(double x) {
				return d*x;
			}
		};
		RombergIntegrator ri = new RombergIntegrator(f);

		org.singinst.uf.common.LogUtil.info(String.format("%g %g", d/2, ri.integrate(0, 1)));
		}
	}
	
	public void testMeh() throws Exception {
		{

		notclearonthis();
		notclearonthis();
			
		MathUtil.EventDiscreteDistributionSchedule s = new MathUtil.EventDiscreteDistributionSchedule();
		double[] s_dot_time = {0, 1, 2, 3, 4, 4.1, 4.2, 4.3, 4.4};
		double[] t = {-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.05,4.1,4.15,4.2,4.25,4.3,4.35,4.4,4.45};
		
		s.time = s_dot_time;
		double[] s_dot_clogProbEvent = {0,1,2,3,4,5,6,7,8};
		s.clogProbEvent = s_dot_clogProbEvent;
		s.logitSubevents = new double[0][];
		MathUtil.EventDiscreteDistributionSchedule r = s.logitnerp(t);
		StringBuilder b = new StringBuilder();
		for (int i=0; i<r.clogProbEvent.length; ++i) {
			b.append(r.clogProbEvent[i]).append(" ");
		}
		org.singinst.uf.common.LogUtil.info(b.toString());
		}
		
		{
		MathUtil.EventDiscreteDistributionSchedule s = heavisideSchedule(1);
		double[] t = {0,0.5,1,1.5,2};
		MathUtil.EventDiscreteDistributionSchedule r = s.logitnerp(t);
		logDoubleList(r.clogProbEvent);
		}
	}	

	static double LOG_DOUBLE_EPS = Math.log(MathUtil.DOUBLE_EPS);
	static double square_clogged(double clog_p) {
		if (clog_p >= -LOG_DOUBLE_EPS) return clog_p - MathUtil.LOG_TWO;
		if (clog_p < LOG_DOUBLE_EPS) return 2*clog_p;
		return -Math.log(Math.exp(-clog_p)*(2-Math.exp(-clog_p)));
	}
	static double[] square_clogged(double[] clog_p) {
		double[] o = new double[clog_p.length];
		for (int i=0; i<clog_p.length; ++i) {
			if (clog_p[i] >= -LOG_DOUBLE_EPS) o[i] = clog_p[i] - MathUtil.LOG_TWO;
			if (clog_p[i] < LOG_DOUBLE_EPS) o[i] = 2*clog_p[i];
			o[i] = -Math.log(Math.exp(-clog_p[i])*(2-Math.exp(-clog_p[i])));			
		}
		return o;
	}
	
	static double PHI = 0.618033988749894848d;
	static double COS_TWO_THIRDS_PI_PHI =  0.272883452047768708d;
	static double SIN_TWO_THIRDS_PI_PHI =  0.962047099469923622d;
	
	interface ScenarioSource { MathUtil.ScienceSpeedScenarioReduced scenario(int i); };
	
	public void testWhatever() throws Exception {
		int iters = 1000;
		MathUtil.ScienceSpeedModelParameters m = nullScienceModel();
		MathUtil.ScienceSpeedScenario sss = new MathUtil.ScienceSpeedScenario();
		MathUtil.ScienceSpeedScenarioReduced sssr = new MathUtil.ScienceSpeedScenarioReduced();
		MathUtil.ScienceSpeedZValues z = new MathUtil.ScienceSpeedZValues();
		
		sss.institutional_base_year = -MathUtil.RANDOM.nextDouble();
		//sss.institutional_base_year = -0.5;
		sss.institutional_slope_log_factor = MathUtil.RANDOM.nextGaussian();
		//sss.institutional_slope_log_factor = Math.log(1.0001)*10;
		sss.sequencedata_init_log_factor = MathUtil.RANDOM.nextGaussian();
		sss.sequencedata_slope_log_factor = MathUtil.RANDOM.nextGaussian();
		//sss.sequencedata_slope_log_factor = Math.log(1000);
		sss.sequencedata_year = MathUtil.RANDOM.nextDouble();
		sss.population_init_log_rate = MathUtil.RANDOM.nextGaussian();
		// //sss.population_init_log_rate = Math.log(10000);
		sss.population_slope_log_rate = MathUtil.RANDOM.nextGaussian();
		sss.population_year = MathUtil.RANDOM.nextDouble();
		
		MathUtil.EventDiscreteDistributionSchedule s = heavisideSchedule(1);
		m.institutional_base_year = sss.institutional_base_year;
		m.population_base_year = 0;
		m.population_mean_log_year = 0;
		m.population_stddev_log_year = 1;
		m.sequencedata_base_year = 0;
		m.sequencedata_mean_log_year = 0;
		m.sequencedata_stddev_log_year = 1;
		m = pointMassScienceSpeedModel(m, sss);
		//m.institutional_base_year = sss.institutional_base_year;
		sss = applyScienceSpeedZValues(m,z);
		sssr = applyScienceSpeedZValuesReduce(m,z);
		double[] t = {0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,2.0};
		double[] clog_weird
		           = {0  ,0  ,1  ,1  ,2  ,2  ,3  ,3  ,4  ,4  ,5  ,5  ,6  ,6  ,7  ,7  ,8  };
		
		logDoubleList(t);
		logDoubleList(projectScienceYearsNumerical(sss,t));
		logDoubleList(projectScienceYearsAnalytic(sss,t));
		logDoubleList(projectScienceYearsAnalytic(sssr,t,m.institutional_base_year));
		s.time = t;
		s.clogProbEvent = t;
		//s.clogProbEvent = clog_weird;
		MathUtil.EventDiscreteDistributionSchedule r = MathUtil.clogMarginalScienceSpeedEventRescheduling(m, s, t, iters);
		logDoubleList(r.clogProbEvent);
		
		r = s;
		MathUtil.EventDiscreteDistributionSchedule r_regular = s;
		MathUtil.EventDiscreteDistributionSchedule r_montecarlo = s;
		MathUtil.ScienceSpeedModelParameters m_v = m.clone();

		double wt = 0, wt_i;
		double[] squared_clog = new double[t.length];
		
		m_v.institutional_stddev_slope_log_factor = 1;
		for (int i=0; i<=200; ++i) {
			z.institutional_rate_Z = (i-100)/10d;
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = Math.exp(-z.institutional_rate_Z*z.institutional_rate_Z/2);
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
		}
		r_montecarlo = MathUtil.clogMarginalScienceSpeedEventRescheduling(m_v, s, t, iters);
		logDoubleList(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year));
		logDoubleList(s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)).clogProbEvent);
		logDoubleList("inst rate Z by regular", r_regular.clogProbEvent);
		logDoubleList("inst rate Z by regular", r_montecarlo.clogProbEvent);
		m_v.institutional_stddev_slope_log_factor = 0;
		z.institutional_rate_Z = 0;
		
		m_v.population_stddev_init_log_rate = 1;
		wt = 0;
		for (int i=0; i<=200; ++i) {
			z.population_rate_Z = (i-100)/10d;
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = Math.exp(-z.population_rate_Z*z.population_rate_Z/2);
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
		}
		r_montecarlo = MathUtil.clogMarginalScienceSpeedEventRescheduling(m_v, s, t, iters);
		logDoubleList("pop init rate Z by regular", r_regular.clogProbEvent);
		logDoubleList("pop init rate Z by MC     ", r_montecarlo.clogProbEvent);
		m_v.population_stddev_init_log_rate = 0;
		z.population_rate_Z = 0;
		
		m_v.population_stddev_slope_log_rate = 1;
		wt = 0;
		for (int i=0; i<=200; ++i) {
			z.population_rate_Z = (i-100)/10d;
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = Math.exp(-z.population_rate_Z*z.population_rate_Z/2);
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
		}
		r_montecarlo = MathUtil.clogMarginalScienceSpeedEventRescheduling(m_v, s, t, iters);
		logDoubleList("pop slope rate Z by regular", r_regular.clogProbEvent);
		logDoubleList("pop slope rate Z by MC     ", r_montecarlo.clogProbEvent);
		m_v.population_stddev_slope_log_rate = 0;
		z.population_rate_Z = 0;

		m_v.population_stddev_log_year = 1;
		wt = 0;
		for (int i=0; i<=200; ++i) {
			z.population_year_Z = (i-100)/10d;
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = Math.exp(-z.population_year_Z*z.population_year_Z/2);
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
		}
		logDoubleList("pop year Z by regular Z", r_regular.clogProbEvent);
		for (int i=0; i<iters; ++i) {
			z.population_year_Z = MathUtil.inverseCumulativeProbability(0, 1, (i+0.5d)/iters);
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = 1;
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
			//org.singinst.uf.common.LogUtil.info(String.format("regular %g leads to progress %g", sssr.sciencetalent_rel_year, rescheduled_t[rescheduled_t.length-1]));
		}
		logDoubleList("pop year Z by regular q", r_regular.clogProbEvent);
		r_montecarlo = MathUtil.clogMarginalScienceSpeedEventRescheduling(m_v, s, t, iters);
		logDoubleList("pop year Z by MC       ", r_montecarlo.clogProbEvent);
		m_v.population_stddev_log_year = 0;
		z.population_year_Z = 0;
		
		m_v.sequencedata_stddev_init_log_factor = 1;
		wt = 0;
		for (int i=0; i<=200; ++i) {
			z.sequencedata_factor_Z = (i-100)/10d;
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = Math.exp(-z.sequencedata_factor_Z*z.sequencedata_factor_Z/2);
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
		}
		r_montecarlo = MathUtil.clogMarginalScienceSpeedEventRescheduling(m_v, s, t, iters);
		logDoubleList("seqdata init factor Z by regular", r_regular.clogProbEvent);
		logDoubleList("seqdata init factor Z by MC     ", r_montecarlo.clogProbEvent);
		m_v.sequencedata_stddev_init_log_factor = 0;
		z.sequencedata_factor_Z = 0;
		
		m_v.sequencedata_stddev_slope_log_factor = 1;
		wt = 0;
		for (int i=0; i<=200; ++i) {
			z.sequencedata_factor_Z = (i-100)/10d;
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = Math.exp(-z.sequencedata_factor_Z*z.sequencedata_factor_Z/2);
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
		}
		r_montecarlo = MathUtil.clogMarginalScienceSpeedEventRescheduling(m_v, s, t, iters);
		logDoubleList("seqdata slope factor Z by regular", r_regular.clogProbEvent);
		logDoubleList("seqdata slope factor Z by MC     ", r_montecarlo.clogProbEvent);
		m_v.sequencedata_stddev_slope_log_factor = 0;
		z.sequencedata_factor_Z = 0;

		m_v.sequencedata_stddev_log_year = 1;
		wt = 0;
		for (int i=0; i<=200; ++i) {
			z.sequencedata_year_Z = (i-100)/10d;
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = Math.exp(-z.sequencedata_year_Z*z.sequencedata_year_Z/2);
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
		}
		r_montecarlo = MathUtil.clogMarginalScienceSpeedEventRescheduling(m_v, s, t, iters);
		logDoubleList("seqdata year Z by regular Z", r_regular.clogProbEvent);
		for (int i=0; i<iters; ++i) {
			z.sequencedata_year_Z = MathUtil.inverseCumulativeProbability(0, 1, (i+0.5d)/iters);
			sssr = applyScienceSpeedZValuesReduce(m_v,z);
			wt_i = 1;
			double[] rescheduled_t = projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year);
			MathUtil.EventDiscreteDistributionSchedule s_i = s.logitnerp(rescheduled_t);
			r_regular = MathUtil.EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s_i);
			//r_regular = EventDiscreteDistributionSchedule.weightedMix(wt, wt_i, r_regular, s.logitnerp(projectScienceYearsAnalytic(sssr,t,m_v.institutional_base_year)));
			wt = wt + wt_i;
			//org.singinst.uf.common.LogUtil.info(String.format("regular %g leads to progress %g", sssr.sciencetalent_rel_year, rescheduled_t[rescheduled_t.length-1]));
		}
		logDoubleList("seqdata year Z by regular q", r_regular.clogProbEvent);
		logDoubleList("seqdata year Z by MC       ", r_montecarlo.clogProbEvent);
		m_v.sequencedata_stddev_log_year = 0;
		z.sequencedata_year_Z = 0;
		
		
		fail();
	}
	
	static double SQRT_EPS = Math.sqrt(MathUtil.DOUBLE_EPS);
	public double u(double x) {
		return x+SQRT_EPS*(Math.abs(x)+MathUtil.DOUBLE_MIN_NORMAL);
	}
	public double d(double x) {
		return x-SQRT_EPS*(Math.abs(x)/*+MathUtil.DOUBLE_MIN_NORMAL*/);
	}
	
	public void testWeightedClogMix() throws Exception {
		int i,j,k,l;
		double o, cl_p;
		for (i=0; i<6; ++i) for (j=0; j<6; ++j) for (k=0; k<6; ++k) for (l=0; l<6; ++l) {
			o = MathUtil.weightedClogMix(i, j, k, l);
			cl_p = MathUtil.clog((i*MathUtil.expc(k)+j*MathUtil.expc(l))/(i+j));
			Assert.assertEquals(cl_p, o, 1e-14d*cl_p);
			//LogUtil.info(String.format("%d %d %d %d %20.14g %20.14g %20.14g %20.14g", i,j,k,l, o, p, o-p, Math.exp(-((o-p)%100)) ));
		}
		double a,b,c,d;
		double w0, w1, clog_p0, clog_p1;
		double cl_p_hi, cl_p_lo;
		for (a=-50; a<=50; a+=25) for (b=-50; b<50; b+=25) for (c=0; c<=100; c+=25) for (d=0; d<=100; d+=25) {
			w0 = Math.exp(a); clog_p0 = c;
			w1 = Math.exp(b); clog_p1 = d;
			o = MathUtil.weightedClogMix(w0, w1, clog_p0, clog_p1);
//			cl_p = clog_p0-Math.log1p((w1*Math.expm1(clog_p0-clog_p1))/(w0+w1));
//			Assert.assertEquals(cl_p, o, 16*MathUtil.DOUBLE_EPS*(
//					  Math.abs( clog_p0 )
//					+ Math.abs( Math.log1p((w1*Math.expm1(clog_p0-clog_p1))/(w0+w1)) )
//					+ Math.abs( 1/(1+(w1*Math.expm1(clog_p0-clog_p1))/(w0+w1)) )
//				) );
			cl_p = clog_p0-Math.log1p((Math.expm1(clog_p0-clog_p1))/(w0/w1+1));
			cl_p_hi = u(u(clog_p0)-d(Math.log1p(Math.max(1,d(d(Math.expm1(d(d(clog_p0)-u(clog_p1))))/(w0/w1+1))))));
			cl_p_lo = d(d(clog_p0)-u(Math.log1p(           u(u(Math.expm1(u(u(clog_p0)-d(clog_p1))))/(w0/w1+1)) )));
			try {
				Assert.assertEquals(cl_p, o, 16*SQRT_EPS*(Math.abs(cl_p_hi-cl_p)+Math.abs(cl_p_lo-cl_p)));
			} catch (junit.framework.AssertionFailedError e) {
				o = MathUtil.weightedClogMix(w0, w1, clog_p0, clog_p1);
				throw(e);
			}
			
//			cl_p = clog_p1-Math.log1p((w0*Math.expm1(clog_p1-clog_p0))/(w0+w1));
//			Assert.assertEquals(cl_p, o, 16*MathUtil.DOUBLE_EPS*(
//					  Math.abs( clog_p1 )
//					+ Math.abs( Math.log1p((w0*Math.expm1(clog_p1-clog_p0))/(w0+w1)) )
//					+ Math.abs( 1/(1+(w0*Math.expm1(clog_p1-clog_p0))/(w0+w1)) )
//				) );
			cl_p = clog_p1-Math.log1p((Math.expm1(clog_p1-clog_p0))/(w1/w0+1));
			cl_p_hi = u(u(clog_p1)-d(Math.log1p(Math.max(1,d(d(Math.expm1(d(d(clog_p1)-u(clog_p0))))/(w1/w0+1))))));
			cl_p_lo = d(d(clog_p1)-u(Math.log1p(           u(u(Math.expm1(u(u(clog_p1)-d(clog_p0))))/(w1/w0+1)) )));
			try {
				Assert.assertEquals(cl_p, o, 16*SQRT_EPS*(Math.abs(cl_p_hi-cl_p)+Math.abs(cl_p_lo-cl_p)));
			} catch (junit.framework.AssertionFailedError e) {
				o = MathUtil.weightedClogMix(w0, w1, clog_p0, clog_p1);
				throw(e);
			}

		}
	}
	
	public void test() throws Exception {
	}
}