text:
# Automatically generated
function eom(u, p, t)
    @unpack ae, af, footang, g, he, hf, ina, inb, inc, ind, ine, inf, ing, inh, ini, k1, k2, k3, k4, k5, k6, k7, k8, ke, kf, l1, l10, l11, l2, l3, l4, l5, l6, l7, l8, l9, ma, mb, mc, md, me, mf, mg, mh, mi, mtpb, mtpk, pop1xi, pop2xi, mt, u8, u9, u10, u11, z, virtual_force = p
    @unpack ea, fa, gs, ha, ia, eap, fap, gsp, hap, iap, eapp, fapp, gspp, happ, iapp = p
    @inbounds q1, q2, q3, q4, q5, q6, q7, u1, u2, u3, u4, u5, u6, u7, θhe, θke, θae, θhf, θkf, θaf = u

    # specified variables
    ea = ea(t)
    fa = fa(t)
    ha = ha(t)
    ia = ia(t)
    eap = eap(t)
    fapp = fapp(t)
    eapp = eapp(t)
    happ = happ(t)
    iapp = iapp(t)
    fap = fap(t)
    hap = hap(t)
    iap = iap(t)
    gs = gs(t)
    gsp = gsp(t)
    gspp = gspp(t)

    # calculated variables
    θh = π + q7
    θk = π - q6
    θa = π + q5
    θm = q4
    ωh = u7
    ωk = -u6
    ωa = u5
    ωm = u4
    pop2y = q2 - l2 * sin(q3 - q4 - q5 - q6 - q7)

    # torques
    hetor = update_torque_generator!(he, t, 2π - θh, -ωh, θhe)
    ketor = update_torque_generator!(ke, t, 2π - θk, -ωk, θke)
    aetor = update_torque_generator!(ae, t, 2π - θa, -ωa, θae)
    hftor = update_torque_generator!(hf, t, θh, ωh, θhf)
    kftor = update_torque_generator!(kf, t, θk, ωk, θkf)
    aftor = update_torque_generator!(af, t, θa, ωa, θaf)

    # passive torques
    hpass, kpass, apass = PEtorque(θh, θk, θa)
    hpass = 0.0

    htor = hetor - hftor - hpass
    ktor = ketor - kftor - kpass
    ator = aetor - aftor + apass
    mtor = mtpk * (π - θm) - mtpb * ωm

    # grf
    if q2 < 0.0
        dp1x = q1 - pop1xi
        ry1 = -k3 * q2 - k4 * abs(q2) * u2
        rx1 = ry1 * (-k1 * dp1x - k2 * u1) # (-K1*(Q1-POP1XI) - K2*U1)*RY1

    else
        rx1 = 0.0
        ry1 = 0.0
    end

    if pop2y < 0.0
        pop2x = q1 - l2 * cos(q3 - q4 - q5 - q6 - q7)
        vop2x = u1 - l2 * sin(q3 - q4 - q5 - q6 - q7) * (u3 - u4 - u5 - u6 - u7)
        vop2y = u2 - l2 * cos(q3 - q4 - q5 - q6 - q7) * (u3 - u4 - u5 - u6 - u7)
        dp2x = pop2x - pop2xi
        ry2 = -k7 * pop2y - k8 * abs(pop2y) * vop2y
        rx2 = ry2 * (-k5 * dp2x - k6 * vop2x)
    else
        rx2 = 0.0
        ry2 = 0.0
    end

    # virtual grf
    vrx1::Float64 = virtual_force.vrx1(t)
    vry1::Float64 = virtual_force.vry1(t)
    vrx2::Float64 = virtual_force.vrx2(t)
    vry2::Float64 = virtual_force.vry2(t)

    # @show t q3 q4 q5 q6 q7

    # z variables
    z[399] = vrx2 + vry2
    z[401] = ator + ktor
    z[402] = -htor - ktor
    z[3] = cos(q4)
    z[5] = cos(q5)
    z[4] = sin(q4)
    z[6] = sin(q5)
    z[27] = z[3] * z[5] - z[4] * z[6]
    z[1] = cos(q3)
    z[9] = cos(q6)
    z[12] = sin(q7)
    z[10] = sin(q6)
    z[11] = cos(q7)
    z[31] = -(z[9]) * z[12] - z[10] * z[11]
    z[2] = sin(q3)
    z[30] = z[9] * z[11] - z[10] * z[12]
    z[34] = z[1] * z[31] + z[2] * z[30]
    z[28] = -(z[3]) * z[6] - z[4] * z[5]
    z[32] = z[9] * z[12] + z[10] * z[11]
    z[36] = z[1] * z[30] + z[2] * z[32]
    z[38] = z[27] * z[34] + z[28] * z[36]
    z[33] = z[1] * z[30] - z[2] * z[31]
    z[35] = z[1] * z[32] - z[2] * z[30]
    z[37] = z[27] * z[33] + z[28] * z[35]
    z[29] = z[3] * z[6] + z[4] * z[5]
    z[40] = z[27] * z[36] + z[29] * z[34]
    z[39] = z[27] * z[35] + z[29] * z[33]
    z[41] = z[7] * z[3] - z[8] * z[4]
    z[42] = z[7] * z[4] + z[8] * z[3]
    z[43] = -(z[7]) * z[4] - z[8] * z[3]
    z[47] = z[38] * z[43] + z[40] * z[41]
    z[48] = z[3] * z[37] + z[4] * z[39]
    z[49] = z[3] * z[38] + z[4] * z[40]
    z[50] = z[3] * z[39] - z[4] * z[37]
    z[51] = z[3] * z[40] - z[4] * z[38]
    z[52] = z[9] * z[33] + z[10] * z[35]
    z[53] = z[9] * z[34] + z[10] * z[36]
    z[54] = z[9] * z[35] - z[10] * z[33]
    z[55] = z[9] * z[36] - z[10] * z[34]
    z[103] = ((-u4 - u5) - u6) - u7
    z[105] = z[1] * z[39] + z[2] * z[40]
    z[106] = z[1] * z[38] - z[2] * z[37]
    z[107] = z[1] * z[40] - z[2] * z[39]
    z[108] = (-u5 - u6) - u7
    z[110] = z[1] * z[50] + z[2] * z[51]
    z[111] = z[1] * z[49] - z[2] * z[48]
    z[112] = z[1] * z[51] - z[2] * z[50]
    z[113] = -u6 - u7
    z[141] = l2 * z[4]
    z[142] = l2 * z[3]
    z[143] = z[142] - l6
    z[144] = l6 - z[142]
    z[145] = z[5] * z[48] + z[6] * z[50]
    z[146] = z[5] * z[49] + z[6] * z[51]
    z[147] = z[5] * z[141] + z[6] * z[142]
    z[148] = z[5] * z[141] + z[6] * z[143]
    z[149] = z[6] * z[144] - z[5] * z[141]
    z[150] = z[5] * z[50] - z[6] * z[48]
    z[151] = z[5] * z[51] - z[6] * z[49]
    z[152] = z[5] * z[142] - z[6] * z[141]
    z[153] = z[5] * z[143] - z[6] * z[141]
    z[154] = z[5] * z[144] + z[6] * z[141]
    z[155] = z[153] - l7
    z[156] = l7 + z[154]
    z[157] = z[153] - l8
    z[158] = l8 + z[154]
    z[159] = z[9] * z[145] + z[10] * z[150]
    z[160] = z[9] * z[146] + z[10] * z[151]
    z[161] = z[9] * z[147] + z[10] * z[152]
    z[162] = z[9] * z[148] + z[10] * z[153]
    z[163] = z[9] * z[148] + z[10] * z[157]
    z[164] = z[9] * z[149] + z[10] * z[158]
    z[165] = z[9] * z[150] - z[10] * z[145]
    z[166] = z[9] * z[151] - z[10] * z[146]
    z[167] = z[9] * z[152] - z[10] * z[147]
    z[168] = z[9] * z[153] - z[10] * z[148]
    z[169] = z[9] * z[157] - z[10] * z[148]
    z[170] = z[9] * z[158] - z[10] * z[149]
    z[171] = z[169] - l9
    z[172] = l9 + z[170]
    z[173] = z[169] - l10
    z[174] = l10 + z[170]
    z[317] = l1 * ((((u3 - u4) - u5) - u6) - u7) * (u3 + z[103])
    z[318] = l2 * ((((u3 - u4) - u5) - u6) - u7) * (u3 + z[103])
    z[319] = l4 * (((u3 - u5) - u6) - u7) * (u3 + z[108])
    z[320] = l3 * (((u3 - u5) - u6) - u7) * (u3 + z[108])
    z[324] = z[3] * z[318]
    z[325] = z[4] * z[318]
    z[326] = z[324] - l6 * (((u3 - u5) - u6) - u7) * (u3 + z[108])
    z[327] = z[5] * z[326] - z[6] * z[325]
    z[328] = -(z[5]) * z[325] - z[6] * z[326]
    z[329] = z[327] - l7 * ((u3 - u6) - u7) * (u3 + z[113])
    z[330] = z[327] - l8 * ((u3 - u6) - u7) * (u3 + z[113])
    z[331] = z[9] * z[330] + z[10] * z[328]
    z[332] = z[9] * z[328] - z[10] * z[330]
    z[333] = z[331] - l9 * (u3 - u7)^2
    z[334] = z[331] - l10 * (u3 - u7)^2
    z[374] = l2 * ((((u3 - u4) - u5) - u6) - u7)
    z[375] = l6 * (((u3 - u5) - u6) - u7)
    z[376] = l8 * ((u3 - u6) - u7)
    z[377] = l10 * (u3 - u7)
    z[380] = ((((u3 - u4) - u5) - u6) - u7) * z[374]
    z[381] = (((u3 - u5) - u6) - u7) * z[375]
    z[382] = ((u3 - u6) - u7) * z[376]
    z[383] = (u3 - u7) * z[377]
    z[400] = mtor - ator
    z[403] = z[9] * z[27] + z[10] * z[28]
    z[404] = z[9] * z[28] - z[10] * z[27]
    z[405] = z[9] * z[29] + z[10] * z[27]
    z[406] = z[9] * z[27] - z[10] * z[29]
    z[408] = z[3] * z[404] + z[4] * z[406]
    z[409] = z[3] * z[405] - z[4] * z[403]
    z[410] = z[3] * z[406] - z[4] * z[404]
    z[13] = cos(ea)
    z[14] = sin(ea)
    z[15] = cos(fa)
    z[16] = sin(fa)
    z[17] = cos(ha)
    z[18] = sin(ha)
    z[19] = cos(ia)
    z[20] = sin(ia)
    z[56] = z[13] * z[1] + z[14] * z[2]
    z[57] = z[13] * z[2] - z[14] * z[1]
    z[58] = z[14] * z[1] - z[13] * z[2]
    z[59] = -(z[15]) * z[56] - z[16] * z[58]
    z[60] = -(z[15]) * z[57] - z[16] * z[56]
    z[61] = z[16] * z[56] - z[15] * z[58]
    z[62] = z[16] * z[57] - z[15] * z[56]
    z[63] = z[18] * z[61] - z[17] * z[59]
    z[64] = z[18] * z[62] - z[17] * z[60]
    z[65] = -(z[17]) * z[61] - z[18] * z[59]
    z[66] = -(z[17]) * z[62] - z[18] * z[60]
    z[70] = z[7] * z[66] - z[8] * z[64]
    z[71] = z[20] * z[65] - z[19] * z[63]
    z[72] = z[20] * z[66] - z[19] * z[64]
    z[73] = -z[19] * z[65] - z[20] * z[63]
    z[74] = -(z[19]) * z[66] - z[20] * z[64]
    z[114] = u8 - eap
    z[116] = fapp - eapp
    z[122] = (fapp - eapp) - happ
    z[128] = ((fapp - eapp) - happ) - iapp
    z[175] = z[13] * z[11] + z[14] * z[12]
    z[176] = z[14] * z[11] - z[13] * z[12]
    z[177] = z[13] * z[12] - z[14] * z[11]
    z[178] = z[159] * z[175] + z[165] * z[177]
    z[179] = z[160] * z[175] + z[166] * z[177]
    z[180] = z[161] * z[175] + z[167] * z[177]
    z[181] = z[162] * z[175] + z[168] * z[177]
    z[182] = z[163] * z[175] + z[169] * z[177]
    z[183] = z[163] * z[175] + z[173] * z[177]
    z[184] = z[164] * z[175] + z[174] * z[177]
    z[185] = z[159] * z[176] + z[165] * z[175]
    z[186] = z[160] * z[176] + z[166] * z[175]
    z[187] = z[161] * z[176] + z[167] * z[175]
    z[188] = z[162] * z[176] + z[168] * z[175]
    z[189] = z[163] * z[176] + z[169] * z[175]
    z[190] = z[163] * z[176] + z[173] * z[175]
    z[191] = z[164] * z[176] + z[174] * z[175]
    z[192] = z[21] * eap
    z[193] = z[21] + z[191]
    z[194] = l10 * eap
    z[195] = l10 + z[191]
    z[197] = -(z[15]) * z[178] - z[16] * z[185]
    z[198] = -(z[15]) * z[179] - z[16] * z[186]
    z[199] = -(z[15]) * z[180] - z[16] * z[187]
    z[200] = -(z[15]) * z[181] - z[16] * z[188]
    z[201] = -(z[15]) * z[182] - z[16] * z[189]
    z[202] = -(z[15]) * z[183] - z[16] * z[190]
    z[203] = -(z[15]) * z[184] - z[16] * z[195]
    z[205] = z[16] * z[178] - z[15] * z[185]
    z[206] = z[16] * z[179] - z[15] * z[186]
    z[207] = z[16] * z[180] - z[15] * z[187]
    z[208] = z[16] * z[181] - z[15] * z[188]
    z[209] = z[16] * z[182] - z[15] * z[189]
    z[210] = z[16] * z[183] - z[15] * z[190]
    z[211] = z[16] * z[184] - z[15] * z[195]
    z[214] = fap - eap
    z[215] = z[22] * z[214]
    z[216] = z[22] + z[211]
    z[219] = l8 * z[214]
    z[220] = l8 + z[211]
    z[223] = (fap - eap) - hap
    z[228] = z[24] * z[223]
    z[229] = z[25] * z[223]
    z[230] = l6 * z[223]
    z[234] = z[18] * z[205] - z[17] * z[197]
    z[235] = z[18] * z[206] - z[17] * z[198]
    z[236] = z[18] * z[207] - z[17] * z[199]
    z[237] = z[18] * z[208] - z[17] * z[200]
    z[238] = z[18] * z[209] - z[17] * z[201]
    z[239] = z[18] * z[210] - z[17] * z[202]
    z[240] = z[18] * z[220] - z[17] * z[203]
    z[243] = -(z[17]) * z[205] - z[18] * z[197]
    z[244] = -(z[17]) * z[206] - z[18] * z[198]
    z[245] = -(z[17]) * z[207] - z[18] * z[199]
    z[246] = -(z[17]) * z[208] - z[18] * z[200]
    z[247] = -(z[17]) * z[209] - z[18] * z[201]
    z[248] = -(z[17]) * z[210] - z[18] * z[202]
    z[249] = -(z[17]) * z[220] - z[18] * z[203]
    z[252] = l6 + z[249]
    z[259] = z[20] * z[243] - z[19] * z[234]
    z[260] = z[20] * z[244] - z[19] * z[235]
    z[261] = z[20] * z[245] - z[19] * z[236]
    z[262] = z[20] * z[246] - z[19] * z[237]
    z[263] = z[20] * z[247] - z[19] * z[238]
    z[264] = z[20] * z[248] - z[19] * z[239]
    z[265] = z[20] * z[252] - z[19] * z[240]
    z[269] = -(z[19]) * z[243] - z[20] * z[234]
    z[270] = -(z[19]) * z[244] - z[20] * z[235]
    z[271] = -(z[19]) * z[245] - z[20] * z[236]
    z[272] = -(z[19]) * z[246] - z[20] * z[237]
    z[273] = -(z[19]) * z[247] - z[20] * z[238]
    z[274] = -(z[19]) * z[248] - z[20] * z[239]
    z[275] = -(z[19]) * z[252] - z[20] * z[240]
    z[278] = ((fap - eap) - hap) - iap
    z[279] = z[26] * z[278]
    z[281] = z[26] + z[275]
    z[287] = l2 + z[275]
    z[335] = z[175] * z[334] + z[177] * z[332]
    z[336] = z[175] * z[332] + z[176] * z[334]
    z[337] = z[21] * eapp
    z[338] = z[335] + ((z[192] - z[21] * u3) - z[21] * u8) * (u3 + z[114])
    z[339] = z[336] - z[337]
    z[340] = l10 * eapp
    z[341] = z[335] + ((z[194] - l10 * u3) - l10 * u8) * (u3 + z[114])
    z[342] = z[336] - z[340]
    z[343] = -(z[15]) * z[341] - z[16] * z[342]
    z[344] = z[16] * z[341] - z[15] * z[342]
    z[345] = z[22] * z[116]
    z[346] = z[214] + u8 + u9
    z[347] = z[343] - (z[215] + z[22] * u3 + z[22] * u8 + z[22] * u9) * (u3 + z[346])
    z[348] = z[345] + z[344]
    z[349] = l8 * z[116]
    z[350] = z[343] - (z[219] + l8 * u3 + l8 * u8 + l8 * u9) * (u3 + z[346])
    z[351] = z[349] + z[344]
    z[352] = z[24] * z[122]
    z[353] = z[25] * z[122]
    z[354] = z[223] + u10 + u8 + u9
    z[355] = (z[228] + z[24] * u10 + z[24] * u3 + z[24] * u8 + z[24] * u9) * (u3 + z[354])
    z[356] = (z[229] + z[25] * u10 + z[25] * u3 + z[25] * u8 + z[25] * u9) * (u3 + z[354])
    z[357] = l6 * z[122]
    z[361] = z[18] * z[351] - z[17] * z[350]
    z[362] = -(z[17]) * z[351] - z[18] * z[350]
    z[363] = z[361] - (z[230] + l6 * u10 + l6 * u3 + l6 * u8 + l6 * u9) * (u3 + z[354])
    z[364] = z[357] + z[362]
    z[365] = z[20] * z[364] - z[19] * z[363]
    z[366] = -(z[19]) * z[364] - z[20] * z[363]
    z[367] = z[26] * z[128]
    z[368] = z[278] + u10 + u11 + u8 + u9
    z[369] = z[365] - (z[279] + z[26] * u10 + z[26] * u11 + z[26] * u3 + z[26] * u8 + z[26] * u9) * (u3 + z[368])
    z[370] = z[367] + z[366]
    z[378] = gs * u3
    z[379] = gsp * u3
    z[384] = gspp - u3 * z[378]
    z[385] = gsp * u3 + z[379]
    z[411] = -(z[7]) * z[17] - z[8] * z[18]
    z[412] = z[8] * z[17] - z[7] * z[18]
    z[413] = z[7] * z[18] - z[8] * z[17]
    z[739] = rx1 + rx2 * z[37]^2 + rx2 * z[39]^2 + ry2 * z[37] * z[38] + ry2 * z[39] * z[40] + vrx1 * z[71] * z[259] + vrx1 * z[73] * z[269] + vry1 * z[72] * z[259] + vry1 * z[74] * z[269] + z[390] * (z[37] * z[38] + z[39] * z[40]) + z[391] * (z[37] * z[38] + z[39] * z[40]) + z[392] * (z[34] * z[145] + z[36] * z[150]) + z[393] * (z[53] * z[159] + z[55] * z[165]) + z[394] * (z[56] * z[185] + z[57] * z[178]) + z[395] * (z[60] * z[197] + z[62] * z[205]) + z[397] * (z[60] * z[197] + z[62] * z[205]) + z[398] * (z[72] * z[259] + z[74] * z[269]) + z[399] * (z[64] * z[234] + z[66] * z[243])
    z[740] = z[396] + ry1 + rx2 * z[37] * z[38] + rx2 * z[39] * z[40] + ry2 * z[38]^2 + ry2 * z[40]^2 + vrx1 * z[71] * z[260] + vrx1 * z[73] * z[270] + vry1 * z[72] * z[260] + vry1 * z[74] * z[270] + z[390] * (z[38]^2 + z[40]^2) + z[391] * (z[38]^2 + z[40]^2) + z[392] * (z[34] * z[146] + z[36] * z[151]) + z[393] * (z[53] * z[160] + z[55] * z[166]) + z[394] * (z[56] * z[186] + z[57] * z[179]) + z[395] * (z[60] * z[198] + z[62] * z[206]) + z[397] * (z[60] * z[198] + z[62] * z[206]) + z[398] * (z[72] * z[260] + z[74] * z[270]) + z[399] * (z[64] * z[235] + z[66] * z[244])
    z[742] = htor + z[400] + z[401] + z[402] + vrx1 * z[71] * z[265] + vrx1 * z[73] * z[287] + vry1 * z[72] * z[265] + vry1 * z[74] * z[287] + z[392] * (z[34] * z[149] + z[36] * z[156]) + z[393] * (z[53] * z[164] + z[55] * z[172]) + z[394] * (z[56] * z[193] + z[57] * z[184]) + z[395] * (z[60] * z[203] + z[62] * z[216]) + z[398] * (z[72] * z[265] + z[74] * z[281]) + z[399] * (z[64] * z[240] + z[66] * z[252]) + z[397] * (z[24] * z[66] + z[25] * z[70] + z[60] * z[203] + z[62] * z[220]) - mtor - z[741] * z[40] - l2 * (rx2 * z[39] + ry2 * z[40]) - 0.5 * z[391] * (2 * l2 * z[40] - l3 * z[47] - l4 * z[51]) - z[396] * (l2 * z[40] - l10 * z[55] - l6 * z[51] - l8 * z[36] - gs * z[1])
    z[745] = mtor + z[741] * z[40] + z[743] * z[40] + z[744] * z[40] + vrx1 * z[71] * z[261] + vrx1 * z[73] * z[271] + vry1 * z[72] * z[261] + vry1 * z[74] * z[271] + l2 * (rx2 * z[39] + ry2 * z[40]) + z[392] * (z[34] * z[147] + z[36] * z[152]) + z[393] * (z[53] * z[161] + z[55] * z[167]) + z[394] * (z[56] * z[187] + z[57] * z[180]) + z[395] * (z[60] * z[199] + z[62] * z[207]) + z[397] * (z[60] * z[199] + z[62] * z[207]) + z[398] * (z[72] * z[261] + z[74] * z[271]) + z[399] * (z[64] * z[236] + z[66] * z[245])
    z[746] = mtor + z[741] * z[40] + vrx1 * z[71] * z[262] + vrx1 * z[73] * z[272] + vry1 * z[72] * z[262] + vry1 * z[74] * z[272] + l2 * (rx2 * z[39] + ry2 * z[40]) + z[392] * (z[34] * z[148] + z[36] * z[153]) + z[393] * (z[53] * z[162] + z[55] * z[168]) + z[394] * (z[56] * z[188] + z[57] * z[181]) + z[395] * (z[60] * z[200] + z[62] * z[208]) + z[397] * (z[60] * z[200] + z[62] * z[208]) + z[398] * (z[72] * z[262] + z[74] * z[272]) + z[399] * (z[64] * z[237] + z[66] * z[246]) + z[396] * (l2 * z[40] - l6 * z[51]) + 0.5 * z[391] * (2 * l2 * z[40] - l3 * z[47] - l4 * z[51]) - z[400]
    z[747] = mtor + z[741] * z[40] + vrx1 * z[71] * z[263] + vrx1 * z[73] * z[273] + vry1 * z[72] * z[263] + vry1 * z[74] * z[273] + l2 * (rx2 * z[39] + ry2 * z[40]) + z[392] * (z[34] * z[148] + z[36] * z[155]) + z[393] * (z[53] * z[163] + z[55] * z[169]) + z[394] * (z[56] * z[189] + z[57] * z[182]) + z[395] * (z[60] * z[201] + z[62] * z[209]) + z[397] * (z[60] * z[201] + z[62] * z[209]) + z[398] * (z[72] * z[263] + z[74] * z[273]) + z[399] * (z[64] * z[238] + z[66] * z[247]) + z[396] * (l2 * z[40] - l6 * z[51] - l8 * z[36]) + 0.5 * z[391] * (2 * l2 * z[40] - l3 * z[47] - l4 * z[51]) - z[400] - z[401]
    z[748] = mtor + z[741] * z[40] + vrx1 * z[71] * z[264] + vrx1 * z[73] * z[274] + vry1 * z[72] * z[264] + vry1 * z[74] * z[274] + l2 * (rx2 * z[39] + ry2 * z[40]) + z[392] * (z[34] * z[148] + z[36] * z[155]) + z[393] * (z[53] * z[163] + z[55] * z[171]) + z[394] * (z[56] * z[190] + z[57] * z[183]) + z[395] * (z[60] * z[202] + z[62] * z[210]) + z[397] * (z[60] * z[202] + z[62] * z[210]) + z[398] * (z[72] * z[264] + z[74] * z[274]) + z[399] * (z[64] * z[239] + z[66] * z[248]) + 0.5 * z[391] * (2 * l2 * z[40] - l3 * z[47] - l4 * z[51]) + z[396] * (l2 * z[40] - l10 * z[55] - l6 * z[51] - l8 * z[36]) - z[400] - z[401] - z[402]
    z[765] = ine * eapp
    z[767] = inf * z[116]
    z[769] = inh * z[122]
    z[771] = ini * z[128]
    z[772] = mg + ma * (z[37]^2 + z[39]^2) + mb * (z[37]^2 + z[39]^2) + mc * (z[145]^2 + z[150]^2) + md * (z[159]^2 + z[165]^2) + me * (z[178]^2 + z[185]^2) + mf * (z[197]^2 + z[205]^2) + mh * (z[197]^2 + z[205]^2) + mi * (z[259]^2 + z[269]^2)
    z[773] = ma * (z[37] * z[38] + z[39] * z[40]) + mb * (z[37] * z[38] + z[39] * z[40]) + mc * (z[145] * z[146] + z[150] * z[151]) + md * (z[159] * z[160] + z[165] * z[166]) + me * (z[178] * z[179] + z[185] * z[186]) + mf * (z[197] * z[198] + z[205] * z[206]) + mh * (z[197] * z[198] + z[205] * z[206]) + mi * (z[259] * z[260] + z[269] * z[270])
    z[774] = ((mc * (z[145] * z[149] + z[150] * z[156]) + md * (z[159] * z[164] + z[165] * z[172]) + me * (z[178] * z[184] + z[185] * z[193]) + mf * (z[197] * z[203] + z[205] * z[216]) + mi * (z[259] * z[265] + z[269] * z[281]) + mg * (((l10 * z[54] + l6 * z[50] + l8 * z[35]) - l2 * z[39]) - gs * z[2]) + 0.5 * mb * (((l3 * z[37] * z[43] + l3 * z[39] * z[41] + l4 * z[3] * z[39]) - 2 * l2 * z[39]) - l4 * z[4] * z[37])) - z[618] * z[39]) - mh * (((((z[24] * z[17] * z[205] + z[24] * z[18] * z[197]) - z[197] * z[203]) - z[205] * z[220]) - z[25] * z[411] * z[205]) - z[25] * z[412] * z[197])
    z[775] = z[618] * z[39] + z[623] * z[39] + z[687] * z[39] + mc * (z[145] * z[147] + z[150] * z[152]) + md * (z[159] * z[161] + z[165] * z[167]) + me * (z[178] * z[180] + z[185] * z[187]) + mf * (z[197] * z[199] + z[205] * z[207]) + mh * (z[197] * z[199] + z[205] * z[207]) + mi * (z[259] * z[261] + z[269] * z[271])
    z[776] = (z[618] * z[39] + mc * (z[145] * z[148] + z[150] * z[153]) + md * (z[159] * z[162] + z[165] * z[168]) + me * (z[178] * z[181] + z[185] * z[188]) + mf * (z[197] * z[200] + z[205] * z[208]) + mh * (z[197] * z[200] + z[205] * z[208]) + mi * (z[259] * z[262] + z[269] * z[272]) + mg * (l2 * z[39] - l6 * z[50])) - 0.5 * mb * (((l3 * z[37] * z[43] + l3 * z[39] * z[41] + l4 * z[3] * z[39]) - 2 * l2 * z[39]) - l4 * z[4] * z[37])
    z[777] = (z[618] * z[39] + mc * (z[145] * z[148] + z[150] * z[155]) + md * (z[159] * z[163] + z[165] * z[169]) + me * (z[178] * z[182] + z[185] * z[189]) + mf * (z[197] * z[201] + z[205] * z[209]) + mh * (z[197] * z[201] + z[205] * z[209]) + mi * (z[259] * z[263] + z[269] * z[273]) + mg * ((l2 * z[39] - l6 * z[50]) - l8 * z[35])) - 0.5 * mb * (((l3 * z[37] * z[43] + l3 * z[39] * z[41] + l4 * z[3] * z[39]) - 2 * l2 * z[39]) - l4 * z[4] * z[37])
    z[778] = (z[618] * z[39] + mc * (z[145] * z[148] + z[150] * z[155]) + md * (z[159] * z[163] + z[165] * z[171]) + me * (z[178] * z[183] + z[185] * z[190]) + mf * (z[197] * z[202] + z[205] * z[210]) + mh * (z[197] * z[202] + z[205] * z[210]) + mi * (z[259] * z[264] + z[269] * z[274]) + mg * (((l2 * z[39] - l10 * z[54]) - l6 * z[50]) - l8 * z[35])) - 0.5 * mb * (((l3 * z[37] * z[43] + l3 * z[39] * z[41] + l4 * z[3] * z[39]) - 2 * l2 * z[39]) - l4 * z[4] * z[37])
    z[779] = ma * z[37] * z[317] + mc * (z[145] * z[329] + z[150] * z[328]) + md * (z[159] * z[333] + z[165] * z[332]) + me * (z[178] * z[338] + z[185] * z[339]) + mf * (z[197] * z[347] + z[205] * z[348]) + mi * (z[259] * z[369] + z[269] * z[370]) + mg * (((((z[1] * z[384] + z[37] * z[380]) - z[2] * z[385]) - z[33] * z[382]) - z[48] * z[381]) - z[52] * z[383]) + 0.5 * mb * ((((2 * z[37] * z[318] - z[3] * z[37] * z[319]) - z[4] * z[39] * z[319]) - z[37] * z[41] * z[320]) - z[39] * z[42] * z[320]) + mh * ((((((z[411] * z[353] * z[205] + z[412] * z[353] * z[197] + z[197] * z[350] + z[205] * z[351] + z[17] * z[197] * z[355]) - z[17] * z[352] * z[205]) - z[18] * z[352] * z[197]) - z[18] * z[205] * z[355]) - z[411] * z[197] * z[356]) - z[413] * z[205] * z[356])
    z[780] = mg + ma * (z[38]^2 + z[40]^2) + mb * (z[38]^2 + z[40]^2) + mc * (z[146]^2 + z[151]^2) + md * (z[160]^2 + z[166]^2) + me * (z[179]^2 + z[186]^2) + mf * (z[198]^2 + z[206]^2) + mh * (z[198]^2 + z[206]^2) + mi * (z[260]^2 + z[270]^2)
    z[781] = (((mc * (z[146] * z[149] + z[151] * z[156]) + md * (z[160] * z[164] + z[166] * z[172]) + me * (z[179] * z[184] + z[186] * z[193]) + mf * (z[198] * z[203] + z[206] * z[216]) + mi * (z[260] * z[265] + z[270] * z[281]) + 0.5 * mb * (((l3 * z[38] * z[43] + l3 * z[40] * z[41] + l4 * z[3] * z[40]) - 2 * l2 * z[40]) - l4 * z[4] * z[38])) - z[618] * z[40]) - mg * ((((l2 * z[40] - l10 * z[55]) - l6 * z[51]) - l8 * z[36]) - gs * z[1])) - mh * (((((z[24] * z[17] * z[206] + z[24] * z[18] * z[198]) - z[198] * z[203]) - z[206] * z[220]) - z[25] * z[411] * z[206]) - z[25] * z[412] * z[198])
    z[782] = z[618] * z[40] + z[623] * z[40] + z[687] * z[40] + mc * (z[146] * z[147] + z[151] * z[152]) + md * (z[160] * z[161] + z[166] * z[167]) + me * (z[179] * z[180] + z[186] * z[187]) + mf * (z[198] * z[199] + z[206] * z[207]) + mh * (z[198] * z[199] + z[206] * z[207]) + mi * (z[260] * z[261] + z[270] * z[271])
    z[783] = (z[618] * z[40] + mc * (z[146] * z[148] + z[151] * z[153]) + md * (z[160] * z[162] + z[166] * z[168]) + me * (z[179] * z[181] + z[186] * z[188]) + mf * (z[198] * z[200] + z[206] * z[208]) + mh * (z[198] * z[200] + z[206] * z[208]) + mi * (z[260] * z[262] + z[270] * z[272]) + mg * (l2 * z[40] - l6 * z[51])) - 0.5 * mb * (((l3 * z[38] * z[43] + l3 * z[40] * z[41] + l4 * z[3] * z[40]) - 2 * l2 * z[40]) - l4 * z[4] * z[38])
    z[784] = (z[618] * z[40] + mc * (z[146] * z[148] + z[151] * z[155]) + md * (z[160] * z[163] + z[166] * z[169]) + me * (z[179] * z[182] + z[186] * z[189]) + mf * (z[198] * z[201] + z[206] * z[209]) + mh * (z[198] * z[201] + z[206] * z[209]) + mi * (z[260] * z[263] + z[270] * z[273]) + mg * ((l2 * z[40] - l6 * z[51]) - l8 * z[36])) - 0.5 * mb * (((l3 * z[38] * z[43] + l3 * z[40] * z[41] + l4 * z[3] * z[40]) - 2 * l2 * z[40]) - l4 * z[4] * z[38])
    z[785] = (z[618] * z[40] + mc * (z[146] * z[148] + z[151] * z[155]) + md * (z[160] * z[163] + z[166] * z[171]) + me * (z[179] * z[183] + z[186] * z[190]) + mf * (z[198] * z[202] + z[206] * z[210]) + mh * (z[198] * z[202] + z[206] * z[210]) + mi * (z[260] * z[264] + z[270] * z[274]) + mg * (((l2 * z[40] - l10 * z[55]) - l6 * z[51]) - l8 * z[36])) - 0.5 * mb * (((l3 * z[38] * z[43] + l3 * z[40] * z[41] + l4 * z[3] * z[40]) - 2 * l2 * z[40]) - l4 * z[4] * z[38])
    z[786] = ma * z[38] * z[317] + mc * (z[146] * z[329] + z[151] * z[328]) + md * (z[160] * z[333] + z[166] * z[332]) + me * (z[179] * z[338] + z[186] * z[339]) + mf * (z[198] * z[347] + z[206] * z[348]) + mi * (z[260] * z[369] + z[270] * z[370]) + mg * ((((z[1] * z[385] + z[2] * z[384] + z[38] * z[380]) - z[34] * z[382]) - z[49] * z[381]) - z[53] * z[383]) + 0.5 * mb * ((((2 * z[38] * z[318] - z[3] * z[38] * z[319]) - z[4] * z[40] * z[319]) - z[38] * z[41] * z[320]) - z[40] * z[42] * z[320]) + mh * ((((((z[411] * z[353] * z[206] + z[412] * z[353] * z[198] + z[198] * z[350] + z[206] * z[351] + z[17] * z[198] * z[355]) - z[17] * z[352] * z[206]) - z[18] * z[352] * z[198]) - z[18] * z[206] * z[355]) - z[411] * z[198] * z[356]) - z[413] * z[206] * z[356])
    z[804] = ((z[787] + mc * (z[149]^2 + z[156]^2) + md * (z[164]^2 + z[172]^2) + me * (z[184]^2 + z[193]^2) + mf * (z[203]^2 + z[216]^2) + mi * (z[265]^2 + z[281]^2) + 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mh * ((((((((2 * z[24] * z[17] * z[220] + 2 * z[24] * z[18] * z[203]) - 2 * z[793]) - z[791]) - z[792]) - z[203]^2) - z[220]^2) - 2 * z[25] * z[411] * z[220]) - 2 * z[25] * z[412] * z[203])) - mg * ((((((((((((2 * z[794] * z[406] + 2 * z[795] * z[3] + 2 * z[796] * z[27] + 2 * l2 * gs * z[107]) - z[797]) - z[798]) - z[799]) - z[800]) - gs^2) - 2 * z[801] * z[410]) - 2 * z[802] * z[9]) - 2 * z[803] * z[5]) - 2 * l10 * gs * z[11]) - 2 * l6 * gs * z[112]) - 2 * l8 * gs * z[30])
    z[806] = (((((mc * (z[147] * z[149] + z[152] * z[156]) + md * (z[161] * z[164] + z[167] * z[172]) + me * (z[180] * z[184] + z[187] * z[193]) + mf * (z[199] * z[203] + z[207] * z[216]) + mi * (z[261] * z[265] + z[271] * z[281])) - ina) - z[805]) - 0.5 * z[623] * ((2l2 - l3 * z[41]) - l4 * z[3])) - z[687] * ((((l2 - l10 * z[406]) - l6 * z[3]) - l8 * z[27]) - gs * z[107])) - mh * (((((z[24] * z[17] * z[207] + z[24] * z[18] * z[199]) - z[199] * z[203]) - z[207] * z[220]) - z[25] * z[411] * z[207]) - z[25] * z[412] * z[199])
    z[807] = (((((mc * (z[148] * z[149] + z[153] * z[156]) + md * (z[162] * z[164] + z[168] * z[172]) + me * (z[181] * z[184] + z[188] * z[193]) + mf * (z[200] * z[203] + z[208] * z[216]) + mi * (z[262] * z[265] + z[272] * z[281]) + mg * ((((((z[794] * z[406] + z[796] * z[27] + 2 * z[795] * z[3] + l2 * gs * z[107]) - z[798]) - z[799]) - z[801] * z[410]) - z[803] * z[5]) - l6 * gs * z[112])) - ina) - inb) - z[805]) - 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mh * (((((z[24] * z[17] * z[208] + z[24] * z[18] * z[200]) - z[200] * z[203]) - z[208] * z[220]) - z[25] * z[411] * z[208]) - z[25] * z[412] * z[200])
    z[808] = ((((((mc * (z[148] * z[149] + z[155] * z[156]) + md * (z[163] * z[164] + z[169] * z[172]) + me * (z[182] * z[184] + z[189] * z[193]) + mf * (z[201] * z[203] + z[209] * z[216]) + mi * (z[263] * z[265] + z[273] * z[281]) + mg * (((((((((z[794] * z[406] + 2 * z[795] * z[3] + 2 * z[796] * z[27] + l2 * gs * z[107]) - z[798]) - z[799]) - z[800]) - 2 * z[803] * z[5]) - z[801] * z[410]) - z[802] * z[9]) - l6 * gs * z[112]) - l8 * gs * z[30])) - ina) - inb) - inc) - z[805]) - 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mh * (((((z[24] * z[17] * z[209] + z[24] * z[18] * z[201]) - z[201] * z[203]) - z[209] * z[220]) - z[25] * z[411] * z[209]) - z[25] * z[412] * z[201])
    z[809] = (((((((mc * (z[148] * z[149] + z[155] * z[156]) + md * (z[163] * z[164] + z[171] * z[172]) + me * (z[183] * z[184] + z[190] * z[193]) + mf * (z[202] * z[203] + z[210] * z[216]) + mi * (z[264] * z[265] + z[274] * z[281]) + mg * (((((((((((2 * z[794] * z[406] + 2 * z[795] * z[3] + 2 * z[796] * z[27] + l2 * gs * z[107]) - z[797]) - z[798]) - z[799]) - z[800]) - 2 * z[801] * z[410]) - 2 * z[802] * z[9]) - 2 * z[803] * z[5]) - l10 * gs * z[11]) - l6 * gs * z[112]) - l8 * gs * z[30])) - ina) - inb) - inc) - ind) - z[805]) - 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mh * (((((z[24] * z[17] * z[210] + z[24] * z[18] * z[202]) - z[202] * z[203]) - z[210] * z[220]) - z[25] * z[411] * z[210]) - z[25] * z[412] * z[202])
    z[816] = (z[767] + z[769] + z[771] + mc * (z[149] * z[329] + z[156] * z[328]) + md * (z[164] * z[333] + z[172] * z[332]) + me * (z[184] * z[338] + z[193] * z[339]) + mf * (z[203] * z[347] + z[216] * z[348]) + mi * (z[265] * z[369] + z[281] * z[370]) + 0.25 * mb * (((z[810] * z[319] + 2 * l2 * z[4] * z[319] + 2 * l2 * z[42] * z[320] + 2 * l3 * z[43] * z[318]) - z[811] * z[320]) - 2 * l4 * z[4] * z[318]) + mh * (((((((((z[24] * z[352] + z[25] * z[353] + z[812] * z[353] + z[813] * z[352] + z[411] * z[353] * z[220] + z[412] * z[353] * z[203] + z[814] * z[355] + z[203] * z[350] + z[220] * z[351] + z[25] * z[411] * z[351] + z[25] * z[412] * z[350] + z[17] * z[203] * z[355]) - z[17] * z[352] * z[220]) - z[18] * z[352] * z[203]) - z[815] * z[356]) - z[24] * z[17] * z[351]) - z[24] * z[18] * z[350]) - z[18] * z[220] * z[355]) - z[411] * z[203] * z[356]) - z[413] * z[220] * z[356]) + mg * ((((((((((gs * z[385] + l10 * z[10] * z[382] + l10 * z[11] * z[385] + l10 * z[12] * z[384] + l10 * z[404] * z[380] + l2 * z[4] * z[381] + l2 * z[29] * z[382] + l2 * z[405] * z[383] + l6 * z[110] * z[384] + l6 * z[112] * z[385] + l8 * z[6] * z[381] + l8 * z[28] * z[380] + l8 * z[30] * z[385] + l8 * z[32] * z[384] + gs * z[12] * z[383] + gs * z[106] * z[380]) - l10 * z[408] * z[381]) - l2 * z[105] * z[384]) - l2 * z[107] * z[385]) - l6 * z[4] * z[380]) - l6 * z[6] * z[382]) - l6 * z[409] * z[383]) - l8 * z[10] * z[383]) - gs * z[31] * z[382]) - gs * z[111] * z[381])) - z[765]
    z[818] = z[817] + mc * (z[147]^2 + z[152]^2) + md * (z[161]^2 + z[167]^2) + me * (z[180]^2 + z[187]^2) + mf * (z[199]^2 + z[207]^2) + mh * (z[199]^2 + z[207]^2) + mi * (z[261]^2 + z[271]^2)
    z[820] = z[819] + mc * (z[147] * z[148] + z[152] * z[153]) + md * (z[161] * z[162] + z[167] * z[168]) + me * (z[180] * z[181] + z[187] * z[188]) + mf * (z[199] * z[200] + z[207] * z[208]) + mh * (z[199] * z[200] + z[207] * z[208]) + mi * (z[261] * z[262] + z[271] * z[272]) + z[687] * (l2 - l6 * z[3]) + 0.5 * z[623] * ((2l2 - l3 * z[41]) - l4 * z[3])
    z[821] = z[819] + mc * (z[147] * z[148] + z[152] * z[155]) + md * (z[161] * z[163] + z[167] * z[169]) + me * (z[180] * z[182] + z[187] * z[189]) + mf * (z[199] * z[201] + z[207] * z[209]) + mh * (z[199] * z[201] + z[207] * z[209]) + mi * (z[261] * z[263] + z[271] * z[273]) + z[687] * ((l2 - l6 * z[3]) - l8 * z[27]) + 0.5 * z[623] * ((2l2 - l3 * z[41]) - l4 * z[3])
    z[822] = z[819] + mc * (z[147] * z[148] + z[152] * z[155]) + md * (z[161] * z[163] + z[167] * z[171]) + me * (z[180] * z[183] + z[187] * z[190]) + mf * (z[199] * z[202] + z[207] * z[210]) + mh * (z[199] * z[202] + z[207] * z[210]) + mi * (z[261] * z[264] + z[271] * z[274]) + 0.5 * z[623] * ((2l2 - l3 * z[41]) - l4 * z[3]) + z[687] * (((l2 - l10 * z[406]) - l6 * z[3]) - l8 * z[27])
    z[823] = (mc * (z[147] * z[329] + z[152] * z[328]) + md * (z[161] * z[333] + z[167] * z[332]) + me * (z[180] * z[338] + z[187] * z[339]) + mf * (z[199] * z[347] + z[207] * z[348]) + mi * (z[261] * z[369] + z[271] * z[370]) + z[687] * ((((z[105] * z[384] + z[107] * z[385]) - z[4] * z[381]) - z[29] * z[382]) - z[405] * z[383]) + mh * ((((((z[411] * z[353] * z[207] + z[412] * z[353] * z[199] + z[199] * z[350] + z[207] * z[351] + z[17] * z[199] * z[355]) - z[17] * z[352] * z[207]) - z[18] * z[352] * z[199]) - z[18] * z[207] * z[355]) - z[411] * z[199] * z[356]) - z[413] * z[207] * z[356])) - 0.5 * z[623] * (z[4] * z[319] + z[42] * z[320])
    z[826] = z[824] + mc * (z[148]^2 + z[153]^2) + md * (z[162]^2 + z[168]^2) + me * (z[181]^2 + z[188]^2) + mf * (z[200]^2 + z[208]^2) + mg * (z[825] - 2 * z[795] * z[3]) + mh * (z[200]^2 + z[208]^2) + mi * (z[262]^2 + z[272]^2) + 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])
    z[827] = (z[824] + mc * (z[148]^2 + z[153] * z[155]) + md * (z[162] * z[163] + z[168] * z[169]) + me * (z[181] * z[182] + z[188] * z[189]) + mf * (z[200] * z[201] + z[208] * z[209]) + mh * (z[200] * z[201] + z[208] * z[209]) + mi * (z[262] * z[263] + z[272] * z[273]) + 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mg * ((((z[796] * z[27] + 2 * z[795] * z[3]) - z[798]) - z[799]) - z[803] * z[5])
    z[828] = (z[824] + mc * (z[148]^2 + z[153] * z[155]) + md * (z[162] * z[163] + z[168] * z[171]) + me * (z[181] * z[183] + z[188] * z[190]) + mf * (z[200] * z[202] + z[208] * z[210]) + mh * (z[200] * z[202] + z[208] * z[210]) + mi * (z[262] * z[264] + z[272] * z[274]) + 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mg * (((((z[794] * z[406] + z[796] * z[27] + 2 * z[795] * z[3]) - z[798]) - z[799]) - z[801] * z[410]) - z[803] * z[5])
    z[829] = ((mc * (z[148] * z[329] + z[153] * z[328]) + md * (z[162] * z[333] + z[168] * z[332]) + me * (z[181] * z[338] + z[188] * z[339]) + mf * (z[200] * z[347] + z[208] * z[348]) + mi * (z[262] * z[369] + z[272] * z[370]) + mh * ((((((z[411] * z[353] * z[208] + z[412] * z[353] * z[200] + z[200] * z[350] + z[208] * z[351] + z[17] * z[200] * z[355]) - z[17] * z[352] * z[208]) - z[18] * z[352] * z[200]) - z[18] * z[208] * z[355]) - z[411] * z[200] * z[356]) - z[413] * z[208] * z[356])) - 0.25 * mb * (((z[810] * z[319] + 2 * l2 * z[4] * z[319] + 2 * l2 * z[42] * z[320] + 2 * l3 * z[43] * z[318]) - z[811] * z[320]) - 2 * l4 * z[4] * z[318])) - mg * ((((((l2 * z[4] * z[381] + l2 * z[29] * z[382] + l2 * z[405] * z[383] + l6 * z[110] * z[384] + l6 * z[112] * z[385]) - l2 * z[105] * z[384]) - l2 * z[107] * z[385]) - l6 * z[4] * z[380]) - l6 * z[6] * z[382]) - l6 * z[409] * z[383])
    z[831] = (z[830] + mc * (z[148]^2 + z[155]^2) + md * (z[163]^2 + z[169]^2) + me * (z[182]^2 + z[189]^2) + mf * (z[201]^2 + z[209]^2) + mh * (z[201]^2 + z[209]^2) + mi * (z[263]^2 + z[273]^2) + 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mg * (((((2 * z[795] * z[3] + 2 * z[796] * z[27]) - z[798]) - z[799]) - z[800]) - 2 * z[803] * z[5])
    z[832] = (z[830] + mc * (z[148]^2 + z[155]^2) + md * (z[163]^2 + z[169] * z[171]) + me * (z[182] * z[183] + z[189] * z[190]) + mf * (z[201] * z[202] + z[209] * z[210]) + mh * (z[201] * z[202] + z[209] * z[210]) + mi * (z[263] * z[264] + z[273] * z[274]) + 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mg * (((((((z[794] * z[406] + 2 * z[795] * z[3] + 2 * z[796] * z[27]) - z[798]) - z[799]) - z[800]) - 2 * z[803] * z[5]) - z[801] * z[410]) - z[802] * z[9])
    z[833] = (mc * (z[148] * z[329] + z[155] * z[328]) + md * (z[163] * z[333] + z[169] * z[332]) + me * (z[182] * z[338] + z[189] * z[339]) + mf * (z[201] * z[347] + z[209] * z[348]) + mi * (z[263] * z[369] + z[273] * z[370]) + mh * ((((((z[411] * z[353] * z[209] + z[412] * z[353] * z[201] + z[201] * z[350] + z[209] * z[351] + z[17] * z[201] * z[355]) - z[17] * z[352] * z[209]) - z[18] * z[352] * z[201]) - z[18] * z[209] * z[355]) - z[411] * z[201] * z[356]) - z[413] * z[209] * z[356]) + mg * ((((((((((l2 * z[105] * z[384] + l2 * z[107] * z[385] + l6 * z[4] * z[380] + l6 * z[6] * z[382] + l6 * z[409] * z[383] + l8 * z[10] * z[383]) - l2 * z[4] * z[381]) - l2 * z[29] * z[382]) - l2 * z[405] * z[383]) - l6 * z[110] * z[384]) - l6 * z[112] * z[385]) - l8 * z[6] * z[381]) - l8 * z[28] * z[380]) - l8 * z[30] * z[385]) - l8 * z[32] * z[384])) - 0.25 * mb * (((z[810] * z[319] + 2 * l2 * z[4] * z[319] + 2 * l2 * z[42] * z[320] + 2 * l3 * z[43] * z[318]) - z[811] * z[320]) - 2 * l4 * z[4] * z[318])
    z[835] = (z[834] + mc * (z[148]^2 + z[155]^2) + md * (z[163]^2 + z[171]^2) + me * (z[183]^2 + z[190]^2) + mf * (z[202]^2 + z[210]^2) + mh * (z[202]^2 + z[210]^2) + mi * (z[264]^2 + z[274]^2) + 0.25 * mb * ((z[788] - 4 * z[789] * z[41]) - 4 * z[790] * z[3])) - mg * ((((((((2 * z[794] * z[406] + 2 * z[795] * z[3] + 2 * z[796] * z[27]) - z[797]) - z[798]) - z[799]) - z[800]) - 2 * z[801] * z[410]) - 2 * z[802] * z[9]) - 2 * z[803] * z[5])
    z[836] = (mc * (z[148] * z[329] + z[155] * z[328]) + md * (z[163] * z[333] + z[171] * z[332]) + me * (z[183] * z[338] + z[190] * z[339]) + mf * (z[202] * z[347] + z[210] * z[348]) + mi * (z[264] * z[369] + z[274] * z[370]) + mh * ((((((z[411] * z[353] * z[210] + z[412] * z[353] * z[202] + z[202] * z[350] + z[210] * z[351] + z[17] * z[202] * z[355]) - z[17] * z[352] * z[210]) - z[18] * z[352] * z[202]) - z[18] * z[210] * z[355]) - z[411] * z[202] * z[356]) - z[413] * z[210] * z[356]) + mg * ((((((((((((((l10 * z[408] * z[381] + l2 * z[105] * z[384] + l2 * z[107] * z[385] + l6 * z[4] * z[380] + l6 * z[6] * z[382] + l6 * z[409] * z[383] + l8 * z[10] * z[383]) - l10 * z[10] * z[382]) - l10 * z[11] * z[385]) - l10 * z[12] * z[384]) - l10 * z[404] * z[380]) - l2 * z[4] * z[381]) - l2 * z[29] * z[382]) - l2 * z[405] * z[383]) - l6 * z[110] * z[384]) - l6 * z[112] * z[385]) - l8 * z[6] * z[381]) - l8 * z[28] * z[380]) - l8 * z[30] * z[385]) - l8 * z[32] * z[384])) - 0.25 * mb * (((z[810] * z[319] + 2 * l2 * z[4] * z[319] + 2 * l2 * z[42] * z[320] + 2 * l3 * z[43] * z[318]) - z[811] * z[320]) - 2 * l4 * z[4] * z[318])
    z[874] = z[739] - z[779]
    z[875] = z[740] - z[786]
    z[876] = z[742] - z[816]
    z[877] = z[745] - z[823]
    z[878] = z[746] - z[829]
    z[879] = z[747] - z[833]
    z[880] = z[748] - z[836]

    # coef
    coef11 = -(z[772])
    coef12 = -(z[773])
    coef13 = -(z[774])
    coef14 = -(z[775])
    coef15 = -(z[776])
    coef16 = -(z[777])
    coef17 = -(z[778])
    coef21 = -(z[773])
    coef22 = -(z[780])
    coef23 = -(z[781])
    coef24 = -(z[782])
    coef25 = -(z[783])
    coef26 = -(z[784])
    coef27 = -(z[785])
    coef31 = -(z[774])
    coef32 = -(z[781])
    coef33 = -(z[804])
    coef34 = -(z[806])
    coef35 = -(z[807])
    coef36 = -(z[808])
    coef37 = -(z[809])
    coef41 = -(z[775])
    coef42 = -(z[782])
    coef43 = -(z[806])
    coef44 = -(z[818])
    coef45 = -(z[820])
    coef46 = -(z[821])
    coef47 = -(z[822])
    coef51 = -(z[776])
    coef52 = -(z[783])
    coef53 = -(z[807])
    coef54 = -(z[820])
    coef55 = -(z[826])
    coef56 = -(z[827])
    coef57 = -(z[828])
    coef61 = -(z[777])
    coef62 = -(z[784])
    coef63 = -(z[808])
    coef64 = -(z[821])
    coef65 = -(z[827])
    coef66 = -(z[831])
    coef67 = -(z[832])
    coef71 = -(z[778])
    coef72 = -(z[785])
    coef73 = -(z[809])
    coef74 = -(z[822])
    coef75 = -(z[828])
    coef76 = -(z[832])
    coef77 = -(z[835])

    # rhs
    rhs1 = -(z[874])
    rhs2 = -(z[875])
    rhs3 = -(z[876])
    rhs4 = -(z[877])
    rhs5 = -(z[878])
    rhs6 = -(z[879])
    rhs7 = -(z[880])

    # set up system of equations
    coef = @SMatrix [coef11 coef12 coef13 coef14 coef15 coef16 coef17; coef21 coef22 coef23 coef24 coef25 coef26 coef27; coef31 coef32 coef33 coef34 coef35 coef36 coef37; coef41 coef42 coef43 coef44 coef45 coef46 coef47; coef51 coef52 coef53 coef54 coef55 coef56 coef57; coef61 coef62 coef63 coef64 coef65 coef66 coef67; coef71 coef72 coef73 coef74 coef75 coef76 coef77]
    rhs = @SVector [rhs1, rhs2, rhs3, rhs4, rhs5, rhs6, rhs7]

    # derivatives
    q1p = u1
    q2p = u2
    q3p = u3
    q4p = u4
    q5p = u5
    q6p = u6
    q7p = u7
    @inbounds u1p, u2p, u3p, u4p, u5p, u6p, u7p = coef \ rhs

    return @SVector [q1p, q2p, q3p, q4p, q5p, q6p, q7p, u1p, u2p, u3p, u4p, u5p, u6p, u7p, he.cc.ω, ke.cc.ω, ae.cc.ω, hf.cc.ω, kf.cc.ω, af.cc.ω]
end


text:
"""
Contains functions for generating matrices with correlated or motif structure
"""

using LinearAlgebra
import Random.MersenneTwister

#------------------------------------------------------------------------------

"""
	generate_nv(N::Int64, T::Type; seed=nothing)

Returns a vector with entries i.i.d. on a zero-mean, unit-variance Gaussian
"""
function generate_nv(N::Int64, T::Type; seed=nothing)

	@assert T<:AbstractFloat "T must be a floating point type"
	rng = MersenneTwister(seed)

	return randn(rng, T, N)
end

"""
	generate_ONB(N::Int64, T::Type; seed=nothing)

Returns a random unitary matrix
"""
function generate_ONB(N::Int64, T::Type; seed=nothing)

	@assert T<:AbstractFloat "T must be a floating point type"
	rng = MersenneTwister(seed)
	Z = randn(rng, T, N, N)
	F = qr!(Z)

	return F.Q * diagm(0 => sign.(diag(F.R)))
end

"""
	generate_J(N::Int64, eta::T; seed=nothing) where T <: AbstractFloat

Returns a `N x N` Gaussian matrix with correlation between entries opposite the
diagonal equal to `eta` and diagonal entries themselves equal to `0`.

Entries of the matrix are scaled `sqrt(N)`.
"""
function generate_J(N::Int64, eta::T; seed=nothing) where T <: AbstractFloat

	rng = MersenneTwister(seed)

	J = zeros(T, N, N)
	M = [convert(T, 1)     convert(T, 0);
		 eta               sqrt(1-eta^2)]
	for j in 1:N-1
		for i in j+1:N
			(J[i,j], J[j,i]) = M * randn(rng, T, 2)
		end
	end

	return J ./ convert(T, sqrt(N))
end

"""
	generate_motif_P(N::Int64, kcon::T, kdiv::T, kchn::T; seed=nothing,
		no_autapses=true) where T <: AbstractFloat

Returns a rank-2 perturbation matrix that introduces converging, diverging, and
chain motifs when added to a dense random Gaussian matrix.

In dense Gaussian networks, motifs are interpreted as covariances between
synapses that share presynaptic or postsynaptic units. Use care when adding
perurbations to partially symmetric or antisymmetric matrices, as the final
symmetry depends on both the eta and kchn parameters.
"""
function generate_motif_P(N::Int64, kcon::T, kdiv::T, kchn::T; seed=nothing,
		no_autapses=true) where T <: AbstractFloat

	@assert kchn^2<=kcon*kdiv "κchn^2 cannot exceed κdiv*κcon"
	@assert kcon>=0 "Negative κcon is not allowed"
	@assert kdiv>=0 "Negative κdiv is not allowed"
	@assert kchn>=0 "Negative κchn is not allowed"

	rng = MersenneTwister(seed)

	if kcon == 0
		M = randn(rng, T, 2, N)
		a = zeros(N)
		b = sqrt(kdiv) .* reshape(M[2,:], (1, N))
	else
		M = [sqrt(kcon)          convert(T, 0);
			 kchn/sqrt(kcon)     sqrt(max(0, kdiv - kchn^2/kcon))]
		M *= randn(rng, T, 2, N)
		a = M[1,:]
		b = reshape(M[2,:], (1, N))
	end

	P = a .+ b

	if no_autapses
		P[diagind(P)] .= 0
	end

	return P ./ convert(T, sqrt(N))
end

"""
	generate_outlier_P(N::Int64, re::T, im::T, alpha::T; seed=nothing,
		no_autapses=true) where T <: AbstractFloat

Returns a rank-2 perturbation matrix that, when added to a random Gaussian
matrix, introduces a conjugate pair of outlier eigenvalues with given expected
real and imaginary parts.

The `alpha` parameter controls the strength of the uncorrelated components of
the outer product perturbations. This does not affect the expected value of the
real and imaginary parts.
"""
function generate_outlier_P(N::Int64, re::T, im::T, alpha::T; seed=nothing,
		no_autapses=true) where T <: AbstractFloat

	@assert re>0 "Real part of outliers must be positive"

	gamma = 2 * sqrt(1 /(1 + im^2 / re^2))
	rho = sqrt(2 * re / gamma)

	rng = MersenneTwister(seed)
	x1 = randn(rng, T, N)
	x2 = randn(rng, T, N)
	x3 = randn(rng, T, N)
	x4 = randn(rng, T, N)
	y1 = randn(rng, T, N)
	y2 = randn(rng, T, N)

	m1 = alpha .* x1 .+ rho .* y1
	m2 = alpha .* x2 .+ rho .* y2
	n1 = alpha .* x3 .+ rho .* y2 .+ rho * gamma .* y1
	n2 = alpha .* x4 .- rho .* y1

	P = m1 * n1' .+ m2 * n2'

	if no_autapses
		P[diagind(P)] .= 0
	end

	return P ./ N
end

ext: jl


text:
using FlowOverCircle
using FluxTraining
using Test

@testset "FlowOverCircle" begin
    learner = FlowOverCircle.train(epochs=5)
    loss = learner.cbstate.metricsepoch[ValidationPhase()][:Loss].values[end]
    @test loss < 0.1
end


text:
dispersion(obj::EconometricModel{<:LinearModelEstimators}) = true
dispersion(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}}) = false
isiv(obj::EconometricModel{<:LinearModelEstimators}) = obj.iv > 0
coef(obj::EconometricModel) = obj.β
coefnames(obj::EconometricModel) = obj.vars
# StatsBase.confint(obj::StatisticalModel) = error("coefint is not defined for $(typeof(obj)).")
deviance(obj::EconometricModel) = -2loglikelihood(obj)
deviance(obj::EconometricModel{<:LinearModelEstimators}) =
	weights(obj) |>
	(wts -> isa(wts, FrequencyWeights) ?
			sum(w * (y - ŷ)^2 for (w, y, ŷ) ∈ zip(wts, response(obj), fitted(obj))) :
			sum((y - ŷ)^2 for (y, ŷ) ∈ zip(response(obj), fitted(obj))))
hasintercept(obj::EconometricModel) = true
hasintercept(obj::EconometricModel{<:LinearModelEstimators}) =
	!any(t -> isa(t, InterceptTerm{false}), terms(obj.f.rhs))
islinear(obj::EconometricModel{<:LinearModelEstimators}) =
	!(isa(obj.estimator, NominalResponse) || isa(obj.estimator, OrdinalResponse))
nulldeviance(obj::EconometricModel{<:LinearModelEstimators}) =
	meanresponse(obj) |> (ȳ -> sum((y - ȳ)^2 for y ∈ response(obj)))
loglikelihood(obj::EconometricModel{<:LinearModelEstimators}) =
	√(deviance(obj) / dof_residual(obj)) |>
	(ϕ -> (sum(w * logpdf(Normal(μ, ϕ), y)
		   for (w, μ, y) ∈ zip(weights(obj), response(obj), fitted(obj)))))
function loglikelihood(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}})
	y = response(obj)
	categories = levels(y)
	b = mapreduce(elem -> (eachindex(categories) .== elem)', vcat, y)
	μ = predict(obj)
	wts = weights(obj)
	sum(wts[idx[1]] * logpdf(Categorical(collect(μ[idx[1],:])), idx[2]) for idx in findall(b))
end
function nullloglikelihood(obj::EconometricModel{<:LinearModelEstimators})
	ϕ = √(nulldeviance(obj) / (dof_residual(obj) + (length(coef(obj)) - hasintercept(obj))))
	μ = meanresponse(obj)
	sum(w * logpdf(Normal(μ, ϕ), y) for (w, y) ∈ zip(weights(obj), response(obj)))
end
function nullloglikelihood(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}})
	y = response(obj)
	categories = levels(y)
	b = mapreduce(elem -> (eachindex(categories) .== elem)', vcat, y)
	μ = collect(vec(mean(b, dims = 1)))
	wts = weights(obj)
	sum(wᵢ * logpdf(Categorical(μ), yᵢ) for (yᵢ, wᵢ) ∈ zip(y, wts))
end
# StatsBase.score(obj::StatisticalModel) = error("score is not defined for $(typeof(obj)).")
nobs(obj::EconometricModel) = sum(obj.w)
dof(obj::EconometricModel{<:LinearModelEstimators}) =
	length(coef(obj)) + dispersion(obj) + obj.iv
dof(obj::EconometricModel{<:ContinuousResponse}) =
	length(coef(obj)) + dispersion(obj) + obj.iv +
	(obj.estimator.groups |>
	 (g -> isempty(g) ? 0 : sum(length(group) - 1 for group ∈ g)))
dof(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}}) =
	length(coef(obj)) + dispersion(obj)
dof_residual(obj::EconometricModel{<:LinearModelEstimators}) =
	nobs(obj) - dof(obj) + dispersion(obj) + hasintercept(obj.f)
# Verify the dof_residual(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}})
dof_residual(obj::EconometricModel{<:Union{NominalResponse,OrdinalResponse}}) =
	nobs(obj) - dof(obj) + dispersion(obj) + (length(obj.estimator.categories) - 1)
mss(obj::EconometricModel{<:LinearModelEstimators}) =
	meanresponse(obj) |> (ȳ -> sum(w * abs2(ŷ - ȳ) for (ŷ, w) ∈ zip(fitted(obj), weights(obj))))
rss(obj::EconometricModel{<:LinearModelEstimators}) =
	residuals(obj) |> (û -> û' * Diagonal(weights(obj)) * û)
r2(obj::EconometricModel) = 1 - loglikelihood(obj) / nullloglikelihood(obj)
r2(obj::EconometricModel{<:LinearModelEstimators}) =
	(isiv(obj) || !hasintercept(obj)) ? NaN : rss(obj) |> (r -> 1 - r / (r + mss(obj)))
function adjr2(obj::EconometricModel{<:LinearModelEstimators})
	(isiv(obj) || !hasintercept(obj)) && return(NaN)
	ℓℓ = loglikelihood(obj)
	ℓℓ₀ = nullloglikelihood(obj)
	k = dof(obj)
	if variant == :McFadden
		1 - (ll - k)/ll0
	else
		error(":McFadden is the only currently supported variant")
	end
end
informationmatrix(obj::EconometricModel; expected::Bool = true) = obj.Ψ
vcov(obj::EconometricModel) = informationmatrix(obj)
vcov(obj::EconometricModel{<:LinearModelEstimators}) =
	deviance(obj) / dof_residual(obj) * informationmatrix(obj)
stderror(obj::EconometricModel) = sqrt.(diag(vcov(obj)))
confint(obj::EconometricModel, α = 0.05) =
    stderror(obj) * quantile(TDist(dof_residual(obj)), 1 - α / 2) |>
    (σ -> coef(obj) |>
        (β -> hcat(β .- σ, β .+ σ)))
weights(obj::EconometricModel) = obj.w
isfitted(obj::EconometricModel) = !isempty(obj.Ψ)
fitted(obj::EconometricModel) = obj.ŷ
response(obj::EconometricModel) = obj.y
meanresponse(obj::EconometricModel{<:LinearModelEstimators}) = mean(response(obj), weights(obj))
modelmatrix(obj::EconometricModel) = obj.X
# leverage(obj::RegressionModel) = error("leverage is not defined for $(typeof(obj)).")
residuals(obj::EconometricModel{<:LinearModelEstimators}) = response(obj) - fitted(obj)
function residuals(obj::EconometricModel{<:NominalResponse})
	@unpack y = obj
	@unpack categories = obj.estimator
	# Change from levels to the categories after StatsModels patch
	l = levels(y)
	y = [ findfirst(isequal(x), l) for x ∈ y ]
	b = mapreduce(elem -> (eachindex(categories) .== elem)', vcat, y)
	b - predict(obj)
end
predict(obj::EconometricModel) = fitted(obj)
predict(obj::EconometricModel{<:NominalResponse}) =
	mapslices(softmax, fitted(obj), dims = 2)
function predict(obj::EconometricModel{<:OrdinalResponse})
	y = response(obj)
	ŷ = fitted(obj)
	β = coef(obj)
	outcomes = obj.estimator.categories
	ζ = vcat(-Inf, β[end - length(outcomes) + 2:end], Inf)
	mapreduce(row -> (cdf.(Logistic(), ζ[2:end] .- ŷ[row]) .-
					  cdf.(Logistic(), ζ[1:end - 1] .- ŷ[row]))',
			  vcat,
			  eachindex(y))
end
function coeftable(obj::EconometricModel)
    β = coef(obj)
    σ = stderror(obj)
    t = β ./ σ
    p = 2ccdf.(TDist(dof_residual(obj)), abs.(t))
    mat = hcat(β, σ, t, p, confint(obj))
    colnms = ["PE   ", "SE   ", "t-value", "Pr > |t|", "2.5%", "97.5%"]
    rownms = obj.vars[2]
    CoefTable(mat, colnms, rownms, 4)
end
function coeftable(obj::EconometricModel{<:NominalResponse})
    β = coef(obj)
    σ = stderror(obj)
    t = β ./ σ
    p = 2ccdf.(TDist(dof_residual(obj)), abs.(t))
    mat = hcat(β, σ, t, p, confint(obj))
    colnms = ["PE   ", "SE   ", "t-value", "Pr > |t|", "2.5%", "97.5%"]
	vars = coefnames(obj)
    rownms = [ string(lhs, " ~ ", rhs) for lhs ∈ vars[1] for rhs ∈ vars[2] ]
    CoefTable(mat, colnms, rownms, 4)
end
function coeftable(obj::EconometricModel{<:OrdinalResponse})
    β = coef(obj)
    σ = stderror(obj)
    t = β ./ σ
    p = 2ccdf.(TDist(dof_residual(obj)), abs.(t))
    mat = hcat(β, σ, t, p, confint(obj))
    colnms = ["PE   ", "SE   ", "t-value", "Pr > |t|", "2.5%", "97.5%"]
	vars = coefnames(obj)
	outcomes = obj.estimator.categories
	rownms = vcat(vars[2],
				  [ string("(Intercept): ", outcomes[i], " | ", outcomes[i + 1])
				    for i ∈ 1:length(outcomes) - 1 ])
    CoefTable(mat, colnms, rownms, 4)
end


text:
function createMap(filename)
    values = readlines(filename*".txt")
    coords = Tuple[]
    instructions = String[]
    maxX = maxY = 1
    for line in values
        if !isempty(line) && line[1] == 'f'
            inst = split(line, " ")
            instructions = push!(instructions, inst[3])
        elseif !isempty(line)
            st = split(line, ",")
            x = parse(Int, st[1])
            y = parse(Int, st[2])
            juliacoords = (y+1, x+1)
            if juliacoords[1] > maxX
                maxX = juliacoords[1]
            end
            if juliacoords[2] > maxY
                maxY = juliacoords[2]
            end
            coords = push!(coords, juliacoords)
        end
    end
    grid = zeros(Int, maxX, maxY)
    for c in coords
        grid[c[1], c[2]] = 1
    end
    return grid, instructions
end

function foldH(loc, grid)
    cols = size(grid)[2]
    rows = loc - 1
    res = zeros(Int, rows, cols)
    top = grid[1:loc-1, :]
    low = grid[loc+1:size(grid)[1], :]
    for i in 1:rows
        for j in 1:cols
            if top[i, j] == 1 || low[rows+1-i, j] == 1
                res[i,j] = 1
            end
        end
    end
    return res
end

function foldV(loc, grid)
    cols = loc - 1
    rows = size(grid)[1]
    res = zeros(Int, rows, cols)
    left = grid[:, 1:loc-1]
    right = grid[:, loc+1:size(grid)[2]]
    for i in 1:rows
        for j in 1:cols
            if left[i,j] == 1 || right[i, cols+1-j] == 1
                res[i,j] = 1
            end
        end
    end
    return res
end

function solve(filename)
    grid, inst = createMap(filename)
    #println(grid)
    inst1 = inst[1]
    if inst1[1] == 'x'
        loc = parse(Int, inst1[3:length(inst1)]) + 1
        grid = foldV(loc, grid)
    else
        loc = parse(Int, inst1[3:length(inst1)]) + 1
        grid = foldH(loc, grid)
    end
    ndots = count(x -> x == 1, grid)
    println("number of dots: "*string(ndots))
end

function visualize(grid)
    for i in 1:size(grid)[1]
        line = ""
        for j in 1:size(grid)[2]
            if grid[i,j] == 1
                line = line * "⬜"
            else
                line = line * "⬛"
            end
        end
        println(line)
    end
end

function solve2(filename)
    grid, inst = createMap(filename)
    for ins in inst
        if ins[1] == 'x'
            println(ins)
            loc = parse(Int, ins[3:length(ins)]) + 1
            grid = foldV(loc, grid)
        else
            println(ins)
            loc = parse(Int, ins[3:length(ins)]) + 1
            grid = foldH(loc, grid)
        end
    end
    visualize(grid)
end
#solve("example")
solve("input")

solve2("input")


text:
immutable NullToNAPipe <: AbstractPipe
    F::AbstractPipe
    input::Input
    output::Output

    NullToNAPipe(F::AbstractPipe) =
        begin
            if isnonempty(F)
                return F
            end
            F = OptPipe(F)
            output = Output(Union{odomain(F),NAType}, omode(F))
            output = Output(output, ltotal=true, runique=false)
            return new(F, input(F), output)
        end
end

show(io::IO, pipe::NullToNAPipe) = show(io, pipe.F)

arms(pipe::NullToNAPipe) = AbstractPipe[pipe.F]

codegen(pipe::NullToNAPipe, X, I) =
    let T = odomain(pipe)
        @gensym Y
        quote
            $Y = $(codegen(pipe.F, X, I))
            isnull($Y) ? Iso{$T}(NA) : Iso{$T}(get($Y))
        end
    end

immutable DataFramePipe <: AbstractPipe
    F::AbstractPipe
    G::TuplePipe
    tags::Vector{Symbol}
    input::Input
    output::Output

    DataFramePipe(F::AbstractPipe, tagged_Gs::Vector{TaggedPipe}) =
        begin
            F = SeqPipe(F)
            G = TuplePipe(AbstractPipe[G0 for (tag, G0) in tagged_Gs])
            @assert(
                odomain(F) <: idomain(G),
                "$(repr(F)) and $(repr(G)) are not composable")
            tags = Symbol[tag for (tag, G0) in tagged_Gs]
            input = RBT.input(F >> G)
            output = Output(DataFrame)
            return new(F, G, tags, input, output)
        end
end

DataFramePipe(F, tagged_Gs) = DataFramePipe(F, TaggedPipe[TaggedPipe(tag, G) for (tag, G) in tagged_Gs])
DataFramePipe(F, tagged_Gs::Pair{Symbol}...) = DataFramePipe(F, tagged_Gs)

show(io::IO, pipe::DataFramePipe) =
    print(io, "DataFrame($(pipe.F), $(pipe.G))")

arms(pipe::DataFramePipe) = AbstractPipe[pipe.F, pipe.G]

codegen(pipe::DataFramePipe, X, I) =
    begin
        @gensym cols Y0 y idx
        Y = codegen_compose(pipe.G, pipe.F, X, I)
        return quote
            $Y0 = $Y
            $cols = Any[
                $([
                    let T = odomain(F)
                        isplural(F) ?
                            :( Vector{Vector{$T}}() ) :
                        ispartial(F) ?
                            :( DataVector{$T}([]) ) :
                            :( Vector{$T}() )
                    end
                    for F in pipe.G.Fs
                ]...)
            ]
            for $y in $Y0
                $([
                    isplural(F) ?
                        :( push!($cols[$idx], $y[$idx]) ) :
                    ispartial(F) ?
                        :( push!($cols[$idx], isnull($y[$idx]) ? NA : get($y[$idx])) ) :
                        :( push!($cols[$idx], $y[$idx]) )
                    for (idx, F) in enumerate(pipe.G.Fs)
                ]...)
            end
            Iso(DataFrame($cols, $(pipe.tags)))
        end
    end

compile(fn::Fn(:dataframe), base::Scope, flow::AbstractSyntax, ops::AbstractSyntax...) =
    let flow = compile(base, flow),
        ops = [compile(flow, op) for op in ops]
        isempty(ops) ? compile(fn, base, flow) : compile(fn, base, compile(Fn{:select}, base, flow, ops...))
    end

function compile(::Fn(:dataframe), base::Scope, flow::Query)
    flow = select(flow)
    if !isnull(flow.scope.items)
        fields = mkdffields(flow)
        pipe = DataFramePipe(flow.pipe, fields)
        scope = nest(base, DataFrame)
        if !isnull(flow.scope.tag)
            scope = settag(scope, get(flow.scope.tag))
        end
        return Query(scope, pipe)
    else
        if ispartial(flow)
            pipe = OptToNAPipe(flow.pipe)
            flow = Query(flow, pipe=pipe)
        end
    end
    return flow
end

function mkdffields(flow::Query)
    fields = ()
    for item in get(flow.scope.items)
        field = compile(item, flow.scope)
        if !issingular(field)
            field = compile(Fn{:dataframe}, flow.scope, field)
        end
        if !isnull(field.scope.items)
            fields = (fields..., mkdffields(field)...)
        else
            name = get(field.scope.tag, Symbol(""))
            F = field.pipe
            fields = (fields..., (name, F))
        end
    end
    return fields
end

function formatdf()
    @eval format(base::Scope, q::Query) =
        Query(compile(Fn{:dataframe}, base, q), syntax=q.syntax)
end

function Base.writemime(io::IO,
                        mime::MIME"text/html",
                        df::DataFrame,
                        sz::Int=0)
    n = size(df, 1)
    cnames = DataFrames._names(df)
    write(io, "
")
    write(io, "
")
    write(io, "	")
    for column_name in cnames
        write(io, "$column_name	")
    end
    write(io, "
")
    if sz <= 0
        tty_rows, tty_cols = Base.tty_size()
        sz = tty_rows
    end
    mxrow = min(n,sz)
    for row in 1:mxrow
        write(io, "
")
        write(io, "$row	")
        for column_name in cnames
            cell = df[row, column_name]
            if isa(cell, DataFrames.DataFrame)
                write(io, "")
                writemime(io, mime, cell, 1+round(Int, sqrt(sz)))
                write(io, "	")
            else
                write(io, "$(DataFrames.html_escape(string(cell)))	")
            end
        end
        write(io, "
")
    end
    if n > mxrow
        write(io, "
")
        write(io, "⋮	")
        for column_name in cnames
            write(io, "⋮	")
        end
        write(io, "
")
    end
    write(io, "
")
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# Test Linear System Optimizaiton
using SwitchTimeOpt
using MathOptInterface
using LinearAlgebra
const MOI = MathOptInterface
using Ipopt

# Define test options
objtol = 1e-04
primaltol = 1e-03

# Define Solver options
maxiter = 50
maxtime = 100.0
verbose = 0
tolerance = 1e-06

#Define solver
solver = Ipopt.Optimizer()
MOI.set(solver, MOI.RawOptimizerAttribute("tol"), tolerance)
MOI.set(solver, MOI.RawOptimizerAttribute("print_level"), verbose)
MOI.set(solver, MOI.RawOptimizerAttribute("max_cpu_time"), maxtime)
MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), maxiter)

# Size of the state space
nx = 2;

# Cost function Matrix
Q = 1.0 * Matrix(I, nx, nx)

# Initial State
x0 = [1.0; 1.0]

# Define initial and final time
t0 = 0.0
tf = 1.0

### Define System Dynamics
A = zeros(nx, nx, 2)
A[:, :, 1] = [-1 0;
              1  2]
A[:, :, 2] = [1 1;
              1 -2]

m = stoproblem(x0, A, solver=solver)
solve!(m)

@testset "Test optimal switching time" begin
  @test isapprox(gettau(m)[1], 0.26486646235103123, atol=primaltol)
end

@testset "Test status and optimal objective value" begin
  @test isapprox(getobjval(m), 5.2545429449272145, atol=objtol)
  @test string(getstat(m)) == "LOCALLY_SOLVED"
end


text:
import Statistics:mean

abstract type MultiObjectiveTestFunction{N, M} <: TestFunction{N,M} end

"""
    -A≤x≤A.
Values of `A` from 10 to 10⁵ have been used successfully.
Higher values of `A` increase the difficulty of the problem.

# References
[1] Schaffer, J. David (1984). Multiple Objective Optimization with Vector Evaluated Genetic Algorithms. Proceedings of the First Int. Conference on Genetic Algortihms, Ed. G.J.E Grefensette, J.J. Lawrence Erlbraum (PhD). Vanderbilt University. OCLC 20004572
[2] https://en.wikipedia.org/wiki/Test_functions_for_optimization#Test_functions_for_multi-objective_optimization
"""
struct Schaffer <: MultiObjectiveTestFunction{1,2} end

(f::Schaffer)(x::AbstractVector) = [x[1]^2, (x[1]-2)^2]

lowerbounds(::Schaffer) = Float64[-10.]
upperbounds(::Schaffer) = Float64[10.]

"""
# References
[1] Zitzler, E., Deb, K., & Thiele, L. (2000). Comparison of multiobjective evolutionary algorithms: Empirical results. Evolutionary Computation, 8(2), 173–195.
"""
struct ZDT1{N} <: MultiObjectiveTestFunction{N,2} end
function (f::ZDT1{N})(x::AbstractVector) where {N}
    f = Vector{Float64}(undef,2)
    f[1] = x[1]
    g = 1.0 + 9.0*mean(x[2:N])
    h = 1.0 - sqrt(f[1]/g)
    f[2] = g*h
    f
end

lowerbounds(::ZDT1{N}) where {N} = zeros(N)
upperbounds(::ZDT1{N}) where {N} = ones(N)

"""
# References
[1] Zitzler, E., Deb, K., & Thiele, L. (2000). Comparison of multiobjective evolutionary algorithms: Empirical results. Evolutionary Computation, 8(2), 173–195.
"""
struct ZDT2{N} <: MultiObjectiveTestFunction{N,2} end

function (f::ZDT2{N})(x::AbstractVector) where {N}
    f = Vector{Float64}(undef,2)
    f[1] = x[1]
    g = 1.0 + 9.0*mean(x[2:N])
    h = 1.0 - (f[1]/g)^2
    f[2] = g*h
    f
end

lowerbounds(::ZDT2{N}) where {N} = zeros(N)
upperbounds(::ZDT2{N}) where {N} = ones(N)

"""
# References
[1] Zitzler, E., Deb, K., & Thiele, L. (2000). Comparison of multiobjective evolutionary algorithms: Empirical results. Evolutionary Computation, 8(2), 173–195.
"""
struct ZDT3{N} <: MultiObjectiveTestFunction{N,2} end

function (f::ZDT3{N})(x::AbstractVector) where {N}
    f = Vector{Float64}(undef,2)
    f[1] = x[1]
    g = 1.0 + 9.0*mean(x[2:N])
    h = 1.0 - sqrt(f[1]/g)-(f[1]/g)*sin(10π*f[1])
    f[2] = g*h
    f
end
lowerbounds(::ZDT3{N}) where {N} = zeros(N)
upperbounds(::ZDT3{N}) where {N} = ones(N)

function ZDT(prob_id::Int, param::Int)
    if param < 2
        error("invalid `param`=$param, must be greater than 1")
    end
    if prob_id == 1
        ZDT1{param}()
    elseif prob_id == 2
        ZDT2{param}()
    elseif prob_id == 3
        ZDT3{param}()
    else
        error("invalid `prob_id`=$prob_id, must be one of (1, 2, 3)")
    end
end


text:
precompile(genspharm, (Int, Int, Float64, Float64))
for LT in [ZeroTo{false}, UnitRange{Int}], MT in [ZeroTo{true}, FullRange{true}, UnitRange{Int}]
    precompile(VSHCache, (Type{Float64}, ML{LT, MT}))
    precompile(VSHCache, (Type{Float64}, LM{LT, MT}))
    lr = LT === UnitRange{Int} ? LT(0:0) : LT(0)
    mr = MT === UnitRange{Int} ? MT(0:0) : MT(0)
    C = VectorSphericalHarmonics.VSHCache(Float64, LM(lr, mr))
    precompile(genspharm!, (typeof(C), Float64, Float64))
    C = VectorSphericalHarmonics.VSHCache(Float64, ML(lr, mr))
    precompile(genspharm!, (typeof(C), Float64, Float64))
    for YT in [PB, Hansen, Irreducible], B in [Cartesian, Polar, SphericalCovariant, HelicityCovariant]
        precompile(VSHCache, (Type{Float64}, YT, B, ML{LT, MT}))
        precompile(VSHCache, (Type{Float64}, YT, B, LM{LT, MT}))
        C = VectorSphericalHarmonics.VSHCache(Float64, YT(), B(), LM(lr, mr))
        precompile(vshbasis!, (typeof(C), YT, B, Float64, Float64))
        C = VectorSphericalHarmonics.VSHCache(Float64, YT(), B(), ML(lr, mr))
        precompile(vshbasis!, (typeof(C), YT, B, Float64, Float64))
    end
end


text:
module TestGilljamRewire
    using BioEnergeticFoodWebs
    using Test
    using Random

    Random.seed!(1)
    S = 3
    A = [0 0 0 ;
         1 0 0 ;
         1 1 1 ]
    biomass = [1.0,1.0,1.0]
    parameters = model_parameters(A, rewire_method = :DS)

#testing Gilljam no extinctions
    GilljamTest = BioEnergeticFoodWebs.Gilljam(S,parameters,biomass)
    @test GilljamTest[1] == A
    @test GilljamTest[2] == parameters
#Testing with extinction
    parameters[:extinctions] = 1
    GilljamTest = BioEnergeticFoodWebs.Gilljam(S,parameters,biomass)
    @test GilljamTest[1] == [0 0 0; 0 1 0; 0 1 1]
    parameters[:costMat][2,2] = 0.0
    @test GilljamTest[2] == parameters
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
using Gadfly
p = plot(y = [1, 2, 3, 4, 5])
img = SVG("plot3.svg", 8cm, 8cm)
draw(img, p)
println("Julia 👩🦰 sagt: I did it, Babe!")

max_stars_repo_name: kantel/julia

max_issues_repo_name: kantel/julia

max_forks_repo_name: kantel/julia


text:
using Oceananigans.Fields: ComputedField

PressureField(model, data=model.pressures.pHY′.data, recompute_safely=false) =
    ComputedField(model.pressures.pHY′ + model.pressures.pNHS, data=data,
                      recompute_safely=recompute_safely)