text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -2.9922365365158297 <= q <= 1.3959867501971397)

                     add_NL_constraint(m, :((-0.7172866311917714 + 0.8530094465195428*(0.49975691401319544 + 0.817850285997455*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.15011340339301915*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.49975691401319544 + 0.817850285997455*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.15011340339301915*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2 + -0.2302663292873568*(0.6056562038214484 + 0.017530899400252764*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.9141090581590978*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.6056562038214484 + 0.017530899400252764*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.9141090581590978*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((-0.7172866311917714 + 0.8530094465195428*(0.49975691401319544 + 0.817850285997455*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.15011340339301915*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.49975691401319544 + 0.817850285997455*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.15011340339301915*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2 + -0.2302663292873568*(0.6056562038214484 + 0.017530899400252764*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.9141090581590978*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.6056562038214484 + 0.017530899400252764*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.9141090581590978*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2 + (0.11802092343800519 + -0.23198020992542823*(0.49975691401319544 + 0.817850285997455*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.15011340339301915*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.49975691401319544 + 0.817850285997455*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.15011340339301915*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2 + 0.3694431087223533*(0.6056562038214484 + 0.017530899400252764*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.9141090581590978*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.6056562038214484 + 0.017530899400252764*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.9141090581590978*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.11802092343800519 + -0.23198020992542823*(0.49975691401319544 + 0.817850285997455*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.15011340339301915*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.49975691401319544 + 0.817850285997455*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + -0.15011340339301915*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)/sqrt(2)))/2 + 0.3694431087223533*(0.6056562038214484 + 0.017530899400252764*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.9141090581590978*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2)*(1 + erf((0.6056562038214484 + 0.017530899400252764*(-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.9044530867886094 + -0.3504940249833268*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + -0.925387475816613*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + 0.9141090581590978*(-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.49866226232151334 + 0.31623027699891626*(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))*(1 + erf((0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))/sqrt(2)))/2 + 0.1823838704363938*(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))*(1 + erf((-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))/sqrt(2)))/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:
s1 = """
123456
abcde
ABCDEFG
"""
d1 = Document("untitled", s1, false)
@test get_text(d1) == s1
@test get_offset(d1, 0, 4) == 4
@test get_offset(d1, 1, 2) == 9
@test get_line_offsets(d1) == [0, 7, 13, 21]
@test get_position_at(d1, 1) == (0, 1)
@test get_position_at(d1, 8) == (1, 1)
@test get_position_at(d1, 15) == (2, 2)
@test get_open_in_editor(d1) == false
set_open_in_editor(d1, true)
@test get_open_in_editor(d1) == true
set_open_in_editor(d1, false)
@test get_open_in_editor(d1) == false
@test is_workspace_file(d1) == false

s2 = """
12μ456
abηde
ABCDEFG"""
d2 = Document("untitled", s2, true)
@test get_offset(d2, 0, 4) == 5
@test get_offset(d2, 1, 2) == 10
@test get_line_offsets(d2) == [0, 8, 15]
@test get_position_at(d2, 1) == (0, 1)
@test get_position_at(d2, 9) == (1, 1)
@test get_position_at(d2, 17) == (2, 2)
@test is_workspace_file(d2) == true

applytextdocumentchanges(d2, LanguageServer.TextDocumentContentChangeEvent(Range(1), 0, "12"))
@test get_line_offsets(d2) == [0, 8, 17]

s4 = "1234\r\nabcd"
d4 = Document("untitled", s4, false)
@test_broken get_line_offsets(d4) == [0, 5]

s5 = "1234\nabcd\n"
d5 = Document("untitled", s5, false)
@test get_line_offsets(d5) == [0, 5, 10]

s6 = "\n"
d6 = Document("untitled", s6, false)
@test get_line_offsets(d6) == [0,1]

@testset "applytextdocumentchanges" begin
    doc = LS.Document("file:///example/path/example.jl", "function foo()", false)
    c1 = LS.TextDocumentContentChangeEvent(LS.Range(LS.Position(0,14), LS.Position(0,14)),
                                        0, "\n")
    c2 = LS.TextDocumentContentChangeEvent(LS.Range(LS.Position(1,0), LS.Position(1,0)),
                                           0, "    ")
    c3 = LS.TextDocumentContentChangeEvent(missing, missing, "println(\"Hello World\")")

    LS.applytextdocumentchanges(doc, c1)
    @test LS.get_text(doc) == "function foo()\n"
    # Implicitly test for issue #403
    LS.applytextdocumentchanges(doc, c2)
    @test LS.get_text(doc) == "function foo()\n    "
    LS.applytextdocumentchanges(doc, c3)
    @test LS.get_text(doc) == "println(\"Hello World\")"
    # doc currently has only one line, applying change to 2nd line should throw
    @test_throws ErrorException LS.applytextdocumentchanges(doc, c2)
end

@testset "UTF16 handling" begin
    doc = LanguageServer.Document("", "aaa", false)
    @test sizeof(LanguageServer.get_text(doc)) == 3
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_at(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 1) == 1
    @test LanguageServer.get_position_at(doc, 1) == (0, 1)
    @test LanguageServer.get_offset(doc, 0, 2) == 2
    @test LanguageServer.get_position_at(doc, 2) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 3) == 3
    @test LanguageServer.get_position_at(doc, 3) == (0, 3)

    doc = LanguageServer.Document("", "ααα", false)
    @test sizeof(LanguageServer.get_text(doc)) == 6
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_at(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 1) == 1
    @test LanguageServer.get_position_at(doc, 1) == (0, 1)
    @test LanguageServer.get_offset(doc, 0, 2) == 3
    @test LanguageServer.get_position_at(doc, 3) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 3) == 5
    @test LanguageServer.get_position_at(doc, 5) == (0, 3)

    doc = LanguageServer.Document("", "ࠀࠀࠀ", false) # 0x0800
    @test sizeof(LanguageServer.get_text(doc)) == 9
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_at(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 1) == 1
    @test LanguageServer.get_position_at(doc, 1) == (0, 1)
    @test LanguageServer.get_offset(doc, 0, 2) == 4
    @test LanguageServer.get_position_at(doc, 4) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 3) == 7
    @test LanguageServer.get_position_at(doc, 7) == (0, 3)

    doc = LanguageServer.Document("", "𐐀𐐀𐐀", false)
    @test sizeof(LanguageServer.get_text(doc)) == 12
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_at(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 2) == 1
    @test LanguageServer.get_position_at(doc, 1) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 4) == 5
    @test LanguageServer.get_position_at(doc, 5) == (0, 4)
    @test LanguageServer.get_offset(doc, 0, 6) == 9
    @test LanguageServer.get_position_at(doc, 9) == (0, 6)

    doc = LanguageServer.Document("", "𐀀𐀀𐀀", false) # 0x010000
    @test sizeof(LanguageServer.get_text(doc)) == 12
    @test LanguageServer.get_offset(doc, 0, 0) == 0
    @test LanguageServer.get_position_at(doc, 0) == (0, 0)
    @test LanguageServer.get_offset(doc, 0, 2) == 1
    @test LanguageServer.get_position_at(doc, 1) == (0, 2)
    @test LanguageServer.get_offset(doc, 0, 4) == 5
    @test LanguageServer.get_position_at(doc, 5) == (0, 4)
    @test LanguageServer.get_offset(doc, 0, 6) == 9
    @test LanguageServer.get_position_at(doc, 9) == (0, 6)
end

ext: jl


text:
"""
$(SIGNATURES)

Compute the centroid or geometric center of the given `ConvexHull` using the
formulas given [here](https://en.wikipedia.org/wiki/Centroid#Of_a_polygon).
"""
function centroid(hull::ConvexHull)
    T = eltype(hull)
    R = arithmetic_closure(T)
    vertices = hull.vertices
    n = length(vertices)
    @inbounds begin
        if n === 0
            error()
        elseif n === 1
            return R.(vertices[1])
        elseif n === 2
            return (vertices[1] + vertices[2]) / 2
        else
            c = cross2(vertices[n], vertices[1])
            centroid = (vertices[n] + vertices[1]) * c
            double_area = c
            @simd for i in Base.OneTo(n - 1)
                c = cross2(vertices[i], vertices[i + 1])
                centroid += (vertices[i] + vertices[i + 1]) * c
                double_area += c
            end
            double_area = abs(double_area)
            centroid /= 3 * double_area
            return centroid
        end
    end
end


text:
module lewd # labour estimators, wrongly done

using Distributions
using LinearAlgebra
using Infiltrator
using StatsBase
using DataFrames
using Plots
using CategoricalArrays
using ShiftedArrays
using Chain
using DataFramesMeta
using LightGraphs
# using TikzGraphs
using Optim
using Random

# Find all movers defined as those with max spell > 0
function find_movers(df::DataFrame)
    move_df = @chain df begin
       groupby(:w_id)
       transform(:spell => maximum)
       subset(:spell_maximum => x -> x .> 0)
    end

   return move_df
end

# Get the next firm the mover is moving to
# by leading the firm column and extracting
# just the firm and next firm as one observation
function find_firm_links(mover_df::DataFrame)
    firm_link_df = @chain mover_df begin
        sort([:w_id, :t])
        groupby(:w_id)
        transform(:j => lead => :j_next)
        transform([:j, :j_next] => .==)
        subset(:j_j_next_BroadcastFunction => x -> x .== false, skipmissing = true)
        select(:j, :j_next)
    end
    return firm_link_df
end

function create_M_flows(link_df)
    M = @chain link_df begin
        groupby([:j, :j_next])
        combine(nrow => :count)
        unique()
    end
    return M
end

# %%

function create_adjacency_matrix(firm_link_df::DataFrame, df::DataFrame; count = false)
    adjacency_matrix = zeros(Int, maximum(df.j), maximum(df.j))
    if count == false
        firm_link_df[!, "count"] .= 1
    end
    for firm_a in unique(df.j), firm_b in unique(df.j)
        subset_a_df = firm_link_df[(firm_link_df.j .== firm_a) .& (firm_link_df.j_next .== firm_b), :]
        if size(subset_a_df)[1] != 0
            adjacency_matrix[firm_b, firm_a] = subset_a_df.count[1]
        end
    end
    return adjacency_matrix
end

function create_fixed_point_matrices(M_df, df)
    sum_flows = @chain M_df begin
        groupby(:j)
        combine(:count => sum)
        sort(:j)
        unique()
    end
    S_kk = Diagonal(sum_flows.count_sum)
    M_0 = create_adjacency_matrix(M_df, df, count = true)
    return S_kk, M_0
end

function create_fixed_point_matrices(df)
    link_df = find_firm_links(find_movers(df))
    M_df = create_M_flows(link_df)

    M_0 = create_adjacency_matrix(M_df, df, count = true)
    S_kk = Diagonal(sum(M_0, dims = 1)[:])
    return S_kk, M_0
end

function create_adjacency_matrix(firm_link_df::DataFrame, df::DataFrame)
    adjacency_matrix = zeros(Int, maximum(df.j), maximum(df.j))
    for firm_a in unique(df.j), firm_b in unique(df.j)
        subset_a_df = firm_link_df[(firm_link_df.j .== firm_a) .& (firm_link_df.j_next .== firm_b), :]
        if size(subset_a_df)[1] != 0
            adjacency_matrix[firm_a, firm_b] = 1
            adjacency_matrix[firm_b, firm_a] = 1
        end
    end
    return adjacency_matrix
end

function create_connected_df(df::DataFrame)
    move_df = find_movers(df)
    firm_df = find_firm_links(move_df)
    adjacency_matrix = create_adjacency_matrix(firm_df, move_df)

    simple_graph = SimpleGraph(adjacency_matrix);

    connected_network = connected_components(simple_graph)

    connected_set = connected_network[1]

    println("We have only $(length(connected_set)) firms fully connected $(length(connected_set)/maximum(df.j)) of the market...") # Double check we might be wrong...

    df_connected = df[in(connected_set).(df.j),:]
    return df_connected
end

function estimate_rank(S_kk, M_0; tol= 1e-6)
    N_connected = size(S_kk, 1)
    S_inv = inv(S_kk)
    initial_V_EE = fill(0.5, N_connected)
    lhs = S_inv * M_0 * exp.(initial_V_EE)
    rhs = exp.(initial_V_EE)
    i = 0
    max_diff = Inf
    while max_diff > tol
        i += 1
        rhs = lhs
        lhs = S_inv*M_0 * rhs
        max_diff = maximum(abs.(lhs .- rhs))
        println("max diff: $max_diff")
    end
    return rhs
end
end


text:
#=
Included in logutils xxx.jl if HTTP is loaded
=#

"HTTP.Response already has a show method, we're not overwriting that.
This metod is called only when logging to an Abstractdevice. The default
show method does not print binary data well as per now."
function _show(d::AbstractDevice, response::HTTP.Messages.Response)
    _log(d, :green, "Response status: ", :bold, response.status," ")
    response.status > 0 && _log(d, HTTP.Messages.STATUS_MESSAGES[response.status], " ")
    if !isempty(response.headers)
        _log(d, :green, " Headers: ", :bold, length(response.headers))
        _log(d, :green, "\n", response.headers)
    end
    if isdefined(response, :cookies)
        if !isempty(response.cookies)
            _log(d, " Cookies: ", :bold, length(response.cookies))
            _log(d, "\n", response.cookies)
        end
    end
    if !isempty(response.body)
        _log(d, "\t", DataDispatch(response.body, HTTP.header(response, "content-type", "")))
    end
    nothing
end

"HTTP.Request already has a show method, we're not overwriting that.
This metod is called only when logging to an Abstractdevice"
function _show(d::AbstractDevice, request::HTTP.Messages.Request)
    _log(d,  :normal, :light_yellow, "Request ", :normal)
    _log(d, :bold, request.method, " ", :cyan, request.target, "\n", :normal)
    if !isempty(request.body)
        _log(d, "\t", DataDispatch(request.body, HTTP.header(request, "content-type", "")))
    end
    if !isempty(request.headers)
        _log(d, "\t", :cyan, " Headers: ", length(request.headers))
        _log(d, :cyan, "\n", request.headers)
    end
    nothing
end


text:
#
# Definitions to move to NDTensors
#

storage(T::Tensor) = store(T)

#function NDTensors.blockdims(inds::IndexSet, block::Tuple)
#  return ntuple(i -> blockdim(inds,block,i), ValLength(block))
#end


text:
using Dyno
using Base.Test

rootdir = Pkg.dir("Dyno")

include("test_model_import.jl")

include("test_trisystem.jl")


text:
# println("Johan/emission_comparison.jl")
# include("Johan/emission_comparison.jl")
# println("Johan/iron_line_comparison_approx.jl")
# include("Johan/iron_line_comparison_approx.jl")
# println("Johan/temperature_renders.jl")
# include("Johan/temperature_renders.jl")
# println("Others/ISCO_spin.jl")
# include("Others/ISCO_spin.jl")
# println("PT/iron_lines/iron_heatmap.jl")
# include("PT/iron_lines/iron_heatmap.jl")
# println("PT/iron_lines/iron_line_basic.jl")
# include("PT/iron_lines/iron_line_basic.jl")
# println("PT/iron_lines/iron_line_combined.jl")
# include("PT/iron_lines/iron_line_combined.jl")
println("PT/iron_lines/iron_line_comparison.jl")
include("PT/iron_lines/iron_line_comparison.jl")
println("PT/renders/renders_parameters.jl")
include("PT/renders/renders_parameters.jl")


text:
"""
    LayerScale(λ, planes::Int)

Creates a `Flux.Scale` layer that performs "`LayerScale`"
([reference](https://arxiv.org/abs/2103.17239)).

# Arguments
- `planes`: Size of channel dimension in the input.
- `λ`: initialisation value for the learnable diagonal matrix.
"""
LayerScale(planes::Int, λ) =
    λ > 0 ? Flux.Scale(fill(Float32(λ), planes), false) : identity

"""
    DropPath(p)

Implements Stochastic Depth - equivalent to `Dropout(p; dims = 4)` when `p` ≥ 0.
([reference](https://arxiv.org/abs/1603.09382))

# Arguments
- `p`: rate of Stochastic Depth.
"""
DropPath(p) = p ≥ 0 ? Dropout(p; dims = 4) : identity


text:
using Knet
using ArgParse
using JLD,MAT,JSON,HDF5

include("util.jl")
include("model.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "CNN-HMM hybrid ASR system in Knet by Ilker Kesen."

    @add_arg_table s begin
        ("--features"; required=true; help="MFSC features file")
        ("--jsonfile"; required=true)
        ("--model"; default=nothing)
        ("--savefile"; default=nothing)
        ("--atype"; default=(gpu()>=0?"KnetArray{Float32}":"Array{Float32}"))
        ("--window"; default=10; arg_type=Int64)
        ("--period"; default=20; arg_type=Int64)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    o[:atype] = eval(parse(o[:atype]))
    atype = o[:atype]
    # sr = o[:seed] > 0 ? srand(o[:seed]) : srand()

    # load data
    jsondata = JSON.parsefile(o[:jsonfile])
    p2i, i2p = build_vocab(jsondata)

    # make data
    trnids, tstids = [], []
    # trn = make_data(o[:features], jsondata; extra=trnids)
    tst = make_data(o[:features], jsondata, "test"; extra=tstids)

    # load model
    w = load(o[:model], "w")

    # convert atype
    w = map(wi->convert(atype,wi), w)

    solutions = Dict()

    data = tst
    getid(i) = tstids[i]
    for (i,sample) in enumerate(data)
        x = make_input(sample)
        # ygold = make_output(sample, p2i)
        ygold = map(si->si[2], sample)
        ygold = make_ygold(ygold)
        ypred = predict(w,convert(atype,x))
        ysoft = softmax(ypred)

        generation = generate(ysoft, o[:window])
        generation = map(t->i2p[t], generation)
        solutions[getid(i)] = Dict("ygold"=>ygold, "ypred"=>generation)

        # println("$i generation so far.")
        if i % o[:period] == 0
            println("$i generation so far ", now())
        end
    end

    ff = open(abspath(o[:savefile]), "w")
    write(ff, JSON.json(solutions))
    close(ff)
    println("done")
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/phoneme-convnet

max_issues_repo_name: ilkerkesen/phoneme-convnet

max_forks_repo_name: ilkerkesen/phoneme-convnet


text:
module TestReducePartitionBy
using ...Utils

_test() = include(joinpath(@__DIR__, "../examples/reduce_partition_by.jl"))

function test()
    try
        _test()
    finally
        Utils.unsafe_free_all!(@__MODULE__)
    end
end

end  # module


text:
# ____________________________________________________________________
# hostgalaxy
#
mutable struct hostgalaxy <: AbstractComponent
    norm::Parameter
    template::String
    base::Vector{Float64}

    function hostgalaxy(template::String)
        filename = qsfit_data() * "/swire/" * template * "_template_norm.sed"
        if !isfile(filename)
            filename = template
        end
        @assert isfile(filename) "File $filename do not exists"
        out = new(Parameter(1), filename, Vector{Float64}())
        out.norm.val = 1
        out.norm.low = 0
        return out
    end
end

function prepare!(comp::hostgalaxy, domain::Domain{1})
    d = readdlm(comp.template, comments=true)
    @assert typeof(d) == Matrix{Float64}
    itp = Spline1D(d[:,1], d[:,2], k=1, bc="error")
    comp.base = itp(domain[:])
    comp.base ./= itp(5500.)
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::hostgalaxy, domain::Domain{1},
                   norm)
    buffer .= norm .* comp.base
end


text:
# ~ bound checking functions ~
checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Integer) =
    (0 < k ≤ size(A, 1) && 0 < j ≤ size(A, 2) || throw(BoundsError(A, (k,j))))

checkbounds(A::AbstractDerivativeOperator, kr::AbstractRange, j::Integer) =
    (checkbounds(A, first(kr), j); checkbounds(A,  last(kr), j))

checkbounds(A::AbstractDerivativeOperator, k::Integer, jr::AbstractRange) =
    (checkbounds(A, k, first(jr)); checkbounds(A, k,  last(jr)))

checkbounds(A::AbstractDerivativeOperator, kr::AbstractRange, jr::AbstractRange) =
    (checkbounds(A, kr, first(jr)); checkbounds(A, kr,  last(jr)))

checkbounds(A::AbstractDerivativeOperator, k::Colon, j::Integer) =
    (0 < j ≤ size(A, 2) || throw(BoundsError(A, (size(A,1),j))))

checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Colon) =
    (0 < k ≤ size(A, 1) || throw(BoundsError(A, (k,size(A,2)))))

# BandedMatrix{A,B,C,D}(A::DerivativeOperator{A,B,C,D}) = BandedMatrix(convert(Array, A, A.stencil_length), A.stencil_length, div(A.stencil_length,2), div(A.stencil_length,2))

# ~~ getindex ~~
@inline function getindex(A::AbstractDerivativeOperator, i::Int, j::Int)
    @boundscheck checkbounds(A, i, j)
    mid = div(A.stencil_length, 2) + 1
    bpc = A.stencil_length - mid
    l = max(1, low(j, mid, bpc))
    h = min(A.stencil_length, high(j, mid, bpc, A.stencil_length, A.dimension))
    slen = h - l + 1
    if abs(i - j) > div(slen, 2)
        return 0
    else
        return A.stencil_coefs[mid + j - i]
    end
end

# scalar - colon - colon
@inline getindex(A::AbstractDerivativeOperator, kr::Colon, jr::Colon) = convert(Array,A)

@inline function getindex(A::AbstractDerivativeOperator, rc::Colon, j)
    T = eltype(A.stencil_coefs)
    v = zeros(T, A.dimension)
    v[j] = one(T)
    copyto!(v, A*v)
    return v
end

# symmetric right now
@inline function getindex(A::AbstractDerivativeOperator, i, cc::Colon)
    T = eltype(A.stencil_coefs)
    v = zeros(T, A.dimension)
    v[i] = one(T)
    copyto!(v, A*v)
    return v
end

# UnitRanges
@inline function getindex(A::AbstractDerivativeOperator, rng::UnitRange{Int}, cc::Colon)
    m = convert(Array,A)
    return m[rng, cc]
end

@inline function getindex(A::AbstractDerivativeOperator, rc::Colon, rng::UnitRange{Int})
    m = convert(Array,A)
    return m[rnd, cc]
end

@inline function getindex(A::AbstractDerivativeOperator, r::Int, rng::UnitRange{Int})
    m = A[r, :]
    return m[rng]
end

@inline function getindex(A::AbstractDerivativeOperator, rng::UnitRange{Int}, c::Int)
    m = A[:, c]
    return m[rng]
end

@inline function getindex(A::AbstractDerivativeOperator{T}, rng::UnitRange{Int}, cng::UnitRange{Int}) where T
    N = A.dimension
    if (rng[end] - rng[1]) > ((cng[end] - cng[1]))
        mat = zeros(T, (N, length(cng)))
        v = zeros(T, N)
        for i = cng
            v[i] = one(T)
            #=
                calculating the effect on a unit vector to get the matrix of transformation
                to get the vector in the new vector space.
            =#
            mul!(view(mat, :, i - cng[1] + 1), A, v)
            v[i] = zero(T)
        end
        return mat[rng, :]

    else
        mat = zeros(T, (length(rng), N))
        v = zeros(T, N)
        for i = rng
            v[i] = one(T)
            #=
                calculating the effect on a unit vector to get the matrix of transformation
                to get the vector in the new vector space.
            =#
            mul!(view(mat, i - rng[1] + 1, :), A, v)
            v[i] = zero(T)
        end
        return mat[:, cng]
    end
end

#=
    This the basic mul! function which is broken up into 3 parts ie. handling the left
    boundary, handling the interior and handling the right boundary. Finally the vector is
    scaled appropriately according to the derivative order and the degree of discretizaiton.
    We also update the time stamp of the DerivativeOperator inside to manage time dependent
    boundary conditions.
=#
function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::Union{DerivativeOperator{T},UpwindOperator{T}}, x::AbstractVector{T}) where T<:Real
    convolve_BC_left!(x_temp, x, A)
    convolve_interior!(x_temp, x, A)
    convolve_BC_right!(x_temp, x, A)
    rmul!(x_temp, @.(1/(A.dx^A.derivative_order)))
    A.t[] += 1 # incrementing the internal time stamp
end

#=
    This definition of the mul! function makes it possible to apply the LinearOperator on
    a matrix and not just a vector. It basically transforms the rows one at a time.
=#
function LinearAlgebra.mul!(x_temp::AbstractArray{T,2}, A::AbstractDerivativeOperator{T}, M::AbstractMatrix{T}) where T<:Real
    if size(x_temp) == reverse(size(M))
        for i = 1:size(M,1)
            mul!(view(x_temp,i,:), A, view(M,i,:))
        end
    else
        for i = 1:size(M,2)
            mul!(view(x_temp,:,i), A, view(M,:,i))
        end
    end
end

# Base.length(A::AbstractDerivativeOperator) = A.stencil_length
Base.ndims(A::AbstractDerivativeOperator) = 2
Base.size(A::AbstractDerivativeOperator) = (A.dimension, A.dimension)
Base.size(A::AbstractDerivativeOperator,i::Integer) = size(A)[i]
Base.length(A::AbstractDerivativeOperator) = reduce(*, size(A))

#=
    For the evenly spaced grid we have a symmetric matrix
=#
Base.transpose(A::Union{DerivativeOperator,UpwindOperator}) = A
Base.adjoint(A::Union{DerivativeOperator,UpwindOperator}) = A
LinearAlgebra.issymmetric(::Union{DerivativeOperator,UpwindOperator}) = true

#=
    This function ideally should give us a matrix which when multiplied with the input vector
    returns the derivative. But the presence of different boundary conditons complicates the
    situation. It is not possible to incorporate the full information of the boundary conditions
    in the matrix as they are additive in nature. So in this implementation we will just
    return the inner stencil of the matrix transformation. The user will have to manually
    incorporate the BCs at the ends.
=#
function Base.convert(::Type{Array}, A::AbstractDerivativeOperator{T}, N::Int=A.dimension) where T
    @assert N >= A.stencil_length # stencil must be able to fit in the matrix
    mat = zeros(T, (N, N))
    v = zeros(T, N)
    for i=1:N
        v[i] = one(T)
        #=
            calculating the effect on a unit vector to get the matrix of transformation
            to get the vector in the new vector space.
        =#
        mul!(view(mat,:,i), A, v)
        v[i] = zero(T)
    end
    return mat
end

#=
    This function ideally should give us a matrix which when multiplied with the input vector
    returns the derivative. But the presence of different boundary conditons complicates the
    situation. It is not possible to incorporate the full information of the boundary conditions
    in the matrix as they are additive in nature. So in this implementation we will just
    return the inner stencil of the matrix transformation. The user will have to manually
    incorporate the BCs at the ends.
=#
function SparseArrays.sparse(A::AbstractDerivativeOperator{T}) where T
    N = A.dimension
    mat = spzeros(T, N, N)
    v = zeros(T, N)
    row = zeros(T, N)
    for i=1:N
        v[i] = one(T)
        #=
            calculating the effect on a unit vector to get the matrix of transformation
            to get the vector in the new vector space.
        =#
        mul!(row, A, v)
        copyto!(view(mat,:,i), row)
        @. row = 0 * row;
        v[i] = zero(T)
    end
    return mat
end

#=
    Fallback methods that use the full representation of the operator
=#
Base.exp(A::AbstractDerivativeOperator{T}) where T = exp(convert(Array,A))
Base.:\(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A \ convert(Array,B)
Base.:\(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = convert(Array,A) \ B
Base.:/(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A / convert(Array,B)
Base.:/(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = convert(Array,A) / B


text:
using EconPDEs, Distributions

# Arbitrage With Holding Costs: A Utility-Based Approach
# Author(s): Bruce Tuckman and Jean-Luc Vila

struct ArbitrageHoldingCosts
    c::Float64
    r::Float64
    ρ::Float64
    σ::Float64
    a::Float64
    T::Float64
end

function ArbitrageHoldingCosts(;c = 0.06, r = 0.09,  ρ = 5.42, σ = 0.88, a = 26.72, T = 100)
    ArbitrageHoldingCosts(c, r, ρ, σ, a, T)
end

function initialize_stategrid(m::ArbitrageHoldingCosts; n = 200)
    d = Normal(0, sqrt(m.σ^2 / (2 * m.ρ)))
    zmin = quantile(d, 0.00001)
    zmax = quantile(d, 0.99999)
    OrderedDict(:z => range(zmin, stop = zmax, length = n))
end

function initialize_y(m::ArbitrageHoldingCosts, stategrid::OrderedDict)
    OrderedDict(:F => zeros(length(stategrid[:z])))
end

function (m::ArbitrageHoldingCosts)(state::NamedTuple, y::NamedTuple, τ::Number)
    c = m.c ; r = m.r ; ρ = m.ρ ; σ = m.σ ; a = m.a ; T = m.T
    z = state.z
    F, Fz, Fzz = y.F, y.Fz, y.Fzz
    ϕ = z * (1 + z^2)^(-1/2)
    ϕz = (1 + z^2)^(-3/2)
    ϕzz = - 3 * z * (1 + z^2)^(-5/2)
    sτ = c / r * (1 - exp(- r * (T - τ)))
    μL = ρ * z - 0.5 * σ^2 * ϕzz / ϕz + c / sτ * (ϕ - 1) / ϕz
    iL = 1 / a * (μL / σ^2 + Fz)
    μS = ρ * z - 0.5 * σ^2 * ϕzz / ϕz + c / sτ * (ϕ + 1) / ϕz
    iS = 1 / a * (μS / σ^2 + Fz)
    μ, i  = 0.0, 0.0
    if iL > 0
        i = iL
        μ = μL
    elseif iS < 0
        i = iS
        μ = μS
    end
    if z >= 0
        i_myopic = max(0, i - Fz / a)
    else
        i_myopic = min(0, i - Fz / a)
    end
    # otherwise i = 0.0 and value of μ does not matter
    Ft = (μ + σ^2 * Fz) * a * i - 0.5 * σ^2 * (a * i)^2 - ρ * z * Fz + 0.5 * σ^2 * (Fzz - Fz^2)
    return (Ft,), (-ρ * z,), (F = F, i = i, I = i / (sτ * ϕz * exp(r * (T - τ))), x  = sτ * ϕ, I_myopic = i_myopic / (sτ * ϕz * exp(r * (T - τ))), I_hedging = Fz / a / (sτ * ϕz * exp(r * (T - τ))))
end

m = ArbitrageHoldingCosts()
τs = range(m.T, stop = 0, length = 100)
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
y, result, distance = pdesolve(m, stategrid, y0, τs)

### reproduce Fig 2
#d = Normal(0, sqrt(m.σ^2 / (2 * m.ρ)))
#zmin = quantile(d, 0.025)
#zmax = quantile(d, 0.975)
#idx = (state[:z] .>= zmin) .& (state[:z] .<= zmax)
#
#using Plots
#plot(result[:x][idx, 20], [result[:I_myopic][idx, 2] result[:I][idx, 2]], label = ["myopic" "all"])
#plot(result[:x][idx, 20], result[:I][idx, 2])
#plot!(result[:x][idx, 60], result[:I][idx, 60])
#plot!(result[:x][idx, 100], result[:I][idx, 100])
#


text:
using Plots
using Random
Random.seed!(1)

include_ddp()

contact_control_path = "/home/taylor/Research/ContactControl.jl/src"

using Parameters
# Utilities
include(joinpath(contact_control_path, "utils.jl"))

# Solver
include(joinpath(contact_control_path, "solver/cones.jl"))
include(joinpath(contact_control_path, "solver/interior_point.jl"))
include(joinpath(contact_control_path, "solver/lu.jl"))

# Environment
include(joinpath(contact_control_path, "simulator/environment.jl"))

# Dynamics
include(joinpath(contact_control_path, "dynamics/model.jl"))

# Simulator
include(joinpath(contact_control_path, "simulation/contact_methods.jl"))
include(joinpath(contact_control_path, "simulation/simulation.jl"))
include(joinpath(contact_control_path, "simulator/trajectory.jl"))

include(joinpath(contact_control_path, "dynamics/code_gen_dynamics.jl"))
include(joinpath(contact_control_path, "dynamics/fast_methods_dynamics.jl"))

# Models
include(joinpath(contact_control_path, "dynamics/quaternions.jl"))
include(joinpath(contact_control_path, "dynamics/mrp.jl"))
include(joinpath(contact_control_path, "dynamics/euler.jl"))

# include("dynamics/particle_2D/model.jl")
# include("dynamics/particle/model.jl")
include(joinpath(contact_control_path, "dynamics/hopper_2D/model.jl"))
# include("dynamics/hopper_3D/model.jl")
# include("dynamics/hopper_3D_quaternion/model.jl")
# include("dynamics/quadruped/model.jl")
# include("dynamics/quadruped_simple/model.jl")
# include("dynamics/biped/model.jl")
# include("dynamics/flamingo/model.jl")
# include("dynamics/pushbot/model.jl")
# include("dynamics/planarpush/model.jl")
# include("dynamics/planarpush_2D/model.jl")
# include("dynamics/rigidbody/model.jl")
# include("dynamics/box/model.jl")

# Simulation
include(joinpath(contact_control_path, "simulation/environments/flat.jl"))
# include("simulation/environments/piecewise.jl")
# include("simulation/environments/quadratic.jl")
# include("simulation/environments/slope.jl")
# include("simulation/environments/sinusoidal.jl")
# include("simulation/environments/stairs.jl")

include(joinpath(contact_control_path, "simulation/residual_approx.jl"))
include(joinpath(contact_control_path, "simulation/code_gen_simulation.jl"))

# Visuals
using MeshCatMechanisms
include(joinpath(contact_control_path, "dynamics/visuals.jl"))
include(joinpath(contact_control_path, "dynamics/visual_utils.jl"))

s = get_simulation("hopper_2D", "flat_2D_lc", "flat")

nq = s.model.dim.q
m = s.model.dim.u

T = 11
h = 0.1

q0 = [0.0; 0.5; 0.0; 0.5]
q1 = [0.0; 0.5; 0.0; 0.5]
qT = [0.0; 0.5; 0.0; 0.5]
q_ref = [0.0; 0.5; 0.0; 0.5]

struct Dynamics{T}
	s::Simulation
	ip_dyn::InteriorPoint
	ip_jac::InteriorPoint
	h::T
end

function gen_dynamics(s::Simulation, h;
		dyn_opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 0.1,
						diff_sol = true),
		jac_opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 0.1,
						diff_sol = true))

	z = zeros(num_var(s.model, s.env))
	θ = zeros(num_data(s.model))

	ip_dyn = interior_point(z, θ,
		idx_ineq = inequality_indices(s.model, s.env),
		r! = s.res.r!,
		rz! = s.res.rz!,
		rθ! = s.res.rθ!,
		rz = s.rz,
		rθ = s.rθ,
		opts = dyn_opts)

	ip_dyn.opts.diff_sol = false

	ip_jac = interior_point(z, θ,
		idx_ineq = inequality_indices(s.model, s.env),
		r! = s.res.r!,
		rz! = s.res.rz!,
		rθ! = s.res.rθ!,
		rz = s.rz,
		rθ = s.rθ,
		opts = jac_opts)

	ip_jac.opts.diff_sol = true

	Dynamics(s, ip_dyn, ip_jac, h)
end

d = gen_dynamics(s, h,
	dyn_opts = InteriorPointOptions{Float64}(κ_tol = 1.0e-4, κ_init = 0.1),
	jac_opts = InteriorPointOptions{Float64}(κ_tol = 1.0e-4, κ_init = 0.1))

function f!(d::Dynamics, q0, q1, u1, mode = :dynamics)
	s = d.s
	ip = (mode == :dynamics ? d.ip_dyn : d.ip_jac)
	h = d.h

	z_initialize!(ip.z, s.model, s.env, copy(q1))
	θ_initialize!(ip.θ, s.model, copy(q0), copy(q1), copy(u1), zeros(s.model.dim.w), s.model.μ_world, h)

	status = interior_point_solve!(ip)

	!status && (@warn "dynamics failure")
end

function f(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :dynamics)
	return copy(d.ip_dyn.z[1:d.s.model.dim.q])
end

f(d, q0, q1, zeros(m))

function fq0(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.s.model.dim.q, 1:d.s.model.dim.q])
end

fq0(d, q0, q1, zeros(m))

function fq1(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.s.model.dim.q, d.s.model.dim.q .+ (1:d.s.model.dim.q)])
end

fq1(d, q0, q1, zeros(m))

function fx1(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.s.model.dim.q, 1:(2 * d.s.model.dim.q)])
end

fx1(d, q0, q1, zeros(m))

function fu1(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.s.model.dim.q, 2 * d.s.model.dim.q .+ (1:d.s.model.dim.u)])
end

fu1(d, q0, q1, zeros(m))

struct HopperCI{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int
	dynamics::Dynamics
end

model = HopperCI{Midpoint, FixedTime}(2 * s.model.dim.q, s.model.dim.u, 0, d)

function fd(model::HopperCI{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.s.model.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]

	q2 = f(model.dynamics, q0, q1, u)

	return [q1; q2]
end

fd(model, x1, ū[1], w[1], h, 1)

function fdx(model::HopperCI{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.s.model.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]
	dq2dx1 = fx1(model.dynamics, q0, q1, u)

	return [zeros(nq, nq) I; dq2dx1]
end

fdx(model, x1, ū[1], w[1], h, 1)

function fdu(model::HopperCI{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.s.model.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]
	dq2du1 = fu1(model.dynamics, q0, q1, u)
	return [zeros(nq, model.m); dq2du1]
end

fdu(model, x1, ū[1], w[1], h, 1)

n = model.n
m = model.m

# Time
T = 11
h = 0.1

# Initial conditions, controls, disturbances
q0 = [0.0; 0.5; 0.0; 0.5]
q1 = [0.0; 0.5; 0.0; 0.5]
qT = [1.0; 0.5; 0.0; 0.5]
q_ref = [0.0; 0.5; 0.0; 0.5]

x1 = [q1; q1]
xT = [qT; qT]

ū = [[0.0; s.model.g * (s.model.mb + s.model.ml) * 0.5 * h] for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

# Objective
Q = [t < T ? 0.1 * Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]) : Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]) for t = 1:T]
q = [-2.0 * Q[t] * xT for t = 1:T]
R = [Diagonal(1.0e-3 * ones(model.m)) for t = 1:T-1]
r = [zeros(model.m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
ul = [-10.0; -10.0]
uu = [10.0; 10.0]
p = [t < T ? 2 * m : n for t = 1:T]
info_t = Dict(:ul => ul, :uu => uu, :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		c .= [ul - u; u - uu]
	else
		c .= x - cons.con[T].info[:xT]
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T,
	analytical_dynamics_derivatives = true)

# Solve
@time constrained_ddp_solve!(prob,
	max_iter = 1000, max_al_iter = 5,
	ρ_init = 1.0, ρ_scale = 10.0,
	con_tol = 1.0e-3)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

q̄ = state_to_configuration(x̄)

vis = Visualizer()
render(vis)
visualize!(vis, s.model, q̄, Δt = h)

function visualize!(vis, model, q;
		Δt = 0.1, scenario = :vertical)

    r_foot = 0.05
    r_leg = 0.5 * r_foot

	default_background!(vis)

    setobject!(vis["body"], Sphere(Point3f0(0),
        convert(Float32, 0.1)),
        MeshPhongMaterial(color = RGBA(0, 1, 0, 1.0)))

    setobject!(vis["foot"], Sphere(Point3f0(0),
        convert(Float32, r_foot)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

    n_leg = 100
    for i = 1:n_leg
        setobject!(vis["leg$i"], Sphere(Point3f0(0),
            convert(Float32, r_leg)),
            MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
    end

    p_leg = [zeros(3) for i = 1:n_leg]
    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        p_body = [q[t][1], 0.0, q[t][2]]
        p_foot = [kinematics(model, q[t])[1], 0.0, kinematics(model, q[t])[2]]

        q_tmp = Array(copy(q[t]))
        r_range = range(0, stop = q[t][4], length = n_leg)
        for i = 1:n_leg
            q_tmp[4] = r_range[i]
            p_leg[i] = [kinematics(model, q_tmp)[1], 0.0, kinematics(model, q_tmp)[2]]
        end
        q_tmp[4] = q[t][4]
        p_foot = [kinematics(model, q_tmp)[1], 0.0, kinematics(model, q_tmp)[2]]

        z_shift = [0.0; 0.0; r_foot]

        MeshCat.atframe(anim, t) do
            settransform!(vis["body"], Translation(p_body + z_shift))
            settransform!(vis["foot"], Translation(p_foot + z_shift))

            for i = 1:n_leg
                settransform!(vis["leg$i"], Translation(p_leg[i] + z_shift))
            end
        end
    end

	if scenario == :vertical
		settransform!(vis["/Cameras/default"],
			compose(Translation(0.0, 0.5, -1.0),LinearMap(RotZ(-pi / 2.0))))
	end

    MeshCat.setanimation!(vis, anim)
end

# # Model
# include_model("double_integrator")
#
# function f(model::DoubleIntegratorContinuous, x, u, w)
#     [x[2]; (1.0 + w[1]) * u[1]]
# end
#
# function fd(model::DoubleIntegratorContinuous{Midpoint, FixedTime}, x, u, w, h, t)
# 	x + h * f(model, x + 0.5 * h * f(model, x, u, w), u, w)
# end
#
# model = DoubleIntegratorContinuous{Midpoint, FixedTime}(2, 1, 1)
# n = model.n
# m = model.m
#
# # Time
# T = 11
# h = 0.1
#
# # Initial conditions, controls, disturbances
# x1 = [1.0; 0.0]
# x_ref = [[0.0; 0.0] for t = 1:T]
# xT = [0.0; 0.0]
# ū = [1.0 * randn(model.m) for t = 1:T-1]
# u_ref = [zeros(model.m) for t = 1:T-1]
# w = [zeros(model.d) for t = 1:T-1]
#
# # Rollout
# x̄ = rollout(model, x1, ū, w, h, T)
#
# # Objective
# Q = [(t < T ? h : 1.0) * (t < T ?
# 	 Diagonal([1.0; 1.0])
# 		: Diagonal([1.0; 1.0])) for t = 1:T]
# q = [-2.0 * Q[t] * x_ref[t] for t = 1:T]
# R = h * [Diagonal(1.0 * ones(model.m)) for t = 1:T-1]
# r = [zeros(model.m) for t = 1:T-1]
#
# obj = StageCosts([QuadraticCost(Q[t], q[t],
# 	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)
#
# function g(obj::StageCosts, x, u, t)
# 	T = obj.T
#     if t < T
# 		Q = obj.cost[t].Q
# 		q = obj.cost[t].q
# 	    R = obj.cost[t].R
# 		r = obj.cost[t].r
#         return x' * Q * x + q' * x + u' * R * u + r' * u
#     elseif t == T
# 		Q = obj.cost[T].Q
# 		q = obj.cost[T].q
#         return x' * Q * x + q' * x
#     else
#         return 0.0
#     end
# end
#
# # Constraints
# ul = [-5.0]
# uu = [5.0]
# p = [t < T ? 2 * m : n for t = 1:T]
# info_t = Dict(:ul => ul, :uu => uu, :inequality => (1:2 * m))
# info_T = Dict(:xT => xT)
# con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]
#
# function c!(c, cons::StageConstraints, x, u, t)
# 	T = cons.T
# 	p = cons.con[t].p
#
# 	if t < T
# 		ul = cons.con[t].info[:ul]
# 		uu = cons.con[t].info[:uu]
# 		c .= [ul - u; u - uu]
# 	else
# 		c .= x - cons.con[T].info[:xT]
# 	end
# end
#
# prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T)
#
# # Solve
# @time constrained_ddp_solve!(prob,
# 	max_iter = 1000, max_al_iter = 10,
# 	ρ_init = 1.0, ρ_scale = 10.0,
# 	con_tol = 1.0e-5)
#
# x, u = current_trajectory(prob)
# x̄, ū = nominal_trajectory(prob)
#
# # Visualize
# using Plots
# # plot(hcat([[0.0; 0.0] for t = 1:T]...)',
# #     width = 2.0, color = :black, label = "")
# plt = plot(hcat(x...)',
# 	width = 2.0,
# 	color = [:cyan :orange],
# 	label = ["x" "ẋ"],
# 	xlabel = "time step",
# 	ylabel = "state")
#
# savefig(plt,
# 	joinpath("/home/taylor/Research/parameter_optimization_manuscript/figures/di_base_state.png"))
#
# plt = plot(hcat(u..., u[end])',
# 	width = 2.0,
# 	color = :magenta,
# 	linetype = :steppost,
# 	xlabel = "time step",
# 	ylabel = "control",
# 	label = "")
#
# savefig(plt,
# 	joinpath("/home/taylor/Research/parameter_optimization_manuscript/figures/di_base_control.png"))
#
# # Simulate policy
# include(joinpath(@__DIR__, "simulate.jl"))
#
# # Model
# model_sim = DoubleIntegratorContinuous{RK3, FixedTime}(model.n, model.m, model.d)
# x1_sim = copy(x1)
# T_sim = 10 * T
#
# # Time
# tf = h * (T - 1)
# t = range(0, stop = tf, length = T)
# t_sim = range(0, stop = tf, length = T_sim)
# dt_sim = tf / (T_sim - 1)
#
# # Policy
# K = [K for K in prob.p_data.K]
# plot(vcat(K...))
# K = [prob.p_data.K[t] for t = 1:T-1]
# # K, _ = tvlqr(model, x̄, ū, h, Q, R)
# # # K = [-k for k in K]
# # K = [-K[1] for t = 1:T-1]
# # plot(vcat(K...))
#
# # Simulate
# N_sim = 1
# x_sim = []
# u_sim = []
# J_sim = []
# Random.seed!(1)
# for k = 1:N_sim
# 	wi_sim = 0.0 * min(0.1, max(-0.1, 1.0e-1 * randn(1)[1]))
# 	w_sim = [wi_sim for t = 1:T-1]
# 	println("sim: $k - w = $(wi_sim[1])")
#
# 	x_ddp, u_ddp, J_ddp, Jx_ddp, Ju_ddp = simulate_linear_feedback(
# 		model_sim,
# 		K,
# 	    x̄, ū,
# 		x_ref, u_ref,
# 		Q, R,
# 		T_sim, h,
# 		x1_sim,
# 		w_sim,
# 		ul = ul,
# 		uu = uu)
#
# 	push!(x_sim, x_ddp)
# 	push!(u_sim, u_ddp)
# 	push!(J_sim, J_ddp)
# end
#
# # Visualize
# idx = (1:2)
# plt = plot(t, hcat(x̄...)[idx, :]',
# 	width = 2.0, color = :black, label = "",
# 	xlabel = "time (s)", ylabel = "state",
# 	title = "double integrator (J_avg = $(round(mean(J_sim), digits = 3)), N_sim = $N_sim)")
#
# for xs in x_sim
# 	plt = plot!(t_sim, hcat(xs...)[idx, :]',
# 	    width = 1.0, color = :magenta, label = "")
# end
# display(plt)
#
# plt = plot(t, hcat(ū..., ū[end])',
# 	width = 2.0, color = :black, label = "",
# 	xlabel = "time (s)", ylabel = "control",
# 	linetype = :steppost,
# 	title = "double integrator (J_avg = $(round(mean(J_sim), digits = 3)), N_sim = $N_sim)")
#
# for us in u_sim
# 	plt = plot!(t_sim, hcat(us..., us[end])',
# 		width = 1.0, color = :magenta, label = "",
# 		linetype = :steppost)
# end
# display(plt)
# # u_sim
# # plot(vcat(K...))

max_issues_repo_name: thowell/DirectMotionPlanning


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 0552a944-cdbe-11eb-2e63-49f12bee1624
using FourierTools, TestImages, Colors, ImageShow, ImageCore, Statistics

# ╔═╡ 9a73fe7d-d298-4791-8f29-5a962d1242d1
gauss(y, x, σ=1) = 1 / σ / √(2π) * exp(-0.5 * (x^2 + y^2) / σ^2)

# ╔═╡ b8b2b314-4fed-420e-9061-6c2716c134c8
begin
	y = fftpos(512, 512)
	x = y'
end;

# ╔═╡ 97dc4b9d-0b32-4057-8a22-85cb40ae1ebf
begin
	kernel = ifftshift_view(gauss.(y, x, Ref(3)))
	kernel ./= sum(kernel)
end;

# ╔═╡ 7500e8d2-16a7-4ea3-9192-78b4acfa327e
function conv_pad(A, B, value=zero(eltype(A)), pad=10)
	A_1 = value .+ FourierTools.select_region(A .- value, new_size=size(A) .+ pad)
	B_1 = ifftshift_view(FourierTools.select_region(fftshift_view(B), new_size=size(A) .+ pad))
	return FourierTools.select_region(conv(A_1, B_1), new_size=size(A))
end

# ╔═╡ 44370313-62d3-4774-aed6-72a51410503d
img = Float64.(Gray.(testimage("house")));

# ╔═╡ d190772c-885f-405b-b0c9-bd11b6ac992c
img_blurry = Gray.(conv(kernel, img))

# ╔═╡ 0e57683e-7af8-4940-8494-032d6190f2cf
img_blurry[end-100:end-51, 1:50]

# ╔═╡ b8a8adfb-af4a-4ff9-b739-a040f05e2896
img_blurry2 = Gray.(conv_pad(img, kernel, 0.7, 10))

# ╔═╡ 94816b91-0ff3-41ef-9146-b71f150be161
img_blurry2[end-100:end-51, 1:50]

# ╔═╡ Cell order:
# ╠═0552a944-cdbe-11eb-2e63-49f12bee1624
# ╠═9a73fe7d-d298-4791-8f29-5a962d1242d1
# ╠═b8b2b314-4fed-420e-9061-6c2716c134c8
# ╠═97dc4b9d-0b32-4057-8a22-85cb40ae1ebf
# ╠═7500e8d2-16a7-4ea3-9192-78b4acfa327e
# ╠═44370313-62d3-4774-aed6-72a51410503d
# ╠═d190772c-885f-405b-b0c9-bd11b6ac992c
# ╠═0e57683e-7af8-4940-8494-032d6190f2cf
# ╠═b8a8adfb-af4a-4ff9-b739-a040f05e2896
# ╠═94816b91-0ff3-41ef-9146-b71f150be161


text:
# ===== AbstractValGraph ==========

abstract type AbstractValGraph{V<:Integer, V_VALS, E_VALS} <: AbstractGraph{V} end

eltype(::Type{<:AbstractValGraph{V}}) where {V} = V
eltype(g::AbstractValGraph) = eltype(typeof(g))

# TODO examples
# TODO might also implement this for SimpleGraph
"""
	edgevals_type(g::AbstractValGraph)

Return the types of the edgevalues of a graph `g`.
"""
function edgevals_type end

edgevals_type(::Type{<:AbstractValGraph{V, V_VALS, E_VALS}}) where {V, V_VALS, E_VALS} = E_VALS
edgevals_type(g::AbstractValGraph) = edgevals_type(typeof(g))

vertices(g::AbstractValGraph) = OneTo{eltype(g)}(nv(g))
has_vertex(g::AbstractValGraph, v) = v ∈ vertices(g)

# ===== AbstractEdgeValGraph ==========

abstract type AbstractEdgeValGraph{V<:Integer, E_VALS} <: AbstractValGraph{V, Tuple{}, E_VALS} end

OneEdgeValGraph{V, E_VAL} = AbstractEdgeValGraph{V, E_VALS} where
    E_VALS <: Union{Tuple{<:E_VAL}, NamedTuple{S, Tuple{<:E_VAL}} where S}

ZeroEdgeValGraph{V} = AbstractEdgeValGraph{V, E_VALS} where
    E_VALS <: Union{Tuple{}, NamedTuple{S, Tuple{}} where S}

const default_edgevals = (Float64,)

const default_edgevals_type = Tuple{Float64}

#= TODO do we still need these?
default_edgeval(E_VAL)::E_VAL = oneunit(E_VAL)
default_edgeval(::Type{Nothing})::Nothing = nothing
default_edgeval(T::Type{<:AbstractTuple})::T = T( default_edgeval(U) for U in T.types )

default_zero_edgeval(::Type{Nothing}) = nothing
default_zero_edgeval(E_VAL) = zero(E_VAL)
default_zero_edgeval(E_VAL::Type{<:AbstractTuple}) = E_VAL( default_zero_edgeval(T) for T in E_VAL.types )
=#

edgetype(G::Type{<:AbstractEdgeValGraph{V, E_VALS}}) where {V, E_VALS} =
    is_directed(G) ? ValDiEdge{V, E_VALS} : ValEdge{V, E_VALS}

edgetype(g::AbstractEdgeValGraph) = edgetype(typeof(g))

edges(g::AbstractEdgeValGraph) = ValEdgeIter(g)
edges_for_key(g::AbstractEdgeValGraph, key) = ValEdgeIter(g, g.edgevals[key])

nv(g::AbstractEdgeValGraph) = eltype(g)(length(g.fadjlist))
ne(g::AbstractEdgeValGraph) = g.ne

# TODO do we need this?
zero(G::AbstractEdgeValGraph{V}) where {V} = G(zero(V))

function is_validkey(
            g::AbstractEdgeValGraph{V, E_VALS},
            key::Symbol) where {V, E_VALS <: NamedTuple}

    return key in E_VALS.names
end

function is_validkey(
            g::AbstractEdgeValGraph{V, E_VALS},
            key::Integer) where {V, E_VALS <: AbstractTuple}

    return key in OneTo(length(E_VALS.types))
end

function validkey_or_throw(g, key)
    is_validkey(g, key) && return nothing

    error("$key is not a valid edge key for this graph.")
end

# === Iterators =====================

struct ValEdgeIter{G<:AbstractEdgeValGraph} <: AbstractEdgeIter
    g::G
end

length(iter::ValEdgeIter) = ne(iter.g)

function eltype(::Type{<:ValEdgeIter{G}}) where
        {V, E_VALS, G <: AbstractEdgeValGraph{V, E_VALS}}

    return (is_directed(G) ? ValDiEdge : ValEdge){V, E_VALS}
end
eltype(iter::ValEdgeIter) = eltype(typeof(iter))

#
# === Display =====================

# TODO output should be actuall graph type
function show(io::IO, g::AbstractEdgeValGraph)
    dir = is_directed(g) ? "directed" : "undirected"
    println(io, "{$(nv(g)), $(ne(g))} $dir $(eltype(g)) ValueGraph with edge values of type $(edgevals_type(g))")
end


text:
using Plots
using jInv.Mesh
using Flux
using Printf
using Statistics
using LinearAlgebra
using JLD
using MFGnet
include("viewers.jl")
include("runOMThelpers.jl")

# @isdefined(d)  ? d = d : d = 2
# @isdefined(level)  ? level = level : level = 2
# @isdefined(iter)  ? iter = iter : iter = 500
for d=[2]
    for level=[1,2]
        for iter=[100 200 300 400 500]

dir = pwd()

# fname = "OMT-Multilevel-nt-2-m-16-d-2-level-2"
fname = "OMT-Multilevel-noHJB-nt-8-d-$d-level-$level-iter-$iter"
imgdir = dir * "/" * fname *"/"
println("fname=$fname")
mkdir(imgdir)
file = fname * ".jld"
res = load(file)
settings = res["settings"]
(m,α,nTrain,R,d,domain,nTh,m,nt,T) = settings

stepper = RK4Step()
if !@isdefined(ar)
    ar = x-> R.(x)
end
Θbest = res["Θbest"]
println("number of model parameters: $(length(MFGnet.param2vec(Θbest)))")

# rho1 = Gaussian(d)
Gs = Array{Gaussian{R,Vector{R}}}(undef,8)
ang = range(0,stop=2*pi, length=length(Gs)+1)[1:end-1]
sig = R.(.3*ones(d))
for k=1:length(ang)
    μk = R.(4*([cos(ang[k]); sin(ang[k]);zeros(d-2)]))
    Gs[k] = Gaussian(d,sig,μk,R(1.0/length(Gs)))
end
rho0 = GaussianMixture(Gs)
rho1 = Gaussian(d,sig,zeros(d),1.0)

M   = getRegularMesh(domain, [128, 128])
X0  = Matrix(getCellCenteredGrid(M)')
X0  = ar([X0; zeros(d-2,size(X0,2))])

Φ = getPotentialResNet(nTh, T, nTh, R)

# setup validation
rho0x = rho0(X0)
rho1x = rho1(X0)

vol     = (domain[2]-domain[1])^d
w         = (vol/M.nc) * rho0x
F = F0()
G = Gkl(rho0,rho1,rho0x,rho1x,R.(1.0))
J = MeanFieldGame(F,G,X0,rho0,w,Φ=Φ,stepper=stepper,nt=nt,α=α)
J.rho0x = rho0(J.X0)
J.G.rho0x = rho0(J.X0)
J.G.rho1x = rho1(J.X0)

Jc  = J(Θbest)

## plots:
p1= viewImage2D(J.rho0x,M,aspect_ratio=:equal)
title!("rho0(x)")
savefig(p1,imgdir * "rho0x-n-$(M.n[1])x$(M.n[2]).png")

p2= viewImage2D(J.G.rho1x,M,aspect_ratio=:equal)
title!("rho1(x)")
savefig(p2,imgdir * "rho1x-n-$(M.n[1])x$(M.n[2]).png")

detDy = exp.(J.UN[d+1,:])
rho1y = (J.G.rho1(J.UN[1:d,:]).*detDy)
p3= viewImage2D(rho1y ,M,aspect_ratio=:equal,clims=(minimum(J.rho0x),maximum(J.rho0x)))
# p3= viewImage2D(rho1y ,M,aspect_ratio=:equal)
title!("rho1(y).*det")
savefig(p3,imgdir * "rho1y-n-$(M.n[1])x$(M.n[2]).png")

res0 = J.rho0x - rho1y
p4 = viewImage2D(abs.(res0),M,aspect_ratio=:equal)
title!("| rho0(x)-rho1(y).*det |") # where's the determinant
savefig(p4,imgdir * "absDiffRho-n-$(M.n[1])x$(M.n[2]).png")

p5 = viewImage2D(detDy,M,aspect_ratio=:equal)
title!("det")
savefig(p5,imgdir * "det-n-$(M.n[1])x$(M.n[2]).png")

xc = [ J.X0; fill(0.0,1,size(J.X0,2))]
Φ0 = J.Φ(xc,Θbest)
p6 = viewImage2D(Φ0,M,aspect_ratio=:equal)
title!("Phi(x,0)")
savefig(p6,imgdir * "phi0-n-$(M.n[1])x$(M.n[2]).png")

xc = [ J.UN[1:d,:]; fill(T,1,size(J.X0,2))]
Φ1z = vec(J.Φ(xc,Θbest) )
p7 = viewImage2D(Φ1z,M,aspect_ratio=:equal)
title!("Phi(z,1)")
savefig(p7,imgdir * "phi1z-n-$(M.n[1])x$(M.n[2]).png")

deltaG = J.α[3]*MFGnet.getDeltaG(J.G,J.UN)
p8 = viewImage2D(deltaG,M,aspect_ratio=:equal)
title!("deltaG(z)")
savefig(p8,imgdir * "deltaG-n-$(M.n[1])x$(M.n[2]).png")

resHJB = abs.(Φ1z-deltaG)
p9 = viewImage2D(resHJB,M,aspect_ratio=:equal)
title!("|Phi1 - deltaG(z)|")
savefig(p9,imgdir * "absDiffHJBfinal-n-$(M.n[1])x$(M.n[2]).png")

# ###### Create new MFG for characteristics
Xc     = load("OT-X0.jld")["Xt"]
if d>2
    Xc = [Xc; zeros(d-2,size(Xc,2))]
    # [3:d,:] .= R(0.0)
end
wc     = 1/size(Xc,2) * ones(size(Xc,2))
rho0xc = rho0(Xc)
rho1xc = rho1(Xc)
Fc = F0()
Gc = Gkl(rho0,rho1,rho0xc,rho1xc,1.0)
Jc = MeanFieldGame(Fc,Gc,Xc,rho0,wc,Φ=Φ,stepper=RK4Step(),nt=4*J.nt,α=J.α,tspan=J.tspan)

p10 = viewImage2D(J.rho0x ,M,aspect_ratio=:equal)
Ut = MFGnet.integrate2(Jc.stepper,MFGnet.odefun,Jc,ar([Xc;zeros(4,size(Xc,2))]),Θbest,[0.0 T],Jc.nt)
charFwd = copy(Ut);
Y0  = [Ut[1:d,:,end];zeros(4,size(Xc,2))]
U0 = MFGnet.integrate2(Jc.stepper,MFGnet.odefun,Jc,Y0,Θbest,[T 0.0],Jc.nt)
charBwd = copy(U0);
for j=1:size(Xc,2)
    uj = Ut[1:d,j,:]
    vj = U0[1:d,j,:]
    plot!(p10,uj[1,:],uj[2,:],legend=false,linecolor=:white,linewidth=1.5)
    plot!(p10,vj[1,:],vj[2,:],legend=false,linecolor=:red,linewidth=1.5)
end
title!("characteristics, fwd+inv")
savefig(p10,imgdir * "characteristics-n-$(M.n[1])x$(M.n[2]).png")

U0 = MFGnet.integrate(J.stepper,MFGnet.odefun,J,ar([J.X0;zeros(4,size(J.X0,2))]),Θbest,[T 0.0],J.nt)
detDyInv = exp.(U0[d+1,:])
rho0z = rho0(U0[1:d,:]).*detDyInv
p11= viewImage2D(rho0z ,M,aspect_ratio=:equal)
title!("rho0(z).*detInv")
savefig(p11,imgdir * "rho0z-n-$(M.n[1])x$(M.n[2]).png")

res1 = rho1x - rho0z
p4 = viewImage2D(abs.(res1),M,aspect_ratio=:equal)
title!("| rho1(x)-rho0(z).*detInv |") # where's the determinant
savefig(p4,imgdir * "absDiffRho1-n-$(M.n[1])x$(M.n[2]).png")

xc = [ J.X0[1:d,:]; fill(T,1,size(J.X0,2))]
Φ1x = vec(J.Φ(xc,Θbest) )
p12 = viewImage2D(Φ1x,M,aspect_ratio=:equal)
title!("Phi(x,1)")
savefig(p12,imgdir * "phi1-n-$(M.n[1])x$(M.n[2]).png")

X1,X2 = getFaceGrids(M)
tt = range(0,stop=T,length=64)
X1t = t -> [Matrix(X1'); zeros(d-2,size(X1,1)); fill(t,1,size(X1,1))]
X2t = t -> [Matrix(X2'); zeros(d-2,size(X2,1)); fill(t,1,size(X2,1))]
XC = getCellCenteredGrid(M)
XCt = t -> [Matrix(XC'); zeros(d-2,size(XC,1)); fill(t,1,size(XC,1))]
V1 = zeros(M.n[1]+1,M.n[2],length(tt))
V2 = zeros(M.n[1],M.n[2]+1,length(tt))
ΦOpt  = zeros(M.n[1],M.n[2],length(tt))

for k=1:length(tt)
    tk = tt[k]
    println("tk=$tk")
    Phi = J.Φ(X1t(tk),Θbest) # run fwd prop to populate N.tmp
    gradPhi = MFGnet.getGradPotential(J.Φ,X1t(tk),Θbest)
    V1[:,:,k] = reshape(-(1/J.α[1])*gradPhi[1,:],M.n[1]+1,M.n[2])
    Phi = J.Φ(X2t(tk),Θbest) # run fwd prop to populate N.tmp
    gradPhi = MFGnet.getGradPotential(J.Φ,X2t(tk),Θbest)
    V2[:,:,k] = reshape(-(1/J.α[1])*gradPhi[2,:],M.n[1],M.n[2]+1)
    ΦOpt[:,:,k] = reshape(J.Φ(XCt(tk),Θbest),M.n[1],M.n[2])
end
using MAT
matwrite( fname * ".mat", Dict(
    "rho0x" => rho0x,
    "rho1x" => rho1x,
    "rho1y" => rho1y,
    "rho0z" => rho0z,
    "domain" => M.domain,
    "detDy" => detDy,
    "detDyInv" => detDyInv,
    "res0" => res0,
    "res1" => res1,
    "charFwd" => charFwd,
    "charBwd" => charBwd,
    "deltaG" => deltaG,
    "n" => M.n,
    "Phi0" => Φ0,
    "Phi1z" => Φ1z,
    "Phi1" => Φ1x,
    "V1" => V1,
    "V2" => V2,
    "PhiOpt" => ΦOpt,
    "His" => res["His"]))
end

end

end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using MCMCBenchmarks,Distributed
Nchains = 4
setprocs(Nchains)

ProjDir = @__DIR__
cd(ProjDir)

isdir("tmp") && rm("tmp", recursive=true)
mkdir("tmp")
!isdir("results") && mkdir("results")
path = pathof(MCMCBenchmarks)
@everywhere begin
  using MCMCBenchmarks
  #Model and configuration patterns for each sampler are located in a
  #seperate model file.
  include(joinpath($path, "../../Models/SDT/SDT.jl"))
  include(joinpath($path, "../../Models/SDT/SDT_Functions.jl"))
end
include(joinpath(path, "../../Models/SDT/SDT.jl"))
include(joinpath(path, "../../Models/SDT/SDT_Functions.jl"))
#Model and configuration patterns for each sampler are located in a
#seperate model file.

@everywhere Turing.turnprogress(false)
#set seeds on each processor
seeds = (939388,39884,28484,495858,544443)
for (i,seed) in enumerate(seeds)
    @fetch @spawnat i Random.seed!(seed)
end

stanSampler = CmdStanNUTS(CmdStanConfig,ProjDir)
#Initialize model files for each instance of stan
initStan(stanSampler)
#Compile stan model

samplers=(CmdStanNUTS(CmdStanConfig,ProjDir),
    AHMCNUTS(AHMC_SDT,AHMCconfig),
    DHMCNUTS(sampleDHMC)
)

#Number of data points
Nd = [10, 100, 1000, 10_000]

#Number of simulations
Nreps = 100

options = (Nsamples=2000, Nadapt=1000, delta=.8, Nd=Nd)
#perform the benchmark
results = pbenchmark(samplers, simulateSDT, Nreps; options...)

#save results
save(results,ProjDir)

# Make plots
include("primary_plots.jl")