text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -4.847398418942631 <= q <= 3.5461678827819525)

                     add_NL_constraint(m, :((-0.19745953189550702 + -0.698583212677272*(-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.07925279857906231*(-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.8658353105986021*(-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((-0.19745953189550702 + -0.698583212677272*(-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.07925279857906231*(-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.8658353105986021*(-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2 + (0.9047173820312868 + 0.09375520961098927*(-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.9282166334997988*(-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.3928411969784107*(-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.9047173820312868 + 0.09375520961098927*(-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.9282166334997988*(-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.3928411969784107*(-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2 + (0.16193188153327842 + 0.4964335766578394*(-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.6165233713312808*(-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.4911670927879337*(-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.16193188153327842 + 0.4964335766578394*(-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.2024495957518173 + -0.1501654471281335*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.7965619538518913*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.8656395810315112*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.6165233713312808*(-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.4774715210476401 + -0.18966514690787006*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + 0.1687755600264067*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.9882336248352352*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.4911670927879337*(-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)*(1 + erf((-0.5967554189497073 + 0.7632066186026765*(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))*(1 + erf((-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2]))/sqrt(2)))/2 + -0.767873838532454*(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))*(1 + erf((-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2]))/sqrt(2)))/2 + 0.02410968604346264*(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))*(1 + erf((-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Autogenerated wrapper script for MMseqs2_jll for powerpc64le-linux-gnu-cxx03
export mmseqs

using CompilerSupportLibraries_jll
using Zlib_jll
using Bzip2_jll
JLLWrappers.@generate_wrapper_header("MMseqs2")
JLLWrappers.@declare_executable_product(mmseqs)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, Zlib_jll, Bzip2_jll)
    JLLWrappers.@init_executable_product(
        mmseqs,
        "bin/mmseqs",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
"""
    SplineApproximations

Module for function approximation using splines.

The general idea is to find the spline ``g(x)`` in a given spline space that
best approximates a known function ``f(x)``.
In other words, given a predefined [`BSplineBasis`](@ref), the objective is to
find some appropriate B-spline coefficients such that the resulting spline
appropriately approximates ``f``.
Different approximation approaches are implemented, trading accuracy and speed.
"""
module SplineApproximations

using ..BSplines
using ..SplineInterpolations
using ..Splines

# For documentation
using ..Galerkin: galerkin_matrix
using ..Recombinations: RecombinedBSplineBasis

export
    VariationDiminishing,
    ApproxByInterpolation,
    MinimiseL2Error,
    approximate,
    approximate!

"""
    AbstractApproxMethod

Abstract type describing a type of approximation method.
"""
abstract type AbstractApproxMethod end

struct SplineApproximation{
        ApproxSpline <: Spline,
        Method <: AbstractApproxMethod,
        Data,
    } <: SplineWrapper{ApproxSpline}
    method :: Method
    spline :: ApproxSpline  # approximating spline
    data   :: Data
end

method(a::SplineApproximation) = a.method
SplineInterpolations.spline(a::SplineApproximation) = a.spline
data(a::SplineApproximation) = a.data

function Base.show(io::IO, A::SplineApproximation)
    print(io, nameof(typeof(A)), " containing the ", spline(A))
    print(io, "\n approximation method: ", method(A))
    nothing
end

"""
    approximate(f, B::AbstractBSplineBasis, [method = ApproxByInterpolation(B)])

Approximate function `f` in the given basis, using the chosen method.

From lower to higher accuracy (and cost), the possible approximation methods are:

- [`VariationDiminishing`](@ref),
- [`ApproxByInterpolation`](@ref),
- [`MinimiseL2Error`](@ref).

See their respective documentations for details.

Note that, once an approximation has been performed, it's possible to
efficiently perform additional approximations (of other functions `f`) by
calling the in-place [`interpolate!`](@ref).
This completely avoids allocations and strongly reduces computation time.

# Examples

```jldoctest
julia> B = BSplineBasis(BSplineOrder(3), -1:0.4:1);

julia> S_interp = approximate(sin, B)
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [-0.841471, -0.731727, -0.39727, 2.85767e-17, 0.39727, 0.731727, 0.841471]
 approximation method: interpolation at [-1.0, -0.8, -0.4, 0.0, 0.4, 0.8, 1.0]

julia> sin(0.3), S_interp(0.3)
(0.29552020666133955, 0.2959895327282942)

julia> approximate!(exp, S_interp)
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.367879, 0.440373, 0.65701, 0.980127, 1.46223, 2.18111, 2.71828]
 approximation method: interpolation at [-1.0, -0.8, -0.4, 0.0, 0.4, 0.8, 1.0]

julia> exp(0.3), S_interp(0.3)
(1.3498588075760032, 1.3491015490105396)

julia> S_fast = approximate(exp, B, VariationDiminishing())
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.367879, 0.449329, 0.67032, 1.0, 1.49182, 2.22554, 2.71828]
 approximation method: VariationDiminishing()

julia> S_opt = approximate(exp, B, MinimiseL2Error())
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.368074, 0.440342, 0.657077, 0.980279, 1.46216, 2.18201, 2.71669]
 approximation method: MinimiseL2Error()

julia> x = 0.34; exp(x), S_opt(x), S_interp(x), S_fast(x)
(1.4049475905635938, 1.4044530324752076, 1.4044149581073813, 1.4328668494041878)
```
"""
approximate(f, B::AbstractBSplineBasis) =
    approximate(f, B, ApproxByInterpolation(B))

"""
    approximate!(f, A::SplineApproximation)

Approximate function `f` reusing a previous Spline approximation in a given
basis.

See [`approximate`](@ref) for details.
"""
approximate!(f, A::SplineApproximation) = _approximate!(f, A, method(A))

include("variation_diminishing.jl")
include("by_interpolation.jl")
include("minimiseL2.jl")

end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
module TestDict

using Referenceables
using Test

@testset begin
    x = Dict()
    y = referenceable(x)
    r = y[1]
    r[] = 10
    @test x[1] == 10
end

end  # module


text:
module UnitfulEquivalences

export @eqrelation, Equivalence, MassEnergy, Spectral, Thermal

import Unitful
using Unitful: AbstractQuantity, AffineQuantity, DimensionlessQuantity, Dimensions, Level,
               NoDims, NoUnits, Quantity, Units, absoluteunit, dimension, uconvert

"""
    Equivalence

Abstract supertype for all equivalences.
"""
abstract type Equivalence end

Base.broadcastable(x::Equivalence) = Ref(x)

"""
    edconvert(d::Dimensions, x::AbstractQuantity, e::Equivalence)

Convert `x` to the equivalent dimension `d` using the equivalence `e`. (not exported)

# Example

```jldoctest
julia> using UnitfulEquivalences: edconvert

julia> edconvert(dimension(u"J"), 1u"kg", MassEnergy()) # E = m*c^2
89875517873681764 kg m^2 s^-2
```
"""
edconvert(d::Dimensions, x::AbstractQuantity, e::Equivalence) =
    throw(ArgumentError("$e defines no equivalence between dimensions $(dimension(x)) and $d."))

_scalarquantity(x::AbstractQuantity) = x
_scalarquantity(x::AffineQuantity)   = uconvert(absoluteunit(x), x)
_scalarquantity(x::Number)           = Quantity{typeof(x),NoDims,typeof(NoUnits)}(x)

"""
    uconvert(u::Units, x::Quantity, e::Equivalence)

Convert `x` to the units `u` (of different dimensions) by using the specified equivalence.

# Examples

```jldoctest
julia> uconvert(u"keV", 1u"me", MassEnergy()) # electron rest mass is equivalent to ≈511 keV
510.9989499961642 keV

julia> uconvert(u"eV", 589u"nm", Spectral()) # photon energy of sodium D₂ line (≈589 nm)
2.104994880020378 eV
```
"""
Unitful.uconvert(u::Units, x, e::Equivalence) =
    uconvert(u, edconvert(dimension(u), _scalarquantity(x), e))

"""
    ustrip([T::Type,] u::Units, x::Quantity, e::Equivalence)

Convert `x` to the units `u` (of different dimensions) by using the specified equivalence
and return the numeric value of the resulting quantity. If `T` is supplied, also convert the
resulting number to type `T`.

# Examples

```jldoctest
julia> ustrip(u"keV", 1u"me", MassEnergy()) # electron rest mass is equivalent to ≈511 keV
510.9989499961642

julia> ustrip(u"eV", 589u"nm", Spectral()) # photon energy (in eV) of sodium D₂ line
2.104994880020378
```
"""
Unitful.ustrip(u::Units, x, e::Equivalence) =
    ustrip(u, edconvert(dimension(u), _scalarquantity(x), e))
Unitful.ustrip(T::Type, u::Units, x, e::Equivalence) =
    ustrip(T, u, edconvert(dimension(u), _scalarquantity(x), e))

"""
    dimtype(x)

For a quantity type alias as created by `Unitful.@dimension` or `Unitful.@derived_dimension`
(e.g., `Unitful.Energy`), return its `Dimensions` type. (not exported)

# Example

```jldoctest
julia> using UnitfulEquivalences: dimtype

julia> dimtype(Unitful.Length)
Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}
```
"""
dimtype(::Type{Union{Quantity{T,D,U}, Level{L,S,Quantity{T,D,U}} where {L,S}} where {T,U}}) where D = typeof(D)
dimtype(::typeof(DimensionlessQuantity)) = typeof(NoDims)

"""
    @eqrelation Name a/b = c
    @eqrelation Name a*b = c

Add a proportional or antiproportional relation between dimensions `a` and `b` to an
existing equivalence type `Name <: Equivalence`. The dimensions `a` and `b` must be
specified as quantity type aliases like `Unitful.Energy`.

# Example

```@julia
using Unitful: Energy, Mass, c0
struct MassEnergy <: Equivalence end
@eqrelation MassEnergy Energy/Mass = c0^2
```
In the rest frame of a particle, its energy is proportional to its mass. Defining the
`MassEnergy` equivalence like above allows conversion between energies and masses via
`uconvert(massunit, energy, MassEnergy())` and
`uconvert(energyunit, massunit, MassEnergy())`.
"""
macro eqrelation(name, relation)
    relation isa Expr && relation.head == :(=) || _eqrelation_error()
    lhs, rhs = relation.args
    lhs isa Expr && lhs.head == :call && length(lhs.args) == 3 || _eqrelation_error()
    op, a, b = lhs.args
    if op == :/
        quote
            UnitfulEquivalences.edconvert(::dimtype($(esc(a))), x::$(esc(b)), ::$(esc(name))) = x * $(esc(rhs))
            UnitfulEquivalences.edconvert(::dimtype($(esc(b))), x::$(esc(a)), ::$(esc(name))) = x / $(esc(rhs))
            nothing
        end
    elseif op == :*
        quote
            UnitfulEquivalences.edconvert(::dimtype($(esc(a))), x::$(esc(b)), ::$(esc(name))) = $(esc(rhs)) / x
            UnitfulEquivalences.edconvert(::dimtype($(esc(b))), x::$(esc(a)), ::$(esc(name))) = $(esc(rhs)) / x
            nothing
        end
    else
        _eqrelation_error()
    end
end

_eqrelation_error() = error("second macro argument must be an (anti-)proportionality relation " *
                            "`a/b = c` or `a*b = c`, cf. the documentation for `@eqrelation`.")

include("pkgdefaults.jl")

end # module UnitfulEquivalences


text:
"""
```
pseudo_measurement(m::Model1002{T}, TTT::Matrix{T}, RRR::Matrix{T},
    CCC::Vector{T}; reg::Int = 1) where {T<:AbstractFloat}
```

Assign pseudo-measurement equation (a linear combination of states):

```
x_t = ZZ_pseudo*s_t + DD_pseudo
```
"""
function pseudo_measurement(m::Model1002{T},
                            TTT::Matrix{T},
                            RRR::Matrix{T},
                            CCC::Vector{T};
                            reg::Int = 1,
                            TTTs::Vector{<: AbstractMatrix{T}} = Matrix{T}[],
                            CCCs::Vector{<: AbstractVector{T}} = Vector{T}[],
                            information_set::UnitRange = reg:reg,
                            memo::Union{ForwardMultipleExpectationsMemo, Nothing} = nothing) where {T <: AbstractFloat}

    endo      = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    pseudo    = m.pseudo_observables

    # Initialize pseudo ZZ and DD matrices
    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)

    # Set parameters
    for para in m.parameters
        if !isempty(para.regimes)
            if (haskey(get_settings(m), :model2para_regime) ? haskey(get_setting(m, :model2para_regime), para.key) : false)
                ModelConstructors.toggle_regime!(para, reg, get_setting(m, :model2para_regime)[para.key])
            else
                ModelConstructors.toggle_regime!(para, reg)
            end
        end
    end

    # Set up for calculating k-periods ahead expectations and expected sums
    permanent_t = length(information_set[findfirst(x -> x == reg, information_set):end]) - 1 + reg
    if information_set[1] == information_set[end]
        # In this case, we do not need to pass the TTTs, CCCs in, so we redefine them as empty.
        # This step is also necessary to ensure the memo is properly used.
        TTTs = Matrix{T}[]
        CCCs = Vector{T}[]

        memo = nothing # see measurement
        # TODO: maybe instead of emptying TTTs, CCCs, we add
        # a step to recompute the memo since we will still likely be recalculating
        # products/powers of TTT multiple times that could be pre-computed
    end
    use_fwd_exp_sum = haskey(get_settings(m), :use_forward_expected_sum_memo) && get_setting(m, :use_forward_expected_sum_memo)
    use_fwd_exp     = haskey(get_settings(m), :use_forward_expectations_memo) && get_setting(m, :use_forward_expectations_memo)

    # Handle integrated series
    no_integ_inds = inds_states_no_integ_series(m)
    if ((haskey(get_settings(m), :add_altpolicy_pgap) && haskey(get_settings(m), :pgap_type)) ?
        (get_setting(m, :add_altpolicy_pgap) && get_setting(m, :pgap_type) == :ngdp) : false) ||
        (haskey(get_settings(m), :ait_Thalf) ? (exp(log(0.5) / get_setting(m, :ait_Thalf)) ≈ 0.) : false)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:pgap_t]])
    end
    if (haskey(get_settings(m), :gdp_Thalf) ? (exp(log(0.5) / get_setting(m, :gdp_Thalf)) ≈ 0.) : false)
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:ygap_t]])
    end
    if haskey(get_settings(m), :ρ_rw) ? (get_setting(m, :ρ_rw) ≈ 1.) : false
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:rw_t]])
    end
    if haskey(get_settings(m), :rw_ρ_smooth) ? (get_setting(m, :rw_ρ_smooth) ≈ 1.) : false
        no_integ_inds = setdiff(no_integ_inds, [m.endogenous_states[:Rref_t]])
    end
    integ_series = length(no_integ_inds) != n_states_augmented(m) # Are the series integrated?

    # Compute TTT^10, used for Expected10YearRateGap, Expected10YearRate, and Expected10YearNaturalRate
    TTT10, CCC10 = k_periods_ahead_expected_sums(TTT, CCC, TTTs, CCCs, reg, 40, permanent_t;
                                                 integ_series = integ_series,
                                                 memo = use_fwd_exp_sum ? memo : nothing)
    TTT10        = TTT10./ 40. # divide by 40 to average across 10 years
    CCC10        = CCC10 ./ 40.

    if get_setting(m, :add_laborproductivity_measurement)
        # Construct pseudo-obs from integrated states first
        ZZ_pseudo[pseudo[:laborproductivity], endo[:y_t]] = 1.
        ZZ_pseudo[pseudo[:laborproductivity], endo[:L_t]] = -1.
        # ZZ_pseudo[pseudo[:laborproductivity], endo_addl[:cum_z_t]] = 1.
        DD_pseudo[pseudo[:laborproductivity]] = 100. * log(m[:ystar] / m[:Lstar])

        # Remove integrated states (e.g. states w/unit roots)
        # RRR and CCC aren't used, so we don't do anything with them
    end

    if get_setting(m, :add_nominalgdp_level)
        ZZ_pseudo[pseudo[:NominalGDPLevel], endo_addl[:cum_y_t]]     = 1.
        ZZ_pseudo[pseudo[:NominalGDPLevel], endo_addl[:cum_z_t]]     = 1.
        ZZ_pseudo[pseudo[:NominalGDPLevel], endo_addl[:cum_e_gdp_t]] = 1.
        ZZ_pseudo[pseudo[:NominalGDPLevel], endo_addl[:cum_π_t]]     = 1.
    end

    if get_setting(m, :add_nominalgdp_growth)
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo[:y_t]]           = 1.
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo_addl[:y_t1]]     = -1.
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo_addl[:z_t]]      = 1.
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo_addl[:e_gdp_t]]  = 1.
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo_addl[:e_gdp_t1]] = -m[:me_level]
        ZZ_pseudo[pseudo[:NominalGDPGrowth], endo[:π_t]]           = 1.
        DD_pseudo[pseudo[:NominalGDPGrowth]]                       = 100. * (exp(m[:z_star] - 1.) + (m[:π_star] - 1.))
    end

    if get_setting(m, :add_cumulative)
        ZZ_pseudo[pseudo[:AccumOutputGap], endo_addl[:cum_y_t]]   = 1.
        ZZ_pseudo[pseudo[:AccumOutputGap], endo_addl[:cum_y_f_t]] = -1.

        ZZ_pseudo[pseudo[:GDPLevel], endo_addl[:cum_y_t]]     = 1.
        ZZ_pseudo[pseudo[:GDPLevel], endo_addl[:cum_z_t]]     = 1.
        ZZ_pseudo[pseudo[:GDPLevel], endo_addl[:cum_e_gdp_t]] = 1.

        ZZ_pseudo[pseudo[:FlexibleGDPLevel], endo_addl[:cum_y_f_t]] = 1.
        ZZ_pseudo[pseudo[:FlexibleGDPLevel], endo_addl[:cum_z_t]]   = 1.

        ZZ_pseudo[pseudo[:ConsumptionLevel], endo_addl[:cum_c_t]]     = 1.
        ZZ_pseudo[pseudo[:ConsumptionLevel], endo_addl[:cum_z_t]]     = 1.

        ZZ_pseudo[pseudo[:FlexibleConsumptionLevel], endo_addl[:cum_c_f_t]] = 1.
        ZZ_pseudo[pseudo[:FlexibleConsumptionLevel], endo_addl[:cum_z_t]]   = 1.

        ZZ_pseudo[pseudo[:InvestmentLevel], endo_addl[:cum_i_t]]     = 1.
        ZZ_pseudo[pseudo[:InvestmentLevel], endo_addl[:cum_z_t]]     = 1.

        ZZ_pseudo[pseudo[:FlexibleInvestmentLevel], endo_addl[:cum_i_f_t]] = 1.
        ZZ_pseudo[pseudo[:FlexibleInvestmentLevel], endo_addl[:cum_z_t]]   = 1.
    end

    if get_setting(m, :add_flexible_price_growth)
        ZZ_pseudo[pseudo[:FlexibleGDPGrowth], endo[:y_f_t]]       = 1.
        ZZ_pseudo[pseudo[:FlexibleGDPGrowth], endo_addl[:y_f_t1]] = -1.
        ZZ_pseudo[pseudo[:FlexibleGDPGrowth], endo[:z_t]]         = 1.
        DD_pseudo[pseudo[:FlexibleGDPGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)

        ZZ_pseudo[pseudo[:FlexibleInvestmentGrowth], endo[:i_f_t]]       = 1.
        ZZ_pseudo[pseudo[:FlexibleInvestmentGrowth], endo_addl[:i_f_t1]] = -1.
        ZZ_pseudo[pseudo[:FlexibleInvestmentGrowth], endo[:z_t]]         = 1.
        DD_pseudo[pseudo[:FlexibleInvestmentGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)

        ZZ_pseudo[pseudo[:FlexibleConsumptionGrowth], endo[:c_f_t]]       = 1.
        ZZ_pseudo[pseudo[:FlexibleConsumptionGrowth], endo_addl[:c_f_t1]] = -1.
        ZZ_pseudo[pseudo[:FlexibleConsumptionGrowth], endo[:z_t]]         = 1.
        DD_pseudo[pseudo[:FlexibleConsumptionGrowth]]                     = 100. * (exp(m[:z_star]) - 1.)
    end

    ##########################################################
    ## PSEUDO-OBSERVABLE EQUATIONS
    ##########################################################

    ## Output
    ZZ_pseudo[pseudo[:y_t],endo[:y_t]] = 1.

    ## Flexible Output
    ZZ_pseudo[pseudo[:y_f_t],endo[:y_f_t]] = 1.

    ## Pseudo GDP Growth
    if haskey(m.settings, :add_pseudo_gdp)
        if get_setting(m, :add_pseudo_gdp) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            ZZ_pseudo[pseudo[:PseudoGDP], endo[:y_t]]          = 1.0
            ZZ_pseudo[pseudo[:PseudoGDP], endo_addl[:y_t1]]     = -1.0
            ZZ_pseudo[pseudo[:PseudoGDP], endo[:z_t]]          = 1.0
            ZZ_pseudo[pseudo[:PseudoGDP], endo_addl[:e_gdp_t]]  = 1.0
            ZZ_pseudo[pseudo[:PseudoGDP], endo_addl[:e_gdp_t1]] = -m[:me_level]
            DD_pseudo[pseudo[:PseudoGDP]]                      = 100*(exp(m[:z_star])-1)
        end
    end

    ## Pseudo Core PCE # TODO
    if haskey(m.settings, :add_pseudo_corepce)
        if get_setting(m, :add_pseudo_corepce) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            ZZ_pseudo[pseudo[:PseudoCorePCE], endo[:π_t]]              = 1.0
            ZZ_pseudo[pseudo[:PseudoCorePCE], endo_addl[:e_corepce_t]] = 1.0
            DD_pseudo[pseudo[:PseudoCorePCE]]                          = 100. * (m[:π_star] - 1.)
        end
    end

    ## Natural Rate
    ZZ_pseudo[pseudo[:NaturalRate],endo[:r_f_t]] = 1.
    DD_pseudo[pseudo[:NaturalRate]]              = 100.0*(m[:rstar]-1.0)

    ## π_t
    ZZ_pseudo[pseudo[:π_t],endo[:π_t]] = 1.
    DD_pseudo[pseudo[:π_t]]            = 100*(m[:π_star]-1);

    ## Output Gap
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_t]] = 1;
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_f_t]] = -1;

    ## Ex Ante Real Rate
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:R_t]]  = 1;
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:Eπ_t]] = -1;
    DD_pseudo[pseudo[:ExAnteRealRate]]             = m[:Rstarn] - 100*(m[:π_star]-1);

    ## Long Run Inflation
    ZZ_pseudo[pseudo[:LongRunInflation],endo[:π_star_t]] = 1.
    DD_pseudo[pseudo[:LongRunInflation]]                 = 100. *(m[:π_star]-1.)

    ## Marginal Cost
    ZZ_pseudo[pseudo[:MarginalCost],endo[:mc_t]] = 1.

    ## Wages
    ZZ_pseudo[pseudo[:Wages],endo[:w_t]] = 1.

    ## Flexible Wages
    ZZ_pseudo[pseudo[:FlexibleWages],endo[:w_f_t]] = 1.

    # ## b Wages
    # ZZ_pseudo[pseudo[:b_t],endo[:b_t]] = 1.

    ## Hours
    ZZ_pseudo[pseudo[:Hours],endo[:L_t]] = 1.

    ## Flexible Hours
    ZZ_pseudo[pseudo[:FlexibleHours],endo[:L_f_t]] = 1.

    ## z_t
    ZZ_pseudo[pseudo[:z_t], endo[:z_t]] = 1.

    ## Expected 10-Year Rate Gap
    # ZZ_pseudo[pseudo[:Expected10YearRateGap], :] = TTT10[endo[:R_t], :] - TTT10[endo[:r_f_t], :] - TTT10[endo[:Eπ_t], :]
    ZZ_pseudo[pseudo[:Expected10YearRateGap], :] = view(TTT10, endo[:R_t], :) - view(TTT10, endo[:r_f_t], :) - view(TTT10, endo[:Eπ_t], :)
    DD_pseudo[pseudo[:Expected10YearRateGap]]    = CCC10[endo[:R_t]] - CCC10[endo[:r_f_t]] - CCC10[endo[:Eπ_t]]

    ## Nominal FFR
    ZZ_pseudo[pseudo[:NominalFFR], endo[:R_t]] = 1.
    DD_pseudo[pseudo[:NominalFFR]] = m[:Rstarn]

    ## Expected 10-Year Interest Rate
    # ZZ_pseudo[pseudo[:Expected10YearRate], :] = TTT10[endo[:R_t], :]
    ZZ_pseudo[pseudo[:Expected10YearRate], :] = view(TTT10, endo[:R_t], :)
    DD_pseudo[pseudo[:Expected10YearRate]]    = m[:Rstarn] + CCC10[endo[:R_t]]

    ## Expected 10-Year Natural Rate
    # ZZ_pseudo[pseudo[:Expected10YearNaturalRate], :] = TTT10[endo[:r_f_t], :] + TTT10[endo[:Eπ_t], :]
    ZZ_pseudo[pseudo[:Expected10YearNaturalRate], :] = view(TTT10, endo[:r_f_t], :) + view(TTT10, endo[:Eπ_t], :)
    DD_pseudo[pseudo[:Expected10YearNaturalRate]]    = m[:Rstarn] + CCC10[endo[:r_f_t]] + CCC10[endo[:Eπ_t]]

    ## Expected Nominal Natural Rate
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:r_f_t]] = 1.
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:Eπ_t]]  = 1.
    DD_pseudo[pseudo[:ExpectedNominalNaturalRate]]               = m[:Rstarn]

    ## Nominal Rate Gap
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:R_t]]   = 1.
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:r_f_t]] = -1.
    ZZ_pseudo[pseudo[:NominalRateGap], endo[:Eπ_t]]  = -1.

    ## Labor Productivity Growth
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:y_t]]           = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:y_t1]]     = -1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:z_t]]           = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t]]  = 1.
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:e_gdp_t1]] = -m[:me_level]
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo[:L_t]]           = -1
    ZZ_pseudo[pseudo[:LaborProductivityGrowth], endo_addl[:L_t1]]     = 1.
    DD_pseudo[pseudo[:LaborProductivityGrowth]]                       = 100*(exp(m[:z_star]) - 1)

    ## u_t
    ZZ_pseudo[pseudo[:u_t], endo[:u_t]] = 1.

    ## Nominal Wage Growth
    if haskey(m.settings, :add_NominalWageGrowth)
        if get_setting(m, :add_NominalWageGrowth)
            ZZ_pseudo[pseudo[:NominalWageGrowth],endo[:w_t]] = 1.
            ZZ_pseudo[pseudo[:NominalWageGrowth],endo_addl[:w_t1]] = -1.
            ZZ_pseudo[pseudo[:NominalWageGrowth],endo[:z_t]] = 1.
            ZZ_pseudo[pseudo[:NominalWageGrowth],endo[:π_t]] = 1.
            DD_pseudo[pseudo[:NominalWageGrowth]]            = 100*(m[:π_star]-1) + 100*(exp(m[:z_star])-1)
        end
    end

    # ## i_f_t
    # ZZ_pseudo[pseudo[:i_f_t], endo[:i_f_t]] = 1.

    # ## R_t
    # ZZ_pseudo[pseudo[:R_t], endo[:R_t]] = 1.
    # DD_pseudo[pseudo[:R_t]] = 100.0*(m[:rstar]-1.0)

    # ## c_f_t
    # ZZ_pseudo[pseudo[:c_f_t], endo[:c_f_t]] = 1.

    # ## c_t
    # ZZ_pseudo[pseudo[:c_t], endo[:c_t]] = 1.

    # ## qk_f_t
    # ZZ_pseudo[pseudo[:qk_f_t], endo[:qk_f_t]] = 1.

    # ## k_f_t
    # ZZ_pseudo[pseudo[:k_f_t], endo[:k_f_t]] = 1.

    # ## r_f_t
    # ZZ_pseudo[pseudo[:r_f_t], endo[:r_f_t]] = 1.

    # ## kbar_f_t
    # ZZ_pseudo[pseudo[:kbar_f_t], endo[:kbar_f_t]] = 1.

    # ## u_f_t
    # ZZ_pseudo[pseudo[:u_f_t], endo[:u_f_t]] = 1.

    # ## rk_f_t
    # ZZ_pseudo[pseudo[:rk_f_t], endo[:rk_f_t]] = 1.

    # ## w_f_t
    # ZZ_pseudo[pseudo[:w_f_t], endo[:w_f_t]] = 1.

    # ## L_f_t
    # ZZ_pseudo[pseudo[:L_f_t], endo[:L_f_t]] = 1.

    # ## rktil_f_t
    # ZZ_pseudo[pseudo[:rktil_f_t], endo[:rktil_f_t]] = 1.

    # ## n_f_t
    # ZZ_pseudo[pseudo[:n_f_t], endo[:n_f_t]] = 1.

    ## labor share
    if haskey(m.settings, :add_laborshare_measurement)
        if get_setting(m, :add_laborshare_measurement)
            ZZ_pseudo[pseudo[:laborshare_t], endo[:w_t]] = 1.
            ZZ_pseudo[pseudo[:laborshare_t], endo[:L_t]] = 1.
            ZZ_pseudo[pseudo[:laborshare_t], endo[:y_t]] = -1.
            DD_pseudo[pseudo[:laborshare_t]] = 100. * log(m[:wstar] * m[:Lstar] / m[:ystar])
        end
    end

    ## Fundamental inflation related pseudo-obs
    if subspec(m) in ["ss13", "ss14", "ss15", "ss16", "ss17", "ss18", "ss19"]
        # Compute coefficient on Sinf
        betabar = exp((1-m[:σ_c] ) * m[:z_star]) * m[:β]
        κ = ((1 - m[:ζ_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*
             (1 - m[:ζ_p]))/(m[:ζ_p]*((m[:Φ]- 1)*m[:ϵ_p] + 1))/
        (1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
        κcoef = κ * (1 + m[:ι_p] * betabar)

        ZZ_pseudo[pseudo[:Sinf_t], endo_addl[:Sinf_t]] = 1.
        ZZ_pseudo[pseudo[:Sinf_w_coef_t], endo_addl[:Sinf_t]] = κcoef
        DD_pseudo[pseudo[:ι_p]] = m[:ι_p]
        ZZ_pseudo[pseudo[:πtil_t], endo_addl[:πtil_t]] = 1.
        DD_pseudo[pseudo[:πtil_t]] = 100 * (m[:π_star] - 1)
        ZZ_pseudo[pseudo[:e_tfp_t], endo_addl[:e_tfp_t]] = 1.
        if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
            ZZ_pseudo[pseudo[:e_tfp_t1], endo_addl[:e_tfp_t1]] = 1.
        end
    end

    ## Exogenous processes
    if subspec(m) == "ss12"
        to_add = [:g_t, :b_t, :μ_t, :z_t, :λ_f_t, :λ_w_t, :rm_t, :σ_ω_t, :μ_e_t,
                  :γ_t, :π_star_t]
        to_add_addl = [:e_lr_t, :e_tfp_t, :e_gdpdef_t, :e_corepce_t, :e_gdp_t, :e_gdi_t]
        for i in to_add
            ZZ_pseudo[pseudo[i], endo[i]] = 1.
        end
        for i in to_add_addl
            ZZ_pseudo[pseudo[i], endo_addl[i]] = 1.
        end
    end

    if haskey(m.settings, :add_covid_pseudoobs)
        if get_setting(m, :add_covid_pseudoobs) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            ZZ_pseudo[pseudo[:ziid], endo[:ziid_t]] = 1.
            ZZ_pseudo[pseudo[:varphiiid], endo[:φ_t]] = 1.
            ZZ_pseudo[pseudo[:biidc], endo[:biidc_t]] = 1.
        end
    end

    if haskey(m.settings, :add_ztil)
        if get_setting(m, :add_ztil)
            ZZ_pseudo[pseudo[:ztil], endo[:ztil_t]] = 1.
        end
    end

    if haskey(m.settings, :add_zp)
        if get_setting(m, :add_zp)
            ZZ_pseudo[pseudo[:zp], endo[:zp_t]] = 1.
        end
    end

    if haskey(m.settings, :add_pgap)
        if get_setting(m, :add_pgap)
            ZZ_pseudo[pseudo[:pgap], endo[:pgap_t]] = 1.
        end
    end

    if haskey(m.settings, :add_ygap)
        if get_setting(m, :add_ygap)
            ZZ_pseudo[pseudo[:ygap], endo[:ygap_t]] = 1.
        end
    end

    if haskey(m.settings, :add_altpolicy_pgap)
        if get_setting(m, :add_altpolicy_pgap)
            ZZ_pseudo[pseudo[:pgap], endo[:pgap_t]] = 1.
        end
    end

    if haskey(m.settings, :add_altpolicy_ygap)
        if get_setting(m, :add_altpolicy_ygap)
            ZZ_pseudo[pseudo[:ygap], endo[:ygap_t]] = 1.
        end
    end

    if haskey(m.settings, :add_urhat)
        if get_setting(m, :add_urhat)
            ZZ_pseudo[pseudo[:urhat], endo[:L_t]] = -0.7366220966444673
            DD_pseudo[pseudo[:urhat]] = -32.257825316004364 + -0.7366220966444673 * m[:Lmean]
        end
    end

    if haskey(m.settings, :add_rw)
        if get_setting(m, :add_rw)
            ZZ_pseudo[pseudo[:rw], endo[:rw_t]]     = 1.
            ZZ_pseudo[pseudo[:Rref], endo[:Rref_t]] = 1.
            DD_pseudo[pseudo[:Rref]]                = m[:Rstarn]
        end
    end

    for para in m.parameters
        if !isempty(para.regimes)
            ModelConstructors.toggle_regime!(para, 1)
        end
    end

    return PseudoMeasurement(ZZ_pseudo, DD_pseudo)
end


text:
Random.seed!(1234)

## METADATA

plain_classifier  = SVMClassifier()
nu_classifier     = SVMNuClassifier()
linear_classifier = SVMLClassifier(max_iter=10000)

plain_regressor  = SVMRegressor()
nu_regressor     = SVMNuRegressor()
linear_regressor = SVMLRegressor(max_iter=10_000)

i = info_dict(plain_classifier)

@testset "SVM-infos" begin
    @test i[:name]          == "SVMClassifier"
    @test i[:package_name]  == "ScikitLearn"
    @test i[:package_url]   == "https://github.com/cstjean/ScikitLearn.jl"
    @test i[:load_path]     == "MLJModels.ScikitLearn_.SVMClassifier"
    @test i[:is_wrapper]    == false
    @test i[:is_pure_julia] == false
    @test i[:is_supervised] == true

    for m in (nu_classifier, linear_classifier, plain_regressor, nu_regressor,
              linear_regressor)
        info_dict(m)
    end
end

## CLASSIFIERS
@test_logs (:warn,"kernel parameter is not valid, setting to default=\"rbf\" \n") SVMClassifier(kernel="wrong")
@test_logs (:warn,"penalty parameter is not valid, setting to default=\"l2\" \n") SVMLClassifier(penalty="wrong")
@test_logs (:warn,"loss parameter is not valid, setting to default=\"epsilon_insensitive\" \n") SVMLRegressor(loss="wrong")

# test preservation of categorical levels:

iris = dataset("datasets", "iris")

X = iris[:, 1:4]
y = iris[:, 5]

train, test = partition(eachindex(y), 0.6) # levels of y are split across split

fitresultC, cacheC, reportC = fit(plain_classifier, 1,
                                          selectrows(X, train), y[train])
fitresultCnu, cacheCnu, reportCnu = fit(nu_classifier, 1,
                                          selectrows(X, train), y[train])
fitresultCL, cacheCL, reportCL = fit(linear_classifier, 1,
                                          selectrows(X, train), y[train])
pcpred  = predict(plain_classifier, fitresultC, selectrows(X, test))
nucpred = predict(nu_classifier, fitresultCnu, selectrows(X, test))
lcpred  = predict(linear_classifier, fitresultCL, selectrows(X, test))

@test Set(classes(pcpred[1])) == Set(classes(y[1]))
@test Set(classes(nucpred[1])) == Set(classes(y[1]))
@test Set(classes(lcpred[1])) == Set(classes(y[1]))

# test with linear data:
x1 = randn(3_000)
x2 = randn(3_000)
x3 = randn(3_000)
X  = (x1=x1, x2=x2, x3=x3)
y  = x1 - x2 -2x3

ycat = map(y) do η
    η > 0 ? "go" : "stop"
end |> categorical

@testset "SVCs" begin
    train, test = partition(eachindex(ycat), 0.8)
    fitresultC, cacheC, reportC = fit(plain_classifier, 1,
                                              selectrows(X, train), ycat[train])
    fitresultCnu, cacheCnu, reportCnu = fit(nu_classifier, 1,
                                              selectrows(X, train), ycat[train])
    fitresultCL, cacheCL, reportCL = fit(linear_classifier, 1,
                                              selectrows(X, train), ycat[train])
    pcpred  = predict(plain_classifier, fitresultC, selectrows(X, test))
    nucpred = predict(nu_classifier, fitresultCnu, selectrows(X, test))
    lcpred  = predict(linear_classifier, fitresultCL, selectrows(X, test))

    @test sum(pcpred  .!= ycat[test])/length(ycat) < 0.05
    @test sum(nucpred .!= ycat[test])/length(ycat) < 0.05
    @test sum(lcpred  .!= ycat[test])/length(ycat) < 0.05
end

## REGRESSORS

# test with linear data:
fitresultR, cacheR, reportR = fit(plain_regressor, 1,
                                          selectrows(X, train), y[train])
fitresultRnu, cacheRnu, reportRnu = fit(nu_regressor, 1,
                                                selectrows(X, train), y[train])
fitresultRL, cacheRL, reportRL = fit(linear_regressor, 1,
                                             selectrows(X, train), y[train])

@testset "SVRs" begin
    rpred   = predict(plain_regressor, fitresultR, selectrows(X, test))
    nurpred = predict(nu_regressor, fitresultRnu, selectrows(X, test))
    Lrpred  = predict(linear_regressor, fitresultRL, selectrows(X, test))

    @test norm(rpred   - y[test])/sqrt(length(y)) < 0.2
    @test norm(nurpred - y[test])/sqrt(length(y)) < 0.2
    @test norm(Lrpred  - y[test])/sqrt(length(y)) < 0.2
end


text:
struct Cycle{T, C}
    complex   ::C
    indices   ::Vector{Int}
    perimeter ::T
    diameter  ::T
end

Cycle(complex, indices) =
    Cycle(complex, indices, perimeter_and_diameter(complex, indices)...)

Base.show(io::IO, α::Cycle) =
    print(io, "$(length(α))-element Cycle with p = $(perimeter(α)), d = $(diameter(α))")

Base.length(α::Cycle) = length(α.indices)

points(α::Cycle) = points(α.complex, α.indices)

perimeter(α::Cycle) = α.perimeter

LightGraphs.diameter(α::Cycle) = α.diameter

"""
    getindexmatrix(dists)

Create index matrix `M`, such that `M[i, j] = k` if the value `dists[i, j]` is the `k`-th
largest. Assumes `dists` is symmetric.
"""
function getindexmatrix(dists)
    M = fill(0, size(dists))
    for (i, (u, v)) in enumerate(sortededges(dists))
        M[u, v] = i
        M[v, u] = i
    end
    M
end

"""
Internal state for intrinsic persistence.
"""
struct PersistenceState{T}
    # Precomputed properties.
    dists       ::Matrix{T}
    parents     ::Matrix{Int}
    indexmatrix ::Matrix{Int}
    triangles   ::EquilateralIterator{T}
    radius      ::T
    # Internal state
    reduced     ::Vector{BitSet}
    # Buffers
    σ           ::BitSet
    cycle       ::Vector{Int}
    # Results
    results     ::Vector{Tuple{Vector{Int}, NTuple{3, Int}, T}}
end

function PersistenceState(gc)
    r = radius(gc)
    # todo: split connected components
    dists, parents = landmark_shortest_paths(gc)
    indexmatrix = getindexmatrix(dists)
    triangles = equilaterals(dists, 2r)

    reduced = fill(BitSet(), binomial(n_landmarks(gc), 2))

    PersistenceState(dists, parents, indexmatrix, triangles, r, reduced,
                     BitSet(), Int[], Tuple{Vector{Int}, NTuple{3, Int}, eltype(dists)}[])
end

Base.show(io::IO, st::PersistenceState) = print(io, "PersistenceState")

"""
    isgeodesic(st, Δ)

Check if `st.cycle` spanned by triangle `Δ` is a geodesic circle.
"""
function isgeodesic(st, Δ)
    length(st.cycle) == 3 && return true
    d = st.dists
    dik = djk = zero(eltype(d))
    # add offset to get it right in the first go.
    j, k, i = Δ
    for x in st.cycle
        if x ∈ Δ
            # make sure x is always between i and j.
            i, j, k = j, k, i
            dik = d[i, k]
            djk = d[j, k]
            continue
        end
        # add st.radius/2 to a void numerical errors.
        # st.radius is the shortest edge in graph.
        if d[x, k] + st.radius/2 < min(d[x, i] + dik,
                                       d[x, j] + djk)
            return false
        end
    end
    true
end

"""
    low(σ::BitSet)

Get the largest index in `σ` / the lowest entry in the column.
"""
low(σ::BitSet) = length(σ) > 0 ? last(σ) : 0

"""
    reduce!(st::PersistenceState)

Use the standard reduction algorithm to reduce `st.σ`.
"""
function reduce!(st::PersistenceState)
    l = low(st.σ)
    while l ≠ 0 && !isempty(st.reduced[l])
        symdiff!(st.σ, st.reduced[l])
        l = low(st.σ)
    end
    l
end

"""
    appendpath!(st::PersistenceState, u, v)

Append the shortest path from `u` to `v` into `st.cycle` and `st.σ`.
"""
function appendpath!(st::PersistenceState, u, v)
    while u ≠ v
        u′ = st.parents[v, u]
        push!(st.σ, st.indexmatrix[u, u′])
        push!(st.cycle, u)
        u = u′
    end
end

"""
    movecycle!(st::PersistenceState, Δ)

Move the cycle spanned by triangle `Δ` into `st.cycle` and `st.σ`.
"""
function movecycle!(st::PersistenceState, (i, j, k))
    empty!(st.σ)
    empty!(st.cycle)
    appendpath!(st, i, j)
    appendpath!(st, j, k)
    appendpath!(st, k, i)
end

"""
TODO
"""
function processtriangle!(st::PersistenceState, triangle, diameter)
    movecycle!(st, triangle)
    # Skip cycles that visit a node more than once and nongeodesic cycles.
    length(st.σ) == length(st.cycle) || return
    isgeodesic(st, triangle) || return

    l = reduce!(st)
    if l ≠ 0
        st.reduced[l] = copy(st.σ)
        push!(st.results, (copy(st.cycle), triangle, diameter))
    end
end

function persistence(gc; showprogress = false, dense = false, shrinked = true)
    showprogress && println("Calculating intrinsic persistence...")
    st = PersistenceState(gc)

    for (Δ, diam) in st.triangles
        processtriangle!(st, Δ, diam)
    end
    showprogress && println("Postprocessing...")
    cycles = map(st.results) do (α, _, _)
        landmarks(gc, α)
    end
    postprocess(gc, cycles, dense = dense, shrinked = shrinked)
end

# ======================================================================================== #

"""
find the index and distance of point in `gc` within `gc.radius` of point(s) `p` that
minimizes the maximum distance to points in `pts`.
"""
function centerpoint(gc, p, pts)
    inball = nearby_points(gc, p)
    dists = pairwise_ambient_distance(gc, inball, pts)
    μ, i = findmin(vec(maximum(dists, dims = 2)))
    μ, inball[i]
end

shrink(gc, α) = shrink!(gc, copy(α))
shrink(α::Cycle) = shrink(α.complex, α.indices)

function shrink!(gc, α)
    changed = true
    while changed
        changed = false
        unique!(α)
        n = length(α)
        for (i, v) in enumerate(α)
            j = mod1(i - 1, n)
            k = mod1(i + 1, n)
            _, new = centerpoint(gc, v, [α[j], α[k]])
            if new ≠ v
                changed = true
                α[i] = new
            end
        end
    end
    α
end
shrink!(α::Cycle) = shrink!(α.complex, α.indices)

function densify(gc::GeodesicComplex, α)
    α = copy(α)
    β = Int[]
    changed = true
    while changed
        changed = false
        for i in eachindex(α)
            u = α[i]
            v = α[mod1(i + 1, length(α))]
            push!(β, u)
            inballs = nearby_points(gc, [u, v])
            candidates = union(inballs...)
            dists = pairwise_ambient_distance(gc, candidates, [u, v])
            new = candidates[argmin(vec(maximum(dists, dims = 2)))]
            if new ≠ u && new ≠ v
                changed = true
                push!(β, new)
            end
        end
        move!(α, β)
        shrink!(gc, α)
    end
    α
end
densify(α::Cycle) = densify(α.complex, α.indices)

function move!(a, b)
    resize!(a, length(b))
    copyto!(a, b)
    empty!(b)
    a
end

function postprocess(gc, cycles; dense = false, shrinked = true)
    res = map(cycles) do α
        if shrinked
            α = shrink(gc, α)
        end
        if dense
            α = densify(gc, α)
        end
        Cycle(gc, α)
    end
    filter!(res) do α
        length(α) > 2 && diameter(α) > 2radius(gc)
    end
end

function perimeter_and_diameter(gc, α)
    dists = pairwise_ambient_distance(gc, α, α)

    perimeter = sum(eachindex(α)) do i
        dists[i, mod1(i+1, length(α))]
    end
    diameter = maximum(dists)

    perimeter, diameter
end

max_stars_repo_name: mtsch/IntrinsicPersistence

max_issues_repo_name: mtsch/IntrinsicPersistence

max_forks_repo_name: mtsch/IntrinsicPersistence


text:
function dynamics_derivatives!(data::ModelData; mode = :nominal)
    if mode == :nominal
        x̄ = data.x̄
        ū = data.ū
    else
        x̄ = data.x
        ū = data.u
    end

    w = data.w
    h = data.h
    T = data.T
    model = data.model

    for t = 1:T-1
       	if data.analytical_derivatives
			data.dyn_deriv.fx[t] .= fdx(model, x̄[t], ū[t], w[t], h, t)
	      	data.dyn_deriv.fu[t] .= fdu(model, x̄[t], ū[t], w[t], h, t)
		else
			fx(z) = fd(model, z, ū[t], w[t], h, t)
	        fu(z) = fd(model, x̄[t], z, w[t], h, t)
	        # fw(z) = fd(model, x̄[t], ū[t], z, h, t)

	        # data.dyn_deriv.fx[t] .= ForwardDiff.jacobian(fx, x̄[t])
	        # data.dyn_deriv.fu[t] .= ForwardDiff.jacobian(fu, ū[t])
	        # # data.dyn_deriv.fw[t] = ForwardDiff.jacobian(fw, w[t])

			ForwardDiff.jacobian!(data.dyn_deriv.fx[t], fx, x̄[t])
	        ForwardDiff.jacobian!(data.dyn_deriv.fu[t], fu, ū[t])
	        # ForwardDiff.jacobian!(data.dyn_deriv.fw[t], fw, w[t])
		end
    end
end

function objective_derivatives!(obj::StageCosts, data::ModelData;
    mode = :nominal)

    if mode == :nominal
        x̄ = data.x̄
        ū = data.ū
    else
        x̄ = data.x
        ū = data.u
    end

    T = data.T
    model = data.model
    n = data.n
    m = data.m

    for t = 1:T-1
		# println("time step $t")
		if obj.cost[t] isa QuadraticCost
			data.obj_deriv.gx[t] .= 2.0 * obj.cost[t].Q * x̄[t] + obj.cost[t].q
			data.obj_deriv.gu[t] .= 2.0 * obj.cost[t].R * ū[t] + obj.cost[t].r
			data.obj_deriv.gxx[t] .= 2.0 * obj.cost[t].Q
			data.obj_deriv.guu[t] .= 2.0 * obj.cost[t].R
			data.obj_deriv.gux[t] .= 0.0
		else
	        gx(z) = g(obj, z, ū[t], t)
	        gu(z) = g(obj, x̄[t], z, t)
	        gz(z) = g(obj, z[1:n[t]], z[n[t] .+ (1:m[t])], t)

	        ForwardDiff.gradient!(data.obj_deriv.gx[t], gx, x̄[t])
	        ForwardDiff.gradient!(data.obj_deriv.gu[t], gu, ū[t])
			ForwardDiff.hessian!(data.obj_deriv.gxx[t], gx, x̄[t])
	        ForwardDiff.hessian!(data.obj_deriv.guu[t], gu, ū[t])
	        data.obj_deriv.gux[t] .= ForwardDiff.hessian(gz,
	            [x̄[t]; ū[t]])[n[t] .+ (1:m[t]), 1:n[t]]
		end
    end

	if obj.cost[T] isa QuadraticCost
		data.obj_deriv.gx[T] .= 2.0 * obj.cost[T].Q * x̄[T] + obj.cost[T].q
		data.obj_deriv.gxx[T] .= 2.0 * obj.cost[T].Q
	else
	    gxT(z) = g(obj, z, nothing, T)
	    ForwardDiff.gradient!(data.obj_deriv.gx[T], gxT, x̄[T])
	    ForwardDiff.hessian!(data.obj_deriv.gxx[T], gxT, x̄[T])
	end
end

function constraints_derivatives!(cons::StageConstraints, data::ModelData;
    mode = :nominal)

    if mode == :nominal
        x̄ = data.x̄
        ū = data.ū
    else
        x̄ = data.x
        ū = data.u
    end

    T = data.T

    for t = 1:T-1
        c = cons.data.c[t]
        cx!(a, z) = c!(a, cons, z, ū[t], t)
        cu!(a, z) = c!(a, cons, x̄[t], z, t)

        ForwardDiff.jacobian!(cons.data.cx[t], cx!, c, x̄[t])
        ForwardDiff.jacobian!(cons.data.cu[t], cu!, c, ū[t])
    end

    c = cons.data.c[T]
    cxT!(a, z) = c!(a, cons, z, nothing, T)
    ForwardDiff.jacobian!(cons.data.cx[T], cxT!, c, x̄[T])
end

function objective_derivatives!(obj::AugmentedLagrangianCosts, data::ModelData;
        mode = :nominal)

    gx = data.obj_deriv.gx
    gu = data.obj_deriv.gu
    gxx = data.obj_deriv.gxx
    guu = data.obj_deriv.guu
    gux = data.obj_deriv.gux

    c = obj.cons.data.c
    cx = obj.cons.data.cx
    cu = obj.cons.data.cu
    ρ = obj.ρ
    λ = obj.λ
    a = obj.a

    T = data.T
    model = data.model

    objective_derivatives!(obj.costs, data, mode = mode)
    constraints_derivatives!(obj.cons, data, mode = mode)

    for t = 1:T-1
        gx[t] .+= cx[t]' * (λ[t] + ρ[t] .* a[t] .* c[t])
        gu[t] .+= cu[t]' * (λ[t] + ρ[t] .* a[t] .* c[t])
        gxx[t] .+= cx[t]' * Diagonal(ρ[t] .* a[t]) * cx[t]
        guu[t] .+= cu[t]' * Diagonal(ρ[t] .* a[t]) * cu[t]
        gux[t] .+= cu[t]' * Diagonal(ρ[t] .* a[t]) * cx[t]
    end

    gx[T] .+= cx[T]' * (λ[T] + ρ[T] .* a[T] .* c[T])
    gxx[T] .+= cx[T]' * Diagonal(ρ[T] .* a[T]) * cx[T]
end

function derivatives!(m_data::ModelData; mode = :nominal)
    dynamics_derivatives!(m_data, mode = mode)
    objective_derivatives!(m_data.obj, m_data, mode = mode)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Documenter, REopt, JuMP

makedocs(
    sitename="REopt.jl Documentation",
    pages = [
        "Home" => "index.md",
        "REopt" => Any[
            "reopt/examples.md",
            "reopt/inputs.md",
            "reopt/outputs.md",
            "reopt/methods.md"
        ],
        "Model Predictive Control" => Any[
            "mpc/examples.md",
            "mpc/inputs.md",
            "mpc/methods.md",
        ],
        "Developer" => Any[
            "developer/concept.md",
            "developer/organization.md",
            "developer/inputs.md",
            "developer/adding_tech.md"
        ]
    ],
    workdir = joinpath(@__DIR__, "..")
)

deploydocs(
    repo = "github.com/NREL/REopt.jl.git",
)

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
using DimensionalPlotRecipes, Test

A = rand(5,2) .+ im.*rand(5,2)
t = range(0,stop=1,length=5)

using Plots
plot(t,A)

t = range(0,stop=1,length=5)
plot(t,A)
plot(t,A,transformation=:split3D)
plot(t,A,transformation=:modulus)
plot(t,A,transformation=:modulus2)


text:
function create_region(obj::ConvexHull)
    return jcall(obj, "createRegion", Region, ())
end

function get_vertices(obj::ConvexHull)
    return jcall(obj, "getVertices", Vector{Point}, ())
end


text:
x = [1.0, 2.0, missing, missing]
y = [3.0, 3.0, 3.0, 3.0]
f(x,y) = x+y
function g(x,y)
    z = f.(x,y)
    if z isa Union
        return typeof(z)
        convert(typeof(z).b.a, z)
    else
        z
    end
end
z = g(x,y)
@code_warntype g(x,y)

typeof(x), typeof(x[1:2]), typeof(x[3:4]) # preserves eltype
typeof(x.*1.0), typeof(x[1:2].*1.0), typeof(x[3:4].*1.0) # does not
typeof(x*1.0), typeof(x[1:2]*1.0), typeof(x[3:4]*1.0) # does not

x = [1.0, 2.0, missing, missing]
y = 1.0
typeof(x.*y), typeof(x[1:2].*y), typeof(x[3:4].*y)
# (Vector{Union{Missing, Float64}}, Vector{Float64}, Vector{Missing})
res = x[1:2].*y
res_eltype = Union{eltype(x),typeof(y)}
convert(Vector{res_eltype}, res)::Vector{res_eltype}

f1(x,y) = x[1:2].*y
@code_warntype f1(x,y)

function f2(x,y)
    res_eltype = Union{eltype(x),typeof(y)}
    convert(Vector{res_eltype}, x[1:2].*y)::Vector{res_eltype}
end
@inferred f2(x,y)


text:
function calc_pairwise_sums(array)
    sums = Dict{Tuple{Int64, Int64}, Int64}()
    for (ix, x) in enumerate(array)
        for (jy, y) in enumerate(array[ix+1:end])
            sums[(ix, jy + ix)] = x + y
        end
    end
    sums
end

function find_encoding_error(numbers, k)
    preamble = numbers[1:k]
    sums = calc_pairwise_sums(preamble)
    for (ix, elem) in enumerate(numbers[k+1:end])
        if elem ∉ values(sums)
            return elem
        end

        filter!((((i, _), _),) -> i != ix, sums)
        for (jy, y) in enumerate(preamble[2:end])
            sums[(ix + jy, ix + k)] = elem + y
        end

        push!(preamble, elem)
        deleteat!(preamble, 1)
    end
end

function calc_encryption_weakness(numbers, n)
    start = 1
    stop = 1
    s = numbers[1]
    while s != n
        if s < n
            stop += 1
            s += numbers[stop]
        elseif s > n
            s -= numbers[start]
            start += 1
        end
    end

    contiguous_array = numbers[start:stop]
    maximum(contiguous_array) + minimum(contiguous_array)
end

function main()
    numbers = map(x -> parse(Int64, x), readlines("../resources/input_09.txt"))
    k = 25
    n = find_encoding_error(numbers, k)
    result = calc_encryption_weakness(numbers, n)
    println(result)
end

main()


text:
# exact k-medoids modeled similar as in Kotzur et al, 2017

"Holds results of kmedoids run"
mutable struct kmedoidsResult
    medoids::Array{Float64}
    assignments::Array{Int}
    totalcost::Float64
end

"""
    kmedoids_exact(
     data::Array{Float64},
     nclust::Int,
     _dist::SemiMetric = SqEuclidean(),
     env::Any;
     )
   results = kmedoids_exact()
   data { HOURS,DAYS }
Performs the exact kmedoids algorithm as in Kotzur et al, 2017

env = Gurobi.Env() # reusing the same gurobi environment for multiple solves
"""
function kmedoids_exact(
   data::Array{Float64},
   nclust::Int,
   env::Any;
   _dist::SemiMetric = SqEuclidean(),
   )
N_i = size(data,2)

# calculate distance matrix
d_mat=pairwise(_dist,data; dims=2)

# create jump model
m = Model(optimizer=GurobiSolver(env)) # GurobiSolver(env,OutputFlag=0)
@variable(m,z[1:N_i,1:N_i],Bin)
@variable(m,y[1:N_i],Bin)
@objective(m,Min,sum(d_mat[i,j]*z[i,j] for i=1:N_i, j=1:N_i))
for j=1:N_i
  @constraint(m,sum(z[i,j] for i=1:N_i)==1)
end
for i=1:N_i
  for j=1:N_i
  @constraint(m,z[i,j]<=y[i])
  end
end
@constraint(m,sum(y[i] for i=1:N_i) == nclust)

# solve jump model
#tic()
status=solve(m)
#toc()

println("status: ",status)
y_opt=round.(Integer,getvalue(y))
z_opt=round.(Integer,getvalue(z))
 #println("y ",y_opt, " z ",z_opt)
# determine centers and cluster mappings
id = zeros(Int,N_i)
ii=0
for i=1:N_i
  if y_opt[i]==1
    ii +=1
    id[i]=ii
  end
end
nz = findall(!iszero,z_opt)
centerids=Int[]
for i=1:length(nz) # just take the first dimension of each of the elements of the array of cartesianCoordinates
  push!(centerids, nz[i][1])
end
clustids = zeros(Int,N_i)
for i=1:N_i
  clustids[i] = id[centerids[i]]
end
centers = data[:,findall(id.!=0.0)]
tot_dist = getobjectivevalue(m)
 # output centers
results = kmedoidsResult(centers, clustids, tot_dist)

return results

end #kmedoids_exact


text:
@doc """ index array data with the first column
""" ->
function index(dt)
    dict = Dict{ASCIIString,Array{ASCIIString,1}}()
    for i = 1:size(dt,1)
        dict[dt[i,1]] = dt[i,2:end]
    end
    dict
end

@doc """ left join data on the first column, assume have the same ids
         And data have no header part
""" ->
function leftjoindata(d1,d2)
    d1_ind = index(d1)
    d2_ind = index(d2)
    @assert issubset(Set(collect(keys(d1_ind))),Set(collect(keys(d2_ind)))) == true
    nrow = size(d1,1)
    ncol = size(d1,2)-1 + size(d2,2)-1 + 1
    data = Array{ASCIIString,2}(nrow,ncol)
    for (i,id) in enumerate(keys(d1_ind))
        data[i,:] = vcat(id,d1_ind[id], d2_ind[id])'
    end
    data
end

function leftjoindata(d1,d2,d3...)
    #left join
    data = leftjoindata(d1,d2)
    for d in d3
        data = leftjoindata(data,d)
    end
    data
end

@doc """ the last column as lable column, ratio is 6:2:2
         assume data doest not contain header part
""" ->
function splitdata{T<:Any}(dt::Array{T,2})
    labels = Set(dt[:,end])
    num_row,num_col = size(dt)
    tr,val,te = Array{T,2}(0,num_col),Array{T,2}(0,num_col),Array{T,2}(0,num_col)

    for lb in labels
        inds = dt[:,end] .== lb
        data = dt[inds,:]
        len = size(data,1)
        trval = Int64(round(len*0.6))
        valte = trval + Int64(round(len*0.2))
        tr  = vcat(tr, data[1:trval,:])
        val = vcat(val,data[trval+1:valte,:])
        te  = vcat(te, data[valte+1:end,:])
    end

    tr,val,te
end

@doc """ shuffle dataset
random indexes of row
""" ->
function shuffle(data)
    m,n = size(data)
    inds = randperm(m)
    data[inds,:]
end

function makedata()
    ev = readcsv("../data/event_type.csv", ASCIIString)
    re = readcsv("../data/resource_type.csv", ASCIIString)
    se = readcsv("../data/severity_type.csv", ASCIIString)
    log = readcsv("../data/log_feature.csv", ASCIIString)
    ev = reformat(ev)
    re = reformat(re)
    se = reformat(se)
    log = reformat3(log)

    tr = readcsv("../data/train.csv", ASCIIString)
    te = readcsv("../data/test.csv", ASCIIString)
    tr_idloc = reformat(tr[:, [1,2]])
    tr_idlb  = tr[:, [1,3]]
    te_idloc = reformat(te[:, [1,2]])

    train = leftjoindata(tr_idloc,ev,re,se,log,tr_idlb)
    test = leftjoindata(te_idloc,ev,re,se,log)

    tr,val,te = map(shuffle,  splitdata(train))
    tr  = tr[:,2:end]
    val = val[:,2:end]
    te  = te[:,2:end]
    JLD.save(joinpath(data_dir,"tr.jld"),"tr", tr)
    JLD.save(joinpath(data_dir,"val.jld"),"val", val)
    JLD.save(joinpath(data_dir,"te.jld"),"te", te)
    JLD.save(joinpath(data_dir,"test.jld"),"test", test)
end


text:
module AxisIndices

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end AxisIndices

using AbstractFFTs
using ArrayInterface
using ChainedFixes
using IntervalSets
using LinearAlgebra
using MappedArrays
using Metadata
using NamedDims
using SparseArrays
using StaticRanges
using Statistics
using SuiteSparse

using EllipsisNotation: Ellipsis
using StaticRanges
using StaticRanges: OneToUnion
using StaticRanges: can_set_first, can_set_last, can_set_length, same_type
using StaticRanges: checkindexlo, checkindexhi
using StaticRanges: grow_first!, grow_last!, grow_first, grow_last
using StaticRanges: shrink_first!, shrink_last!, shrink_first, shrink_last
using StaticRanges: resize_last, resize_last!, resize_first, resize_first!
using StaticRanges: is_static, is_fixed, similar_type

using Base: @propagate_inbounds, tail, LogicalIndex, Slice, OneTo, Fix2, ReinterpretArray
using Base: IdentityUnitRange, setindex
using Base.Broadcast: Broadcasted, BroadcastStyle, DefaultArrayStyle, AbstractArrayStyle, Unknown
using ArrayInterface: OptionallyStaticUnitRange
using ArrayInterface: known_length, known_first, known_step, known_last, can_change_size
using ArrayInterface: static_length, static_first, static_step, static_last, StaticInt, Zero, One
using ArrayInterface: indices, offsets, to_index, to_axis, to_axes, unsafe_reconstruct, parent_type
using ArrayInterface: can_setindex
using MappedArrays: MultiMappedArray, ReadonlyMultiMappedArray
using Base: print_array, print_matrix_row, print_matrix, alignment

const to_indices = ArrayInterface.to_indices

export
    AbstractAxis,
    AxisArray,
    AxisMatrix,
    AxisVector,
    Axis,
    AxisArray,
    CartesianAxes,
    LinearAxes,
    NamedAxisArray,
    SimpleAxis,
    StructAxis,
    center,
    closest,
    circular_pad,
    idaxis,
    permuteddimsview,
    offset,
    one_pad,
    reflect_pad,
    replicate_pad,
    struct_view,
    symmetric_pad,
    zero_pad

const ArrayInitializer = Union{UndefInitializer, Missing, Nothing}

# Val wraps the number of axes to retain
naxes(A::AbstractArray, v::Val) = naxes(axes(A), v)
naxes(axs::Tuple, v::Val{N}) where {N} = _naxes(axs, N)
@inline function _naxes(axs::Tuple, i::Int)
    if i === 0
        return ()
    else
        return (first(axs), _naxes(tail(axs), i - 1)...)
    end
end

@inline function _naxes(axs::Tuple{}, i::Int)
    if i === 0
        return ()
    else
        return (SimpleAxis(1), _naxes((), i - 1)...)
    end
end

include("errors.jl")
include("abstract_axis.jl")
include("axis_array.jl")

"""
    AxisInitializer <: Function

Supertype for functions that assist in initialization of `AbstractAxis` subtypes.
"""
abstract type AxisInitializer <: Function end

(init::AxisInitializer)(x) = Base.Fix2(init, x)
function (init::AxisInitializer)(collection::AbstractRange, x)
    if known_step(collection) === 1
        return axis_method(init, x, collection)
    else
        return AxisArray(collection, ntuple(_ -> init(x), Val(ndims(collection))))
    end
end
function (init::AxisInitializer)(collection, x)
    return AxisArray(collection, ntuple(_ -> init(x), Val(ndims(collection))))
end
function (init::AxisInitializer)(collection::AbstractRange, x::Tuple)
    if ndims(collection) !== length(x)
        throw(DimensionMismatch("Number of axis arguments provided ($(length(x))) does " *
                                "not match number of collections's axes ($(ndims(collection)))."))
    end
    if known_step(collection) === 1
        return axis_method(init, first(x), collection)
    else
        return AxisArray(collection, map(init, x))
    end
end
function (init::AxisInitializer)(collection, x::Tuple)
    if ndims(collection) !== length(x)
        throw(DimensionMismatch("Number of axis arguments provided ($(length(x))) does " *
                                "not match number of collections's axes ($(ndims(collection)))."))
    end
    return AxisArray(collection, map(init, x))
end

include("simple_axis.jl")
include("axis.jl")
include("offset_axis.jl")
include("centered_axis.jl")
include("identity_axis.jl")
include("padded_axis.jl")
include("struct_axis.jl")

# TODO assign_indices tests
function assign_indices(axis, inds)
    if can_change_size(axis) && !((known_length(inds) === nothing) || known_length(inds) === known_length(axis))
        return unsafe_reconstruct(axis, inds)
    else
        return axis
    end
end

"""
    is_key([collection,] arg) -> Bool

Whether `arg` refers to a key of `axis`.
"""
is_key(arg) = is_key(IndexLinear(), typeof(arg))
is_key(collection, arg) = is_key(IndexStyle(collection), typeof(arg))
is_key(collection, ::Type{I}) where {I} = is_key(IndexStyle(collection), I)
is_key(::IndexStyle, ::Type{T}) where {T<:Colon} = false
is_key(::IndexStyle, ::Type{T}) where {T<:Integer} = false
is_key(S::IndexStyle, ::Type{T}) where {T<:AbstractArray} = is_key(S, eltype(T))
is_key(::IndexStyle, ::Type{T}) where {T} = true

const MetaAxisArray{T,N,P,Axs,M} = Metadata.MetaArray{T,N,AxisArray{T,N,P,Axs},M}

const NamedMetaAxisArray{L,T,N,P,M,Axs} = NamedDimsArray{L,T,N,MetaAxisArray{T,N,P,Axs,M}}

include("to_index.jl")
include("checkindex.jl")
include("getindex.jl")
include("permutedims.jl")
include("axes_methods.jl")
include("combine.jl")
include("arrays.jl")
include("alias_arrays.jl")
include("linear_algebra.jl")
include("deprecations.jl")
include("abstractarray.jl")
include("promotion.jl")
include("named.jl")
include("show.jl")
include("fft.jl")

# TODO move this to ArrayInterface
ArrayInterface._multi_check_index(axs::Tuple, arg::LogicalIndex{<:Any,<:AxisArray}) = axs == axes(arg.mask)

###
### offsets
###
@inline function apply_offset(axis, arg)
    if arg isa Integer
        return Int(arg)
    else
        return arg
    end
end
apply_offset(axis::OffsetAxis, arg) = _apply_offset(getfield(axis, :offset), arg)
apply_offset(axis::IdentityAxis, arg) = _apply_offset(getfield(axis, :offset), arg)
function apply_offset(axis::CenteredAxis, arg)
    p = parent(axis)
    return _apply_offset(_origin_to_offset(first(p), length(p), origin(axis)), arg)
end
_apply_offset(f, arg::Integer) = arg - f
_apply_offset(f, arg::AbstractArray) = arg .- f
function _apply_offset(f, arg::AbstractRange)
    if known_step(arg) === 1
        return (first(arg) - f):(last(arg) - f)
    else
        return (first(arg) - f):step(arg):(last(arg) - f)
    end
end

# add offsets
_add_offset(axis, x) = x
_add_offset(axis::OffsetAxis, arg) = __add_offset(getfield(axis, :offset), arg)
_add_offset(axis::IdentityAxis, arg) = __add_offset(getfield(axis, :offset), arg)
function _add_offset(axis::CenteredAxis, arg)
    p = parent(axis)
    return __add_offset(_origin_to_offset(first(p), length(p), origin(axis)), arg)
end
__add_offset(f, arg::Integer) = arg + f
__add_offset(f, arg::AbstractArray) = arg .+ f
function __add_offset(f, arg::AbstractRange)
    if known_step(arg) === 1
        return (first(arg) + f):(last(arg) + f)
    else
        return (first(arg) + f):step(arg):(last(arg) + f)
    end
end

# subtract offsets
apply_offsets(::Tuple{}, ::Tuple{}) = ()
apply_offsets(::Tuple{}, ::Tuple) = ()
apply_offsets(::Tuple, ::Tuple{}) = ()
@inline function apply_offsets(axs::Tuple{A}, inds::Tuple{<:Integer}) where {A}
    return (_sub_offset(first(axs), first(inds)),)
end
@inline function apply_offsets(axs::Tuple, inds::Tuple)
    return (_sub_offset(first(axs), first(inds)), apply_offsets(tail(axs), tail(inds))...)
end
_sub_offset(axis, x) = x
_sub_offset(axis::OffsetAxis, arg) = __sub_offset(getfield(axis, :offset), arg)
_sub_offset(axis::IdentityAxis, arg) = __sub_offset(getfield(axis, :offset), arg)
function _sub_offset(axis::CenteredAxis, arg)
    p = parent(axis)
    return __sub_offset(_origin_to_offset(first(p), length(p), origin(axis)), arg)
end
__sub_offset(f, arg::Integer) = arg - f
__sub_offset(f, arg::AbstractArray) = arg .- f
function __sub_offset(f, arg::AbstractRange)
    if known_step(arg) === 1
        return (first(arg) - f):(last(arg) - f)
    else
        return (first(arg) - f):step(arg):(last(arg) - f)
    end
end

# Metadata stuff
@inline function Metadata.metadata(x::AxisArray; dim=nothing, kwargs...)
    if dim === nothing
        return metadata(parent(x); kwargs...)
    else
        return metadata(axes(x, dim); kwargs...)
    end
end

include("closest.jl")

export ..

end


text:
using Random

function makesim(sts::Union{STS{:CoupleHenon1D}, STS{:henon1d}})
    @unpack N, T, S, ic = sts
    @assert T isa Integer
    @assert S isa Integer
    # prepare initial condition...
    if ic isa Union{Nothing, Integer}
        rng = MersenneTwister(isnothing(ic) ? 27 : ic)
        u = 0.1rand(rng, N)
        v = 0.1rand(rng, N)
    elseif ic isa Tuple
        u, v = ic
    else
        error(ICERROR)
    end
    t, U, V = coupled_henon_1D(u, v, S, T; sts.p...)
    return @strdict t sts U V
end

"""
    coupled_henon_1D(u, v, S, T; a = 1.4, b = 0.3)
Coupled henon maps in one dimension.
"""
function coupled_henon_1D(u, v, S, T; a = 1.4, b = 0.3)
    @assert length(u) == length(v)
    N = length(u)
    u = copy(u); v = copy(v)

    function henon(U,V)
        Un = copy(U)
        Vn = copy(V)
        Un[1] = Un[N] = 0.5
        Vn[1] = Vn[N] = 0
        for m=2:N-1
            @inbounds Un[m] = 1 - a*(.5*U[m]+ .25*U[m-1] + .25*U[m+1])^2 + b*V[m]
            @inbounds Vn[m] = U[m]
        end
        return Un, Vn
    end

    for n in 1:S
        u, v = henon(u, v)
    end
    U = Vector{Vector{Float32}}(undef,T)
    V = Vector{Vector{Float32}}(undef,T)
    U[1] = u; V[1] = v
    for n = 2:T
        U[n], V[n] = henon(U[n-1],V[n-1])
    end
    return 1:T, U,V
end


text:
workspace()

include("../src/cpu/file_io/read_key.jl")
include("../src/cpu/file_io/read_input.jl")
include("../src/cpu/file_io/read_data.jl")
include("../src/cpu/file_io/write_key.jl")
include("../src/cpu/file_io/write_input.jl")
include("../src/cpu/file_io/write_data.jl")
include("../src/cpu/rand_weighted_index.jl")
include("../src/cpu/generate_experiment.jl")
include("../src/cpu/rand_poisson.jl")
include("../src/cpu/position_periodic.jl")

function test_part_conc_spt_abc()
	# Inititalization of random number generation device.
	random_seed::Int64 = 1#convert(Int64, time_ns())
	srand(random_seed)

	# Timing.
	time_start::Int64 = convert(Int64, time_ns())

	# Experimental parameters.
	ax::Float64 = 40.0 # µm.
	ay::Float64 = 40.0 # µm.
	Lx::Float64 = 60.0 # µm.
	Ly::Float64 = 60.0 # µm.
	Lz::Float64 = 10.0 # µm.
	number_of_frames::Array{Int64, 1} = 250 * ones(50)
	deltat::Float64 = 0.05 # seconds
	kmin::Int64 = 2
	model::String = "discrete-fixed-depth"
	m_real::Array{Float64, 1} = [1.0, 3.0] # µm^2/s.
	c_real::Array{Float64, 1} = [1e8, 1e8] # part/ml.
	az_real::Array{Float64, 1} = [2.0, 2.0] # µm.

	# Simulate experiment.
	(K::Array{Int64, 1}, DE::Array{Float64, 1}) = generate_experiment(	model,
																		m_real,
																		c_real,
																		ax,
																		ay,
																		az_real,
																		Lx,
																		Ly,
																		Lz,
																		number_of_frames,
																		deltat,
																		kmin)

	# Write data to file.
	data_file_path::String = abspath("data.xml")
	write_data(	data_file_path,
				ax,
				ay,
				number_of_frames,
				deltat,
				K,
				DE)

	# Write input to file.
	input_file_path::String = abspath("input.xml")
	number_of_components::Int64 = 2
	number_of_de_bins::Int64 = 200
	ub_de::Float64 = 4.0 * maximum(m_real)
	lb_m::Float64 = 0.25 * minimum(m_real)
	ub_m::Float64 = 4.0 * maximum(m_real)
	lb_c::Float64 = 0.25 * minimum(c_real)
	ub_c::Float64 = 4.0 * maximum(c_real)
	lb_az::Float64 = 0.25 * minimum(az_real)
	ub_az::Float64 = 4.0 * maximum(az_real)
	number_of_abc_samples::Int64 = 128
	gamma_initial::Float64 = 15.0
	gamma_adaptive::Bool = true#false
	delta_gamma::Float64 = 0.01
	#weighting_scheme::String = "pmc-standard"
	weighting_scheme::String = "inverse-distance-squared"
	ub_average_number_of_trials::Int64 = 500
	output_file_path::String = abspath("output.xml")
	write_input(input_file_path,
				data_file_path,
				model,
				number_of_components,
				Lx,
				Ly,
				Lz,
				kmin,
				number_of_de_bins,
				ub_de,
				lb_m,
				ub_m,
				lb_c,
				ub_c,
				lb_az,
				ub_az,
				number_of_abc_samples,
				gamma_initial,
				gamma_adaptive,
				delta_gamma,
				weighting_scheme,
				ub_average_number_of_trials,
				output_file_path)

	# Run inference.
	program_path::String = abspath("../src/cpu/run_part_conc_spt_abc.jl")
	#cmd::Cmd = `julia $program_path $input_file_path`
	str_number_of_cores::String = string(Sys.CPU_CORES)
	cmd::Cmd = `julia -p $str_number_of_cores $program_path $input_file_path`
	#cmd::Cmd = `julia -p 1 $program_path $input_file_path`
	#cmd::Cmd = `julia --math-mode=fast --check-bounds=no --optimize=3 -p $str_number_of_cores $program_path $input_file_path`
	run(cmd)

	time_end::Int64 = convert(Int64, time_ns())
	time_exec::Float64 = convert(Float64, time_end - time_start) / 1e9
	println(time_exec)
	nothing
end

test_part_conc_spt_abc()


text:
using BinaryProvider
using BinaryBuilder
using Base.Test
using SHA
using Compat

# The platform we're running on
const platform = platform_key()

# On windows, the `.exe` extension is very important
const exe_ext = is_windows() ? ".exe" : ""

# We are going to build/install libfoo a lot, so here's our function to make sure the
# library is working properly
function check_foo(fooifier_path = "fooifier$(exe_ext)",
                   libfoo_path = "libfoo.$(Libdl.dlext)")
    # We know that foo(a, b) returns 2*a^2 - b
    result = 2*2.2^2 - 1.1

    # Test that we can invoke fooifier
    @test !success(`$fooifier_path`)
    @test success(`$fooifier_path 1.5 2.0`)
    @test parse(Float64,readchomp(`$fooifier_path 2.2 1.1`)) ≈ result

    # Test that we can dlopen() libfoo and invoke it directly
    libfoo = Libdl.dlopen_e(libfoo_path)
    @test libfoo != C_NULL
    foo = Libdl.dlsym_e(libfoo, :foo)
    @test foo != C_NULL
    @test ccall(foo, Cdouble, (Cdouble, Cdouble), 2.2, 1.1) ≈ result
    Libdl.dlclose(libfoo)
end

@testset "File Collection" begin
    temp_prefix() do prefix
        # Create a file and a link, ensure that only the one file is returned by collect_files()
        f = joinpath(prefix, "foo")
        f_link = joinpath(prefix, "foo_link")
        touch(f)
        symlink(f, f_link)

        files = collect_files(prefix)
        @test length(files) == 1
        @test files[1] == abspath(f)
    end
end

# This file contains tests that require our cross-compilation environment
@testset "Builder Dependency" begin
    temp_prefix() do prefix
        # First, let's create a Dependency that just installs a file
        cd("build_tests") do
            # Our simple executable file, generated by bash
            test_exe_path = joinpath(bindir(prefix),"test_exe")
            test_exe = ExecutableProduct(test_exe_path)
            results = [test_exe]

            # Remember that we don't need to use `bash()` here because these commands
            # will be happening within the cross-compilation environment
            cmds = Cmd[]
            push!(cmds, `mkdir -p $(bindir(prefix))`)
            push!(cmds, `bash -c "printf '#!/bin/bash\necho test' > $(test_exe_path)"`)
            push!(cmds, `chmod 775 $(test_exe_path)`)

            dep = Dependency("bash_test", results, cmds, platform, prefix)

            @test build(dep; verbose=true)
            @test satisfied(dep)
            @test readstring(`$(test_exe_path)`) == "test\n"
        end

        # Next, build a simple shared library and an executable
        cd("build_tests/libfoo") do
            libfoo = LibraryProduct(prefix, "libfoo")
            fooifier = ExecutableProduct(prefix, "fooifier")
            steps = [`make clean`, `make install`]
            dep = Dependency("foo", [libfoo, fooifier], steps, platform, prefix)

            # Build it
            @test build(dep; verbose=true)
            @test satisfied(dep; verbose=true)

            # Test the binaries
            check_foo(locate(fooifier), locate(libfoo))

            # Also test the binaries through `activate()`
            activate(prefix)
            check_foo()
            deactivate(prefix)

            # Test that `collect_files()` works:
            all_files = collect_files(prefix)
            @test locate(libfoo) in all_files
            @test locate(fooifier) in all_files
        end
    end
end

@testset "Builder Packaging" begin
    # Clear out previous build products
    for f in readdir(".")
        if !endswith(f, ".tar.gz")
            continue
        end
        rm(f; force=true)
    end

    # Gotta set this guy up beforehand
    tarball_path = nothing
    tarball_hash = nothing

    temp_prefix() do prefix
        cd("build_tests/libfoo") do
            # First, build libfoo
            libfoo = LibraryProduct(prefix, "libfoo")
            fooifier = ExecutableProduct(prefix, "fooifier")
            steps = [`make clean`, `make install`]
            dep = Dependency("foo", [libfoo, fooifier], steps, platform, prefix)

            @test build(dep)
        end

        # Next, package it up as a .tar.gz file
        tarball_path, tarball_hash = package(prefix, "./libfoo"; verbose=true)
        @test isfile(tarball_path)
    end

    # Test that we can inspect the contents of the tarball
    contents = list_tarball_files(tarball_path)
    @test "bin/fooifier" in contents
    @test "lib/libfoo.$(Libdl.dlext)" in contents

    # Install it within a new Prefix
    temp_prefix() do prefix
        # Install the thing
        @test install(tarball_path, tarball_hash; prefix=prefix, verbose=true)

        # Ensure we can use it
        fooifier_path = joinpath(bindir(prefix), "fooifier")
        libfoo_path = joinpath(libdir(prefix), "libfoo.$(Libdl.dlext)")
        check_foo(fooifier_path, libfoo_path)
    end

    rm(tarball_path; force=true)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VpnClientConnectionHealthDetailListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{VpnClientConnectionHealthDetail} } # spec name: value

    function VpnClientConnectionHealthDetailListResult(;value=nothing)
        o = new()
        validate_property(VpnClientConnectionHealthDetailListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type VpnClientConnectionHealthDetailListResult

const _property_map_VpnClientConnectionHealthDetailListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_VpnClientConnectionHealthDetailListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{VpnClientConnectionHealthDetail}")
Base.propertynames(::Type{ VpnClientConnectionHealthDetailListResult }) = collect(keys(_property_map_VpnClientConnectionHealthDetailListResult))
Swagger.property_type(::Type{ VpnClientConnectionHealthDetailListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VpnClientConnectionHealthDetailListResult[name]))}
Swagger.field_name(::Type{ VpnClientConnectionHealthDetailListResult }, property_name::Symbol) =  _property_map_VpnClientConnectionHealthDetailListResult[property_name]

function check_required(o::VpnClientConnectionHealthDetailListResult)
    true
end

function validate_property(::Type{ VpnClientConnectionHealthDetailListResult }, name::Symbol, val)
end


text:
function update_gamma!(state::State, data::Data, prior::Prior)
  update_gammaC!(state, data, prior)
  update_gammaT!(state, data, prior)
end

function update_gammaC!(state::State, data::Data, prior::Prior)
  a, b = params(prior.gamma)
  anew = a + data.ZC + data.ZT * (1 - state.beta)
  bnew = b + data.NC - data.ZC + (data.NT - data.ZT) * (1 - state.beta)
  state.gammaC = rand(Beta(anew, bnew))
end

function update_gammaT!(state::State, data::Data, prior::Prior)
  a, b = params(prior.gamma)
  anew = a + data.ZT * state.beta
  bnew = b + (data.NT - data.ZT) * state.beta
  state.gammaT = rand(Beta(anew, bnew))
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
using Econometrics, LinearAlgebra, Statistics
function main()
# this illustrates the Cholesky decomposition
# create a random p.d. covariance matrix
z = randn(3,3)
V = z'*z # make it pos. def.
println("true covariance matrix")
prettyprint(V)

# take Cholesky decomp.
P = cholesky(V).U
println("P=chol(V): P'P-V")
prettyprint(P'*P - V)

# how to sample from N(0,V)
errs = randn(10000,3)*P
println("sample covariance matrix of 10000 draws from N(0,V)")
prettyprint(cov(errs))

# verify that GLS transformation works
P = inv(cholesky(cov(errs)).U)
errs *= P
println("sample covariance matrix of transformed errors")
prettyprint(cov(errs))
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics