How to use the gpkit.SignomialsEnabled function in gpkit

" link because %s was set to %s but"
                                         " %s was set to %s" % (
                                             subbed_vk, sub,
                                             vk, self.substitutions[vk]))
            # vks is a set, so it's convenient to use the loop variable here
            # since we've already verified above that vks is not null
            descr = dict(vk.descr)  # pylint: disable=undefined-loop-variable
            descr.pop("value", None)
            descr.pop("models", None)
            descr.pop("modelnums", None)
            newvk = VarKey(**descr)
            if sub:
                self.substitutions[newvk] = sub
            self.linked.update(dict(zip(vks, len(vks)*[newvk])))
            self.reverselinks[newvk] = vks
        with SignomialsEnabled():  # since we're just substituting varkeys.
            self.subinplace(self.linked)
        for constraint in self:
            if hasattr(constraint, "reset_varkeys"):
                constraint.reset_varkeys()
        self.reset_varkeys()

Lval = 0.0065   # [K/m]
        th = GRAVITATIONAL_ACCEL/(GAS_CONSTANT*Lval) # [-]

        L    = Variable('L', Lval, 'K/m', 'Temperature lapse rate')
        p_0  = Variable('p_0', 101325, 'Pa', 'Pressure at sea level')
        T_0  = Variable('T_0', 288.15, 'K', 'Temperature at sea level')

        if min_rho:
            objective = rho  # minimize density
        else:
            objective = 1/rho  # maximize density

        # Temperature lapse rate constraint
        if min_rho:
            with SignomialsEnabled():
                constraints = TCS([T_0 <= T + L*h])
        else:
            constraints = TCS([T_0 >= T + L*h])

        constraints += [h <= 11000*units.m,
                        h >= 1E-6*units.m,

                        # Pressure-altitude relation
                        (p/p_0)**(1/th) == T/T_0,

                        # Ideal gas law
                        rho == p/(R*T),
                        ]
        su = Sutherland()
        lc = su.link(constraints)

po      = Variable('p_{o}','N/m','Center section theoretical wing loading')
        fLo     = Variable('f_{L_{o}}',0.5,'-','Center wing lift reduction coefficient')

        # Wing tip lift reduction variables
        dLt = Variable('\\Delta L_{t}','N','Wing tip lift loss')
        fLt = Variable('f_{L_{t}}',0.05,'-','Wing tip lift reduction coefficient')

        #wing moment variables -- need a good way to model this, currently using TAT
        cmw = Variable('c_{m_{w}}', '-', 'Wing Pitching Moment Coefficient')

        amax    = Variable('\\alpha_{max,w}', '-', 'Max angle of attack')
        
        #make constraints
        constraints = []

        with SignomialsEnabled():
            constraints.extend([
                0.5*state['\\rho']*state['V']**2*self.wing['S']*CLw >= Lw + dLo
                                                                     + 2.*dLt,
                dLo == etao*fLo*self.wing['b']/2*po,
                dLt == fLt*po*self.wing['c_{root}']*self.wing['taper']**2,
                # TODO improve approximations croot~co and taper~gammat

                # DATCOM formula (Mach number makes it SP)
                # Swept wing lift curve slope constraint
                SignomialEquality((self.wing['AR']/self.wing['\\eta'])**2*(1
                                   + self.wing['\\tan(\\Lambda)']**2
                                   - state['M']**2) + 8*pi*self.wing['AR']/CLaw,
                                  (2*pi*self.wing['AR']/CLaw)**2),
 
                CLw == CLaw*alpha,
                alpha <= amax,

th = (g*M)/(R*L)  # dimensionless
        self.g, self.L, self.M, self.R, self.th = g, L, M, R, th

        # Free variables
        h   = Variable('h', 'm', 'Altitude')
        rho = Variable('\\rho', 'kg/m^3', 'Density')
        p   = Variable('p', 'Pa', 'Pressure')
        T   = Variable('T', 'K', 'Temperature')

        # Constants
        L   = Variable('L', L, 'K/m', 'Temperature lapse rate')
        MoR = Variable('(M/R)', M/R, 'kg*K/J', 'Air property')
        p_sl = Variable('p_{sl}', 101325, 'Pa', 'Pressure at sea level')
        T_sl = Variable('T_{sl}', 288.15, 'K', 'Temperature at sea level')

        with SignomialsEnabled():
            objective = rho  # minimize density
            constraints = [  # Model only valid up to top of the troposphere
                             h <= 11000*units.m,
                             # Temperature decreases with height at a rate of L
                             T_sl <= T + L*h,
                             # Pressure-altitude relation and ideal gas law
                             (p/p_sl)**(1/th) == T/T_sl,
                             rho == p*MoR/T]
        return objective, constraints

TCS([W_startCruise2[i] == W_endCruise2[i-1]]),
                ])

        for i in range(0, Nclimb2):
            constraints.extend([
                W_avgClimb2[i] == (W_startClimb2[i]*W_endClimb2[i])**.5,
                TCS([W_startClimb2[i] >= W_endClimb2[i] + W_fuelClimb2[i]]),
                ])

        for i in range(0, Ncruise2):
            constraints.extend([
                W_avgCruise2[i] == (W_startCruise2[i]*W_endCruise2[i])**.5,
                TCS([W_startCruise2[i] >= W_endCruise2[i] + W_fuelCruise2[i]]),
                ])

        with gpkit.SignomialsEnabled():
            if signomial == True:
                constraints.extend([
                    #range constraints
                    TCS([sum(RngClimb1) + sum(RngClimb2) + ReqRngCruise >= ReqRng]),
                    ])
                
            if signomial == False:
                 constraints.extend([
                     TCS([ReqRngCruise >= ReqRng]),
                    ])
                 
            constraints.extend([
                #compute hdclimb2
                hftCruise2 <= dhClimb2 + alt10k,
                ])

T_flight >= Range / V,
                    LoD == C_L/C_D]

        # Thrust and drag model
        C_D_fuse = CDA0 / S
        C_D_wpar = k * C_f * S_wetratio
        C_D_ind  = C_L ** 2 / (np.pi * A * e)
        constraints += [W_f >= TSFC * T_flight * D,
                    D >= 0.5 * rho * S * C_D * V ** 2,
                    C_D >= C_D_fuse + C_D_wpar + C_D_ind,
                    V_f_fuse <= 10*units('m')*CDA0,
                    Re <= (rho / mu) * V * (S / A) ** 0.5,
                    C_f >= 0.074 / Re ** 0.2]

        # Fuel volume model
        with SignomialsEnabled():
            constraints +=[V_f == W_f / g / rho_f,
                    V_f_wing**2 <= 0.0009*S**3/A*tau**2, # linear with b and tau, quadratic with chord
                    V_f_avail <= V_f_wing + V_f_fuse, #[SP]
                    V_f_avail >= V_f
                    ]

        # Wing weight model
        constraints += [W_w_surf >= W_W_coeff2 * S,
                    W_w_strc**2. >= W_W_coeff1**2. / tau**2. * (N_ult**2. * A ** 3. * ((W_0+V_f_fuse*g*rho_f) * W * S)),
                    W_w >= W_w_surf + W_w_strc]

        return constraints

tanL   = Variable('\\tan(\\Lambda_{vt})', '-',
                          'Tangent of leading edge sweep (40 deg)')
        taper  = Variable('\\lambda_{vt}', '-', 'Vertical tail taper ratio')
        tau    = Variable('\\tau_{vt}', '-', 'Vertical tail thickness/chord ratio')
        xCGvt  = Variable('x_{CG_{vt}}', 'm', 'x-location of tail CG')
        y_eng  = Variable('y_{eng}', 'm', 'Engine moment arm')
        zmac   = Variable('z_{\\bar{c}_{vt}}', 'm',
                          'Vertical location of mean aerodynamic chord')

        #engine values
        Te     = Variable('T_e', 'N', 'Thrust per engine at takeoff')

        #constraints
        constraints = []

        with SignomialsEnabled():
            #non vectorized constraints
            constraints.extend([
                LvtEO*lvt >= Te*y_eng + Dwm*y_eng,
                # Force moment balance for one engine out condition
                # TASOPT 2.0 p45

                TCS([dxlead + zmac*tanL + 0.25*cma >= lvt]), # [SP]
                # Tail moment arm

                LvtEO == 0.5*rho0*V1**2*Svt*CLvtEO,
                # Vertical tail force (y-direction) for engine out

                TCS([CLvtEO*(1 + clvtEO/(np.pi*e*Avt)) <= clvtEO]),
                #engine out CL computation

                Avt == bvt**2/Svt,

def __init__(self, **kwargs):
        #set up the overeall model for an on design solve
        m6opt = 1
        m8opt = 0
        
        lpc = FanAndLPC()
        combustor = CombustorCooling()
        turbine = Turbine()
        thrust = ExhaustAndThrust()
        size = OnDesignSizing(m6opt, m8opt)

        self.submodels = [lpc, combustor, turbine, thrust, size]
            
        with SignomialsEnabled():

            substitutions = {
##                'P_0': 30,  
                'M_0': 0.8,
                'T_{t_4}': 1250,
                '\pi_f': 1.5,
                '\pi_{lc}': 3,
                '\pi_{hc}': 10,
                'alpha': 8,
                'alphap1': 9,
                'M_{4a}': 1,    #choked turbines
                'M_2': .4,
                'M_{2.5}': .5,
                'hold_{2}': 1+.5*(1.398-1)*.4**2,
                'hold_{2.5}': 1+.5*(1.354-1)*.5**2,
                }

models = [self.compressor , self. combustor, self. turbine, self. thrust, self.fanmap, self.lpcmap, self.hpcmap, self.sizing, self.state, self.engineP]

        #engine weight
        W_engine = Variable('W_{engine}', 'N', 'Weight of a Single Turbofan Engine')

        #engine fan diameter
        df = Variable('d_{f}', 'm', 'Fan Diameter')
        dlpc = Variable('d_{LPC}', 'm', 'LPC Diameter')
        HTRfSub = Variable('HTR_{f_SUB}', '-', '1 - HTRf^2')
        HTRlpcSub = Variable('HTR_{lpc_SUB}', '-', '1 - HTRlpc^2')

        #make the constraints
        constraints = []

        with SignomialsEnabled():

            weight = [
                  W_engine >= ((self.engineP['m_{total}']/(self.engineP['alphap1']*self.engineP['m_{core}'])*self.engineP['m_{core}'])*.220462)*(1684.5+17.7*(self.engineP['\pi_f']*self.engineP['\pi_{lc}']*self.engineP['\pi_{hc}'])/30+1662.2*(self.engineP['\\alpha']/5)**1.2)*self.engineP['dum2'],
                  ]

            diameter = [
                df == (4 * self.sizing['A_2']/(np.pi * HTRfSub))**.5,
                dlpc == (4 * self.sizing['A_{2.5}']/(np.pi * HTRlpcSub))**.5,
                ]

            fmix = [
                #compute f with mixing
                TCS([self.combustor['eta_{B}'] * self.engineP['f'] * self.combustor['h_f'] >= (1-self.combustor['\\alpha_c'])*self.engineP['h_{t_4}']-(1-self.combustor['\\alpha_c'])*self.engineP['h_{t_3}']+self.combustor['Cp_{fuel}']*self.engineP['f']*(self.engineP['T_{t_4}']-self.combustor['T_{t_f}'])]),
                #compute Tt41...mixing causes a temperature drop
                #had to include Tt4 here to prevent it from being pushed down to zero
                SignomialEquality(self.engineP['h_{t_{4.1}}']*self.engineP['fp1'], ((1-self.combustor['\\alpha_c']+self.engineP['f'])*self.engineP['h_{t_4}'] + self.combustor['\\alpha_c']*self.engineP['h_{t_3}'])),

h = Variable("h", "m", "Altitude")

        """
        Dynamic viscosity (mu) as a function of temperature
        References:
        http://www-mdp.eng.cam.ac.uk/web/library/enginfo/aerothermal_dvd_only/aero/
            atmos/atmos.html
        http://www.cfd-online.com/Wiki/Sutherland's_law
        """
        mu  = Variable('\\mu', 'kg/(m*s)', 'Dynamic viscosity')

        T_s = Variable('T_s', 110.4, "K", "Sutherland Temperature")
        C_1 = Variable('C_1', 1.458E-6, "kg/(m*s*K^0.5)",
                       'Sutherland coefficient')

        with SignomialsEnabled():
            constraints = [
                # Pressure-altitude relation
                (p_atm/p_sl)**(1/5.257) == T_atm/T_sl,

                # Ideal gas law
                rho == p_atm/(R_atm/M_atm*T_atm),

                #temperature equation
                SignomialEquality(T_sl, T_atm + L_atm*alt['h']),

                #constraint on mu
                SignomialEquality((T_atm + T_s) * mu, C_1 * T_atm**1.5),
                ]

        #like to use a local subs here in the future
        substitutions = None