247.sagews

var('T1, T2, T3, T13, T12, T23, T123')
assume(T12 >0); assume(T13 >0); assume(T23 >0); assume(T123>0)

#  only main effects
p1 = [1/(1+T1), 1/(1+T1*T13), 1/(1+T1*T12), 1/(1+T1*T12*T13), T1/(1+T1), (T1*T13)/(1+T1*T13), (T1*T12)/(1+T1*T12),(T1*T12*T13)/(1+T1*T12*T13)]
p2 = [1/(1+T2), 1/(1+T2*T23), T2/(1+T2), (T2*T23)/(1+T2*T23),
1/(1+T2), 1/(1+T2*T23), T2/(1+T2), (T2*T23)/(1+T2*T23)]
p3 = [1/(1+T3), T3/(1+T3),1/(1+T3), T3/(1+T3), 1/(1+T3), T3/(1+T3),1/(1+T3), T3/(1+T3)]

P = [p1[i] * p2[i] * p3[i]  for i in range(len(p1))]
P = Matrix(P, nrows=8)
P

[ 1/((T1 + 1)*(T2 + 1)*(T3 + 1))] [ T3/((T1*T13 + 1)*(T2*T23 + 1)*(T3 + 1))] [ T2/((T1*T12 + 1)*(T2 + 1)*(T3 + 1))] [ T2*T23*T3/((T1*T12*T13 + 1)*(T2*T23 + 1)*(T3 + 1))] [ T1/((T1 + 1)*(T2 + 1)*(T3 + 1))] [ T1*T13*T3/((T1*T13 + 1)*(T2*T23 + 1)*(T3 + 1))] [ T1*T12*T2/((T1*T12 + 1)*(T2 + 1)*(T3 + 1))] [T1*T12*T13*T2*T23*T3/((T1*T12*T13 + 1)*(T2*T23 + 1)*(T3 + 1))]

I = identity_matrix(2); u = Matrix([1,1], ncols=2)
M = I.tensor_product(u); C = M.tensor_product(I)

P13 = C * P

odr = (P13[0]*P13[3])/((P13[1]*P13[2]))
odr = odr.substitute(T1 = 1, T2 = 1, T3 = 1)
odr = (factor(odr.numerator())/factor(odr.denominator()))

f =solve(odr==1, T13)
m = f[1]
m = m.rhs()
m

(sqrt((T12^4 + 12*T12^3 + 38*T12^2 + 12*T12 + 1)*T23^2 + (-6*T12^4 + 24*T12^3 + 92*T12^2 + 24*T12 - 6)*T23 + 9*T12^4 + 12*T12^3 + 22*T12^2 + 12*T12 + 9) + (1 - T12^2)*T23 + 3*T12^2 - 3)/((2*T12^2 + 6*T12)*T23 + 2*T12^2 + 6*T12)

max = 2
a = parametric_plot3d([T12, T23, m], (T12,0.001, max), (T23,0.001,max), plot_points=[50,50]) 
b = parametric_plot3d([T12,T23,1], (T12,0.01, max), (T23, 0.01,max),plot_points=[50,50], color="grey", opacity=.7)
x,y,z=var('x,y,z')
c = parametric_plot3d([1,y,z], (y,0.01, max), (z, 0,3),plot_points=[50,50], color="grey", opacity=.7)
d = parametric_plot3d([x,1,z], (x,0.01, max), (z, 0,3),plot_points=[50,50], color="grey",  opacity=.7)
p = a+b+c+d
p.show(figsize=[8,8])