reset;
option snopt_options "major_optimality_tolerance=1e-16**0.5 outlev=1";
option output_precision 16;
option display_precision 16;
option objective_precision 16;
option solution_precision 16;
option presolve_eps ((1e-16)**2);

param MIN_PARLAY_SIZE > 0 integer;
param MAX_PARLAY_SIZE >= MIN_PARLAY_SIZE integer;
param COMMISSION < 1 default 0;
param KELLYMULT > 0 default 1;
param DISPLAY_INCR > 0, <= 1 default 1e-012;
param DISPLAY_MIN >= 0, <= 1 default 0;
param OPTIMIZE_ROUND_LOTS binary, default 0;
param OPTIMIZE_MIN_BETS binary, default 0;

param n_events > 0  integer;

data;

param COMMISSION = 0;
param MIN_PARLAY_SIZE = 2;
param MAX_PARLAY_SIZE = 15;
param KELLYMULT = 1;
param DISPLAY_INCR = 1e-022;
param DISPLAY_MIN = 0;
param n_events = 9;


model;
param n_outcomes := 2**n_events ;
set OUTCOMES := { 0 .. n_outcomes-1 } ;
set EVENTS   :=  { 0 .. n_events-1 };
set POW {j in OUTCOMES} := {i in EVENTS: (j div 2**i) mod 2 = 1};

param pr { i in EVENTS } >= 0, <= 1;	# single probabilities
param od { i in EVENTS } > 1;		# single odds

param lb { j in OUTCOMES } >= 0, <= 1, default 0;	# lower bounds
param ub { j in OUTCOMES } >= lb[j], <= 1, default 0;	# upper bounds

param outc_p { j in OUTCOMES } >= 0, <= 1 default 1;	# outcome (compund event) probabilities
param par_od { j in OUTCOMES } >= 1 default 1;		# parlay (compound bet) odds

param outc_p_mult { j in OUTCOMES } default 1;			# outcome (compund event) probabilities
#param custom_par_od { j in OUTCOMES } default -1;		# parlay (compound bet) odds

param bet_names { i in OUTCOMES } symbolic;

data;
param :		lb	ub	outc_p_mult	:=
	3	0	1	1
	5	0	1	1
	6	0	1	1
	7	0	1	1
	9	0	1	1
	10	0	1	1
	11	0	1	1
	12	0	1	1
	13	0	1	1
	14	0	1	1
	15	0	1	1
	17	0	1	1
	18	0	1	1
	19	0	1	1
	20	0	1	1
	21	0	1	1
	22	0	1	1
	23	0	1	1
	24	0	1	1
	25	0	1	1
	26	0	1	1
	27	0	1	1
	28	0	1	1
	29	0	1	1
	30	0	1	1
	31	0	1	1
	33	0	1	1
	34	0	1	1
	35	0	1	1
	36	0	1	1
	37	0	1	1
	38	0	1	1
	39	0	1	1
	40	0	1	1
	41	0	1	1
	42	0	1	1
	43	0	1	1
	44	0	1	1
	45	0	1	1
	46	0	1	1
	47	0	1	1
	48	0	1	1
	49	0	1	1
	50	0	1	1
	51	0	1	1
	52	0	1	1
	53	0	1	1
	54	0	1	1
	55	0	1	1
	56	0	1	1
	57	0	1	1
	58	0	1	1
	59	0	1	1
	60	0	1	1
	61	0	1	1
	62	0	1	1
	63	0	1	1
	65	0	1	1
	66	0	1	1
	67	0	1	1
	68	0	1	1
	69	0	1	1
	70	0	1	1
	71	0	1	1
	72	0	1	1
	73	0	1	1
	74	0	1	1
	75	0	1	1
	76	0	1	1
	77	0	1	1
	78	0	1	1
	79	0	1	1
	80	0	1	1
	81	0	1	1
	82	0	1	1
	83	0	1	1
	84	0	1	1
	85	0	1	1
	86	0	1	1
	87	0	1	1
	88	0	1	1
	89	0	1	1
	90	0	1	1
	91	0	1	1
	92	0	1	1
	93	0	1	1
	94	0	1	1
	95	0	1	1
	96	0	1	1
	97	0	1	1
	98	0	1	1
	99	0	1	1
	100	0	1	1
	101	0	1	1
	102	0	1	1
	103	0	1	1
	104	0	1	1
	105	0	1	1
	106	0	1	1
	107	0	1	1
	108	0	1	1
	109	0	1	1
	110	0	1	1
	111	0	1	1
	112	0	1	1
	113	0	1	1
	114	0	1	1
	115	0	1	1
	116	0	1	1
	117	0	1	1
	118	0	1	1
	119	0	1	1
	120	0	1	1
	121	0	1	1
	122	0	1	1
	123	0	1	1
	124	0	1	1
	125	0	1	1
	126	0	1	1
	127	0	1	1
	129	0	1	1
	130	0	1	1
	131	0	1	1
	132	0	1	1
	133	0	1	1
	134	0	1	1
	135	0	1	1
	136	0	1	1
	137	0	1	1
	138	0	1	1
	139	0	1	1
	140	0	1	1
	141	0	1	1
	142	0	1	1
	143	0	1	1
	144	0	1	1
	145	0	1	1
	146	0	1	1
	147	0	1	1
	148	0	1	1
	149	0	1	1
	150	0	1	1
	151	0	1	1
	152	0	1	1
	153	0	1	1
	154	0	1	1
	155	0	1	1
	156	0	1	1
	157	0	1	1
	158	0	1	1
	159	0	1	1
	160	0	1	1
	161	0	1	1
	162	0	1	1
	163	0	1	1
	164	0	1	1
	165	0	1	1
	166	0	1	1
	167	0	1	1
	168	0	1	1
	169	0	1	1
	170	0	1	1
	171	0	1	1
	172	0	1	1
	173	0	1	1
	174	0	1	1
	175	0	1	1
	176	0	1	1
	177	0	1	1
	178	0	1	1
	179	0	1	1
	180	0	1	1
	181	0	1	1
	182	0	1	1
	183	0	1	1
	184	0	1	1
	185	0	1	1
	186	0	1	1
	187	0	1	1
	188	0	1	1
	189	0	1	1
	190	0	1	1
	191	0	1	1
	192	0	1	1
	193	0	1	1
	194	0	1	1
	195	0	1	1
	196	0	1	1
	197	0	1	1
	198	0	1	1
	199	0	1	1
	200	0	1	1
	201	0	1	1
	202	0	1	1
	203	0	1	1
	204	0	1	1
	205	0	1	1
	206	0	1	1
	207	0	1	1
	208	0	1	1
	209	0	1	1
	210	0	1	1
	211	0	1	1
	212	0	1	1
	213	0	1	1
	214	0	1	1
	215	0	1	1
	216	0	1	1
	217	0	1	1
	218	0	1	1
	219	0	1	1
	220	0	1	1
	221	0	1	1
	222	0	1	1
	223	0	1	1
	224	0	1	1
	225	0	1	1
	226	0	1	1
	227	0	1	1
	228	0	1	1
	229	0	1	1
	230	0	1	1
	231	0	1	1
	232	0	1	1
	233	0	1	1
	234	0	1	1
	235	0	1	1
	236	0	1	1
	237	0	1	1
	238	0	1	1
	239	0	1	1
	240	0	1	1
	241	0	1	1
	242	0	1	1
	243	0	1	1
	244	0	1	1
	245	0	1	1
	246	0	1	1
	247	0	1	1
	248	0	1	1
	249	0	1	1
	250	0	1	1
	251	0	1	1
	252	0	1	1
	253	0	1	1
	254	0	1	1
	255	0	1	1
	257	0	1	1
	258	0	1	1
	259	0	1	1
	260	0	1	1
	261	0	1	1
	262	0	1	1
	263	0	1	1
	264	0	1	1
	265	0	1	1
	266	0	1	1
	267	0	1	1
	268	0	1	1
	269	0	1	1
	270	0	1	1
	271	0	1	1
	272	0	1	1
	273	0	1	1
	274	0	1	1
	275	0	1	1
	276	0	1	1
	277	0	1	1
	278	0	1	1
	279	0	1	1
	280	0	1	1
	281	0	1	1
	282	0	1	1
	283	0	1	1
	284	0	1	1
	285	0	1	1
	286	0	1	1
	287	0	1	1
	288	0	1	1
	289	0	1	1
	290	0	1	1
	291	0	1	1
	292	0	1	1
	293	0	1	1
	294	0	1	1
	295	0	1	1
	296	0	1	1
	297	0	1	1
	298	0	1	1
	299	0	1	1
	300	0	1	1
	301	0	1	1
	302	0	1	1
	303	0	1	1
	304	0	1	1
	305	0	1	1
	306	0	1	1
	307	0	1	1
	308	0	1	1
	309	0	1	1
	310	0	1	1
	311	0	1	1
	312	0	1	1
	313	0	1	1
	314	0	1	1
	315	0	1	1
	316	0	1	1
	317	0	1	1
	318	0	1	1
	319	0	1	1
	320	0	1	1
	321	0	1	1
	322	0	1	1
	323	0	1	1
	324	0	1	1
	325	0	1	1
	326	0	1	1
	327	0	1	1
	328	0	1	1
	329	0	1	1
	330	0	1	1
	331	0	1	1
	332	0	1	1
	333	0	1	1
	334	0	1	1
	335	0	1	1
	336	0	1	1
	337	0	1	1
	338	0	1	1
	339	0	1	1
	340	0	1	1
	341	0	1	1
	342	0	1	1
	343	0	1	1
	344	0	1	1
	345	0	1	1
	346	0	1	1
	347	0	1	1
	348	0	1	1
	349	0	1	1
	350	0	1	1
	351	0	1	1
	352	0	1	1
	353	0	1	1
	354	0	1	1
	355	0	1	1
	356	0	1	1
	357	0	1	1
	358	0	1	1
	359	0	1	1
	360	0	1	1
	361	0	1	1
	362	0	1	1
	363	0	1	1
	364	0	1	1
	365	0	1	1
	366	0	1	1
	367	0	1	1
	368	0	1	1
	369	0	1	1
	370	0	1	1
	371	0	1	1
	372	0	1	1
	373	0	1	1
	374	0	1	1
	375	0	1	1
	376	0	1	1
	377	0	1	1
	378	0	1	1
	379	0	1	1
	380	0	1	1
	381	0	1	1
	382	0	1	1
	383	0	1	1
	384	0	1	1
	385	0	1	1
	386	0	1	1
	387	0	1	1
	388	0	1	1
	389	0	1	1
	390	0	1	1
	391	0	1	1
	392	0	1	1
	393	0	1	1
	394	0	1	1
	395	0	1	1
	396	0	1	1
	397	0	1	1
	398	0	1	1
	399	0	1	1
	400	0	1	1
	401	0	1	1
	402	0	1	1
	403	0	1	1
	404	0	1	1
	405	0	1	1
	406	0	1	1
	407	0	1	1
	408	0	1	1
	409	0	1	1
	410	0	1	1
	411	0	1	1
	412	0	1	1
	413	0	1	1
	414	0	1	1
	415	0	1	1
	416	0	1	1
	417	0	1	1
	418	0	1	1
	419	0	1	1
	420	0	1	1
	421	0	1	1
	422	0	1	1
	423	0	1	1
	424	0	1	1
	425	0	1	1
	426	0	1	1
	427	0	1	1
	428	0	1	1
	429	0	1	1
	430	0	1	1
	431	0	1	1
	432	0	1	1
	433	0	1	1
	434	0	1	1
	435	0	1	1
	436	0	1	1
	437	0	1	1
	438	0	1	1
	439	0	1	1
	440	0	1	1
	441	0	1	1
	442	0	1	1
	443	0	1	1
	444	0	1	1
	445	0	1	1
	446	0	1	1
	447	0	1	1
	448	0	1	1
	449	0	1	1
	450	0	1	1
	451	0	1	1
	452	0	1	1
	453	0	1	1
	454	0	1	1
	455	0	1	1
	456	0	1	1
	457	0	1	1
	458	0	1	1
	459	0	1	1
	460	0	1	1
	461	0	1	1
	462	0	1	1
	463	0	1	1
	464	0	1	1
	465	0	1	1
	466	0	1	1
	467	0	1	1
	468	0	1	1
	469	0	1	1
	470	0	1	1
	471	0	1	1
	472	0	1	1
	473	0	1	1
	474	0	1	1
	475	0	1	1
	476	0	1	1
	477	0	1	1
	478	0	1	1
	479	0	1	1
	480	0	1	1
	481	0	1	1
	482	0	1	1
	483	0	1	1
	484	0	1	1
	485	0	1	1
	486	0	1	1
	487	0	1	1
	488	0	1	1
	489	0	1	1
	490	0	1	1
	491	0	1	1
	492	0	1	1
	493	0	1	1
	494	0	1	1
	495	0	1	1
	496	0	1	1
	497	0	1	1
	498	0	1	1
	499	0	1	1
	500	0	1	1
	501	0	1	1
	502	0	1	1
	503	0	1	1
	504	0	1	1
	505	0	1	1
	506	0	1	1
	507	0	1	1
	508	0	1	1
	509	0	1	1
	510	0	1	1
	511	0	1	1
;

param :		pr	od	:=
	0	0.552240923509213	1.90909090909091
	1	0.56928132367734	1.90909090909091
	2	0.538501235732385	1.90909090909091
	3	0.551422443698211	1.90909090909091
	4	0.543298226633478	1.90909090909091
	5	0.528363912577382	1.90909090909091
	6	0.538932400236479	1.90909090909091
	7	0.557452323834035	1.90909090909091
	8	0.572117534309842	1.90909090909091
;

model;
param KELLYEXP = 1 - 1/KELLYMULT;
param parlay_string symbolic;

#var x { j in OUTCOMES }  >= lb[j], <= ub[j] default 0;	# parlay weightings outcome 0 = null	
var y >= 1e-016 default 1;
var rr { 0 .. n_events } >= 0, <= 1 default 0;	# single rr stakes 

param par_indices { k in OUTCOMES } >= 0, <= n_outcomes - 1, integer;	# outcome number associated with k
param parlay_size { k in OUTCOMES } >= 0, <= n_events, integer;		# size of parlay associated with outcome k
#param n_parlays >= 0, <= n_outcomes - 1, integer, default 0;		# number of parlays
param n_parlays >= 0, <= n_outcomes, integer, default 0;		# number of parlays

param temp_num_wins integer, default 0;
param EV default -1;
param EG default 0;
#param total_prob >= 0 default 0;

for {j in 0..n_outcomes-1} {
	let temp_num_wins := 0;
	for {i in 0..n_events-1} {
		if j div 2**(n_events - 1 - i) mod 2 = 0 then {
			let outc_p[j] := (1-pr[i]) * outc_p[j];
		} else {
			let outc_p[j] := pr[i] * outc_p[j];
			let par_od[j] := par_od[j] * od[i];
			let temp_num_wins := temp_num_wins + 1;
		}
	}
	let outc_p[j] := outc_p_mult[j] * outc_p[j];
#	let total_prob := total_prob + outc_p[j];
	if temp_num_wins > MAX_PARLAY_SIZE || temp_num_wins < MIN_PARLAY_SIZE then {
		let ub[j] := lb[j];
	}
	if ub[j] >= 0 then {
		let n_parlays := n_parlays + 1;
		let par_indices[n_parlays-1] := j ;
		let parlay_size[n_parlays-1] := temp_num_wins;
		let par_od[j] := (par_od[j]-1)*(1-COMMISSION)+1;
	}
}


var x {j in 0..n_outcomes-1} = rr[parlay_size[j]];
var z {j in 0..n_outcomes-1} = 
  if KELLYMULT = 1 then (
      log(
        y + sum {k in 0..n_parlays-1} (
	   if card(POW[par_indices[k]] diff POW[j]) = 0 then  par_od[par_indices[k]] * x[par_indices[k]]
        )
      )
  ) else (
      (
        y + sum {k in 0..n_parlays-1} (
	   if card(POW[par_indices[k]] diff POW[j]) = 0 then  par_od[par_indices[k]] * x[par_indices[k]]
        )
      ) ** KELLYEXP
  )
;

maximize Util:
  (if KELLYMULT = 1 then 0 else -1/KELLYEXP) +
  (if KELLYMULT = 1 then 1 else 1/KELLYEXP) *
    sum {j in 0..n_outcomes-1} outc_p[j] * z[j]
#    sum {j in 0..n_outcomes-1} outc_p[j]/total_prob * z[j]
;

subj to Budget_Contraint:
  y = 1 - sum {k in 0..n_parlays-1} x[par_indices[k]];
subj to RR_Contraint2:
  rr[0]=0;
subj to BoundsUp {k in 0..n_parlays-1}:
  x[par_indices[k]] <= ub[k];
subj to BoundsDown {k in 0..n_parlays-1}:
  x[par_indices[k]] >= lb[k];


solve;

let y := 1;
printf "%-16s\t%s\n",":Bet", "Weight";
for { k in 0 .. n_parlays-1 } {
	let parlay_string := "";
	for {i in 0..n_events-1} {
		let parlay_string := parlay_string & ( if par_indices[k] div 2**(n_events - 1 - i) mod 2 = 0 then "" else i+1 & "_" );
	}
	if (0) then {
		let bet_names[par_indices[k]] := parlay_string;
	}

	if  x[par_indices[k]] > 1E-16 then {
		printf "%-16s\t%.016f\t%d\n", parlay_string, x[par_indices[k]], par_indices[k];
#		let x[par_indices[k]] :=
		let rr[parlay_size[k]] :=
			if not(OPTIMIZE_MIN_BETS) and x[par_indices[k]] < DISPLAY_MIN then (
				DISPLAY_MIN*round(x[par_indices[k]]/DISPLAY_MIN )
			) else if not(OPTIMIZE_ROUND_LOTS) then (
				DISPLAY_INCR*round(x[par_indices[k]]/DISPLAY_INCR )
			) else (
				x[par_indices[k]]
			)
		;
		let y := y - x[par_indices[k]];
	}
}

#let {k in 0 .. n_parlays-1}  x[par_indices[k]] :=
#	if not(OPTIMIZE_MIN_BETS) and x[par_indices[k]] < DISPLAY_MIN then (
#		DISPLAY_MIN*round(x[par_indices[k]]/DISPLAY_MIN )
#	) else if not(OPTIMIZE_ROUND_LOTS) then (
#		DISPLAY_INCR*round(x[par_indices[k]]/DISPLAY_INCR )
#	) else (
#		x[par_indices[k]]
#	)
#;
#let y := 1 - sum {k in 0..n_parlays-1} x[par_indices[k]];

for {j in 0..n_outcomes-1} {
#	let EV := EV + outc_p[j]/total_prob * if KELLYMULT = 1 then exp(z[j]) else (z[j]**(1/KELLYEXP));
#	let EG := EG + outc_p[j]/total_prob * if KELLYMULT = 1 then z[j] else log(z[j]**(1/KELLYEXP));
	let EV := EV + outc_p[j] * if KELLYMULT = 1 then exp(z[j]) else (z[j]**(1/KELLYEXP));
	let EG := EG + outc_p[j] * if KELLYMULT = 1 then z[j] else log(z[j]**(1/KELLYEXP));
}
let EV := EV;
let EG := exp(EG) - 1;

printf "EV: %.016f\n", EV;
printf "EG: %.016f\n", EG;
printf "Util: %.016f\n", Util;
printf "Total bet: %.016f\n", 1 - y;
printf "Solver time: %.02f\n", _total_solve_time;
printf "AMPL time: %.02f\n", _ampl_time + _total_solve_time;

option show_stats 1;