Thanks, transform=pm.transforms.ordered in the constructor works, but for anyone else that uses this also values for testval are also needed when you provide this constraint.
One more related question, what about adding constraints on what is provided pm.Normal.dist()? This question really stems from not fully understanding what is provided by the dist() method. I wanted to do something like this:
# Here depth is the predictor variable
dist1 = pm.Normal(mu=(c1[0] + c2[0] * depth), sigma=std).dist()
dist2 = pm.Normal(mu=(c1[1] + c2[1] * depth), sigma=std).dist()
# Want dist2 to always been greater than dist1. This doesn't work
ordered = pm.Potential('ordered', tt.switch(dists1 - dists2 < 0, -np.inf, 0))
I am essentially try to perform both a classification and regression simultaneously.