runIntegration ::
  forall x p deg n
  . ( Vectorize x, Vectorize p, KnownNat deg, KnownNat n )
  => Proxy n -> Proxy deg
  -> DirCollOptions
  -> (forall a . Floating a => x a -> p a -> a -> x a)
  -> x Double -> p Double -> Double
  -> IO (Either String (x Double))
runIntegration _ _ dirCollOpts ode x0 p tf = do
  let ocp :: OcpPhase' (IntegrationOcp x p)
      ocp =
        OcpPhase
        { ocpMayer = \_ _ _ _ _ _ -> 0
        , ocpLagrange = \_ _ _ _ _ _ _ _ -> 0
        , ocpDae = \x' x _ _ pp _ t -> ((ode x pp t) `minus` x', None)
        , ocpQuadratures = \_ _ _ _ _ _ _ _ -> None
        , ocpQuadratureOutputs = \_ _ _ _ _ _ _ _ -> None
        , ocpBc = \x0' _ _ _ _ _ -> x0'
        , ocpPathC = \_ _ _ _ _ _ _ -> None
        , ocpPlotOutputs = \_ _ _ _ _ _ _ _ _ _ _ -> None
        , ocpObjScale      = Nothing
        , ocpTScale        = Nothing
        , ocpXScale        = Nothing
        , ocpZScale        = Nothing
        , ocpUScale        = Nothing
        , ocpPScale        = Nothing
        , ocpResidualScale = Nothing
        , ocpBcScale       = Nothing
        , ocpPathCScale    = Nothing
        }
      ocpInputs :: OcpPhaseInputs' (IntegrationOcp x p)
      ocpInputs =
        OcpPhaseInputs
        { ocpPathCBnds = None
        , ocpBcBnds =  fmap (\x -> (Just x, Just x)) x0
        , ocpXbnd = pure (Nothing, Nothing)
        , ocpUbnd = None
        , ocpZbnd = None
        , ocpPbnd = fmap (\x -> (Just x, Just x)) p
        , ocpTbnd = (Just tf, Just tf)
        , ocpFixedP = None
        }
  let guess :: J (CollTraj x None None p n deg) (Vector Double)
      guess = cat $ makeGuessSim (collocationRoots dirCollOpts) tf x0 (\_ x _ -> ode x p 0) (\_ _ -> None) p
  cp  <- makeCollProblem dirCollOpts ocp ocpInputs guess :: IO (CollProblem x None None p x None x None None None None None n deg)
  (_, eopt) <- solveNlp "test_ocp_integrator" solver (cpNlp cp) Nothing
  return $ case eopt of
    Left m -> Left m
    Right opt -> Right (toXf (xOpt opt))
runIntegration ::
  forall x p deg n
  . ( Vectorize x, Vectorize p, KnownNat deg, KnownNat n )
  => Proxy n -> Proxy deg
  -> DirCollOptions
  -> (forall a . Floating a => x a -> p a -> a -> x a)
  -> x Double -> p Double -> Double
  -> IO (Either String (x Double))
runIntegration _ _ dirCollOpts ode x0 p tf = do
  let ocp :: OcpPhase' (IntegrationOcp x p)
      ocp =
        OcpPhase
        { ocpMayer = \_ _ _ _ _ _ -> 0
        , ocpLagrange = \_ _ _ _ _ _ _ _ -> 0
        , ocpDae = \x' x _ _ pp _ t -> ((ode x pp t) `minus` x', None)
        , ocpQuadratures = \_ _ _ _ _ _ _ _ -> None
        , ocpQuadratureOutputs = \_ _ _ _ _ _ _ _ -> None
        , ocpBc = \x0' _ _ _ _ _ -> x0'
        , ocpPathC = \_ _ _ _ _ _ _ -> None
        , ocpPlotOutputs = \_ _ _ _ _ _ _ _ _ _ _ -> None
        , ocpObjScale      = Nothing
        , ocpTScale        = Nothing
        , ocpXScale        = Nothing
        , ocpZScale        = Nothing
        , ocpUScale        = Nothing
        , ocpPScale        = Nothing
        , ocpResidualScale = Nothing
        , ocpBcScale       = Nothing
        , ocpPathCScale    = Nothing
        }
      ocpInputs :: OcpPhaseInputs' (IntegrationOcp x p)
      ocpInputs =
        OcpPhaseInputs
        { ocpPathCBnds = None
        , ocpBcBnds =  fmap (\x -> (Just x, Just x)) x0
        , ocpXbnd = pure (Nothing, Nothing)
        , ocpUbnd = None
        , ocpZbnd = None
        , ocpPbnd = fmap (\x -> (Just x, Just x)) p
        , ocpTbnd = (Just tf, Just tf)
        , ocpFixedP = None
        }
  let guess :: J (CollTraj x None None p n deg) (Vector Double)
      guess = cat $ makeGuessSim (collocationRoots dirCollOpts) tf x0 (\_ x _ -> ode x p 0) (\_ _ -> None) p
  cp  <- makeCollProblem dirCollOpts ocp ocpInputs guess :: IO (CollProblem x None None p x None x None None None None None n deg)
  (_, eopt) <- solveNlp "test_ocp_integrator" solver (cpNlp cp) Nothing
  return $ case eopt of
    Left m -> Left m
    Right opt -> Right (toXf (xOpt opt))
runIntegration _ _ dirCollOpts ode x0 p tf = do
  let ocp :: OcpPhase' (IntegrationOcp x p)
      ocp =
        OcpPhase
        { ocpMayer = \_ _ _ _ _ _ -> 0
        , ocpLagrange = \_ _ _ _ _ _ _ _ -> 0
        , ocpDae = \x' x _ _ pp _ t -> ((ode x pp t) `minus` x', None)
        , ocpQuadratures = \_ _ _ _ _ _ _ _ -> None
        , ocpQuadratureOutputs = \_ _ _ _ _ _ _ _ -> None
        , ocpBc = \x0' _ _ _ _ _ -> x0'
        , ocpPathC = \_ _ _ _ _ _ _ -> None
        , ocpPlotOutputs = \_ _ _ _ _ _ _ _ _ _ _ -> None
        , ocpObjScale      = Nothing
        , ocpTScale        = Nothing
        , ocpXScale        = Nothing
        , ocpZScale        = Nothing
        , ocpUScale        = Nothing
        , ocpPScale        = Nothing
        , ocpResidualScale = Nothing
        , ocpBcScale       = Nothing
        , ocpPathCScale    = Nothing
        }
      ocpInputs :: OcpPhaseInputs' (IntegrationOcp x p)
      ocpInputs =
        OcpPhaseInputs
        { ocpPathCBnds = None
        , ocpBcBnds =  fmap (\x -> (Just x, Just x)) x0
        , ocpXbnd = pure (Nothing, Nothing)
        , ocpUbnd = None
        , ocpZbnd = None
        , ocpPbnd = fmap (\x -> (Just x, Just x)) p
        , ocpTbnd = (Just tf, Just tf)
        , ocpFixedP = None
        }
  let guess :: J (CollTraj x None None p n deg) (Vector Double)
      guess = cat $ makeGuessSim (collocationRoots dirCollOpts) tf x0 (\_ x _ -> ode x p 0) (\_ _ -> None) p
  cp  <- makeCollProblem dirCollOpts ocp ocpInputs guess :: IO (CollProblem x None None p x None x None None None None None n deg)
  (_, eopt) <- solveNlp "test_ocp_integrator" solver (cpNlp cp) Nothing
  return $ case eopt of
    Left m -> Left m
    Right opt -> Right (toXf (xOpt opt))