text:
const pentatope1an = [
    SVector{5}([0.2 0.2 0.2 0.2 0.2] ),
]

const pentatope1wn = [
    1.0
]

# x = [0.11835034190712865]
# w = [0.2]

const pentatope5an = [
  SVector{5}([0.118350341907226 0.118350341907226 0.118350341907226 0.118350341907226 0.526598632371096] ),
  SVector{5}([0.118350341907226 0.118350341907226 0.118350341907226 0.526598632371096 0.118350341907226] ),
  SVector{5}([0.118350341907226 0.118350341907226 0.526598632371096 0.118350341907226 0.118350341907226] ),
  SVector{5}([0.118350341907226 0.526598632371096 0.118350341907226 0.118350341907226 0.118350341907226] ),
  SVector{5}([0.526598632371096 0.118350341907226 0.118350341907226 0.118350341907226 0.118350341907226] )
]

const pentatope5wn = [
  0.2
  0.2
  0.2
  0.2
  0.2
]

# x = [0.057076947867029514, 0.08285226049309079]
# w = [0.019898619367251447, 0.09005069031637428]

const pentatope15an = [
  SVector{5}([0.0570769479251544 0.0570769479251544 0.0570769479251544 0.0570769479251544 0.771692208299383]),
  SVector{5}([0.0570769479251544 0.0570769479251544 0.0570769479251544 0.771692208299383 0.0570769479251544]),
  SVector{5}([0.0570769479251544 0.0570769479251544 0.771692208299383 0.0570769479251544 0.0570769479251544]),
  SVector{5}([0.0570769479251544 0.771692208299383 0.0570769479251544 0.0570769479251544 0.0570769479251544]),
  SVector{5}([0.771692208299383 0.0570769479251544 0.0570769479251544 0.0570769479251544 0.0570769479251544]),
  SVector{5}([0.0828522604971176 0.0828522604971176 0.0828522604971176 0.375721609254324 0.375721609254324] ),
  SVector{5}([0.0828522604971176 0.0828522604971176 0.375721609254324 0.0828522604971176 0.375721609254324] ),
  SVector{5}([0.0828522604971176 0.0828522604971176 0.375721609254324 0.375721609254324 0.0828522604971176] ),
  SVector{5}([0.0828522604971176 0.375721609254324 0.0828522604971176 0.0828522604971176 0.375721609254324] ),
  SVector{5}([0.0828522604971176 0.375721609254324 0.0828522604971176 0.375721609254324 0.0828522604971176] ),
  SVector{5}([0.0828522604971176 0.375721609254324 0.375721609254324 0.0828522604971176 0.0828522604971176] ),
  SVector{5}([0.375721609254324 0.0828522604971176 0.0828522604971176 0.0828522604971176 0.375721609254324] ),
  SVector{5}([0.375721609254324 0.0828522604971176 0.0828522604971176 0.375721609254324 0.0828522604971176] ),
  SVector{5}([0.375721609254324 0.0828522604971176 0.375721609254324 0.0828522604971176 0.0828522604971176] ),
  SVector{5}([0.375721609254324 0.375721609254324 0.0828522604971176 0.0828522604971176 0.0828522604971176] )
]

const pentatope15wn = [
  0.0198986193930243,
  0.0198986193930243,
  0.0198986193930243,
  0.0198986193930243,
  0.0198986193930243,
  0.09005069030348785,
  0.09005069030348785,
  0.09005069030348785,
  0.09005069030348785,
  0.09005069030348785,
  0.09005069030348785,
  0.09005069030348785,
  0.09005069030348785,
  0.09005069030348785,
  0.09005069030348785,
]

#x = [0.08639272923224803, 0.02401496720069173, 0.29381800402891267, 0.06247517556325446]
#w = [0.05144687284128575, 0.010758106723214462, 0.031759228428071336]
const pentatope35an = [

    SVector{5}([0.086392729232226 0.086392729232226 0.086392729232226 0.086392729232226 0.654429083071094]),
    SVector{5}([0.086392729232226 0.086392729232226 0.086392729232226 0.654429083071094 0.086392729232226]),
    SVector{5}([0.086392729232226 0.086392729232226 0.654429083071094 0.086392729232226 0.086392729232226]),
    SVector{5}([0.086392729232226 0.654429083071094 0.086392729232226 0.086392729232226 0.086392729232226]),
    SVector{5}([0.654429083071094 0.086392729232226 0.086392729232226 0.086392729232226 0.086392729232226]),
    SVector{5}([0.024014967200824 0.024014967200824 0.024014967200824 0.463977549198764 0.463977549198764]),
    SVector{5}([0.024014967200824 0.024014967200824 0.463977549198764 0.024014967200824 0.463977549198764]),
    SVector{5}([0.024014967200824 0.024014967200824 0.463977549198764 0.463977549198764 0.024014967200824]),
    SVector{5}([0.024014967200824 0.463977549198764 0.024014967200824 0.024014967200824 0.463977549198764]),
    SVector{5}([0.024014967200824 0.463977549198764 0.024014967200824 0.463977549198764 0.024014967200824]),
    SVector{5}([0.024014967200824 0.463977549198764 0.463977549198764 0.024014967200824 0.024014967200824]),
    SVector{5}([0.463977549198764 0.024014967200824 0.024014967200824 0.024014967200824 0.463977549198764]),
    SVector{5}([0.463977549198764 0.024014967200824 0.024014967200824 0.463977549198764 0.024014967200824]),
    SVector{5}([0.463977549198764 0.024014967200824 0.463977549198764 0.024014967200824 0.024014967200824]),
    SVector{5}([0.463977549198764 0.463977549198764 0.024014967200824 0.024014967200824 0.024014967200824]),
    SVector{5}([0.293818004028897 0.293818004028897 0.293818004028897 0.056070812349935 0.062475175563373]),
    SVector{5}([0.293818004028897 0.293818004028897 0.293818004028897 0.062475175563373 0.056070812349935]),
    SVector{5}([0.293818004028897 0.293818004028897 0.056070812349935 0.293818004028897 0.062475175563373]),
    SVector{5}([0.293818004028897 0.293818004028897 0.056070812349935 0.062475175563373 0.293818004028897]),
    SVector{5}([0.293818004028897 0.293818004028897 0.062475175563373 0.293818004028897 0.056070812349935]),
    SVector{5}([0.293818004028897 0.293818004028897 0.062475175563373 0.056070812349935 0.293818004028897]),
    SVector{5}([0.293818004028897 0.056070812349935 0.293818004028897 0.293818004028897 0.062475175563373]),
    SVector{5}([0.293818004028897 0.056070812349935 0.293818004028897 0.062475175563373 0.293818004028897]),
    SVector{5}([0.293818004028897 0.056070812349935 0.062475175563373 0.293818004028897 0.293818004028897]),
    SVector{5}([0.293818004028897 0.062475175563373 0.293818004028897 0.293818004028897 0.056070812349935]),
    SVector{5}([0.293818004028897 0.062475175563373 0.293818004028897 0.056070812349935 0.293818004028897]),
    SVector{5}([0.293818004028897 0.062475175563373 0.056070812349935 0.293818004028897 0.293818004028897]),
    SVector{5}([0.056070812349935 0.293818004028897 0.293818004028897 0.293818004028897 0.062475175563373]),
    SVector{5}([0.056070812349935 0.293818004028897 0.293818004028897 0.062475175563373 0.293818004028897]),
    SVector{5}([0.056070812349935 0.293818004028897 0.062475175563373 0.293818004028897 0.293818004028897]),
    SVector{5}([0.056070812349935 0.062475175563373 0.293818004028897 0.293818004028897 0.293818004028897]),
    SVector{5}([0.062475175563373 0.293818004028897 0.293818004028897 0.293818004028897 0.056070812349935]),
    SVector{5}([0.062475175563373 0.293818004028897 0.293818004028897 0.056070812349935 0.293818004028897]),
    SVector{5}([0.062475175563373 0.293818004028897 0.056070812349935 0.293818004028897 0.293818004028897]),
    SVector{5}([0.062475175563373 0.056070812349935 0.293818004028897 0.293818004028897 0.293818004028897])
]

const pentatope35wn = [
    0.05144687284124748,
    0.05144687284124748,
    0.05144687284124748,
    0.05144687284124748,
    0.05144687284124748,
    0.010758106723251771,
    0.010758106723251771,
    0.010758106723251771,
    0.010758106723251771,
    0.010758106723251771,
    0.010758106723251771,
    0.010758106723251771,
    0.010758106723251771,
    0.010758106723251771,
    0.010758106723251771,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246,
    0.031759228428062246
]

#x = [0.10992516562792191, 0.2415199352789164, 0.07254405441102943, 0.02167246479934806, 0.160496134843098, 0.005100356073981291, 0.22307374869322513]
#w = [ 0.03901739660684611, 0.055608078377922185, 0.030645207859279433, 0.0025320024042530293, 0.0056593499466101186]
const pentatope70an = [
  SVector{5}([0.109925165627922 0.109925165627922 0.109925165627922 0.109925165627922 0.560299337488312]    ),
  SVector{5}([0.109925165627922 0.109925165627922 0.109925165627922 0.560299337488312 0.109925165627922]    ),
  SVector{5}([0.109925165627922 0.109925165627922 0.560299337488312 0.109925165627922 0.109925165627922]    ),
  SVector{5}([0.109925165627922 0.560299337488312 0.109925165627922 0.109925165627922 0.109925165627922]    ),
  SVector{5}([0.560299337488312 0.109925165627922 0.109925165627922 0.109925165627922 0.109925165627922]    ),
  SVector{5}([0.241519935278916 0.241519935278916 0.241519935278916 0.241519935278916 0.0339202588843344]   ),
  SVector{5}([0.241519935278916 0.241519935278916 0.241519935278916 0.0339202588843344 0.241519935278916]   ),
  SVector{5}([0.241519935278916 0.241519935278916 0.0339202588843344 0.241519935278916 0.241519935278916]   ),
  SVector{5}([0.241519935278916 0.0339202588843344 0.241519935278916 0.241519935278916 0.241519935278916]   ),
  SVector{5}([0.0339202588843344 0.241519935278916 0.241519935278916 0.241519935278916 0.241519935278916]   ),
  SVector{5}([0.0725440544110294 0.0725440544110294 0.0725440544110294 0.391183918383456 0.391183918383456] ),
  SVector{5}([0.0725440544110294 0.0725440544110294 0.391183918383456 0.0725440544110294 0.391183918383456] ),
  SVector{5}([0.0725440544110294 0.0725440544110294 0.391183918383456 0.391183918383456 0.0725440544110294] ),
  SVector{5}([0.0725440544110294 0.391183918383456 0.0725440544110294 0.0725440544110294 0.391183918383456] ),
  SVector{5}([0.0725440544110294 0.391183918383456 0.0725440544110294 0.391183918383456 0.0725440544110294] ),
  SVector{5}([0.0725440544110294 0.391183918383456 0.391183918383456 0.0725440544110294 0.0725440544110294] ),
  SVector{5}([0.391183918383456 0.0725440544110294 0.0725440544110294 0.0725440544110294 0.391183918383456] ),
  SVector{5}([0.391183918383456 0.0725440544110294 0.0725440544110294 0.391183918383456 0.0725440544110294] ),
  SVector{5}([0.391183918383456 0.0725440544110294 0.391183918383456 0.0725440544110294 0.0725440544110294] ),
  SVector{5}([0.391183918383456 0.391183918383456 0.0725440544110294 0.0725440544110294 0.0725440544110294] ),
  SVector{5}([0.0216724647993481 0.0216724647993481 0.0216724647993481 0.160496134843098 0.774486470758858] ),
  SVector{5}([0.0216724647993481 0.0216724647993481 0.0216724647993481 0.774486470758858 0.160496134843098] ),
  SVector{5}([0.0216724647993481 0.0216724647993481 0.160496134843098 0.0216724647993481 0.774486470758858] ),
  SVector{5}([0.0216724647993481 0.0216724647993481 0.160496134843098 0.774486470758858 0.0216724647993481] ),
  SVector{5}([0.0216724647993481 0.0216724647993481 0.774486470758858 0.0216724647993481 0.160496134843098] ),
  SVector{5}([0.0216724647993481 0.0216724647993481 0.774486470758858 0.160496134843098 0.0216724647993481] ),
  SVector{5}([0.0216724647993481 0.160496134843098 0.0216724647993481 0.0216724647993481 0.774486470758858] ),
  SVector{5}([0.0216724647993481 0.160496134843098 0.0216724647993481 0.774486470758858 0.0216724647993481] ),
  SVector{5}([0.0216724647993481 0.160496134843098 0.774486470758858 0.0216724647993481 0.0216724647993481] ),
  SVector{5}([0.0216724647993481 0.774486470758858 0.0216724647993481 0.0216724647993481 0.160496134843098] ),
  SVector{5}([0.0216724647993481 0.774486470758858 0.0216724647993481 0.160496134843098 0.0216724647993481] ),
  SVector{5}([0.0216724647993481 0.774486470758858 0.160496134843098 0.0216724647993481 0.0216724647993481] ),
  SVector{5}([0.160496134843098 0.0216724647993481 0.0216724647993481 0.0216724647993481 0.774486470758858] ),
  SVector{5}([0.160496134843098 0.0216724647993481 0.0216724647993481 0.774486470758858 0.0216724647993481] ),
  SVector{5}([0.160496134843098 0.0216724647993481 0.774486470758858 0.0216724647993481 0.0216724647993481] ),
  SVector{5}([0.160496134843098 0.774486470758858 0.0216724647993481 0.0216724647993481 0.0216724647993481] ),
  SVector{5}([0.774486470758858 0.0216724647993481 0.0216724647993481 0.0216724647993481 0.160496134843098] ),
  SVector{5}([0.774486470758858 0.0216724647993481 0.0216724647993481 0.160496134843098 0.0216724647993481] ),
  SVector{5}([0.774486470758858 0.0216724647993481 0.160496134843098 0.0216724647993481 0.0216724647993481] ),
  SVector{5}([0.774486470758858 0.160496134843098 0.0216724647993481 0.0216724647993481 0.0216724647993481] ),
  SVector{5}([0.00510035607398129 0.00510035607398129 0.223073748693225 0.223073748693225 0.543651790465587]),
  SVector{5}([0.00510035607398129 0.00510035607398129 0.223073748693225 0.543651790465587 0.223073748693225]),
  SVector{5}([0.00510035607398129 0.00510035607398129 0.543651790465587 0.223073748693225 0.223073748693225]),
  SVector{5}([0.00510035607398129 0.223073748693225 0.00510035607398129 0.223073748693225 0.543651790465587]),
  SVector{5}([0.00510035607398129 0.223073748693225 0.00510035607398129 0.543651790465587 0.223073748693225]),
  SVector{5}([0.00510035607398129 0.223073748693225 0.223073748693225 0.00510035607398129 0.543651790465587]),
  SVector{5}([0.00510035607398129 0.223073748693225 0.223073748693225 0.543651790465587 0.00510035607398129]),
  SVector{5}([0.00510035607398129 0.223073748693225 0.543651790465587 0.00510035607398129 0.223073748693225]),
  SVector{5}([0.00510035607398129 0.223073748693225 0.543651790465587 0.223073748693225 0.00510035607398129]),
  SVector{5}([0.00510035607398129 0.543651790465587 0.00510035607398129 0.223073748693225 0.223073748693225]),
  SVector{5}([0.00510035607398129 0.543651790465587 0.223073748693225 0.00510035607398129 0.223073748693225]),
  SVector{5}([0.00510035607398129 0.543651790465587 0.223073748693225 0.223073748693225 0.00510035607398129]),
  SVector{5}([0.223073748693225 0.00510035607398129 0.00510035607398129 0.223073748693225 0.543651790465587]),
  SVector{5}([0.223073748693225 0.00510035607398129 0.00510035607398129 0.543651790465587 0.223073748693225]),
  SVector{5}([0.223073748693225 0.00510035607398129 0.223073748693225 0.00510035607398129 0.543651790465587]),
  SVector{5}([0.223073748693225 0.00510035607398129 0.223073748693225 0.543651790465587 0.00510035607398129]),
  SVector{5}([0.223073748693225 0.00510035607398129 0.543651790465587 0.00510035607398129 0.223073748693225]),
  SVector{5}([0.223073748693225 0.00510035607398129 0.543651790465587 0.223073748693225 0.00510035607398129]),
  SVector{5}([0.223073748693225 0.223073748693225 0.00510035607398129 0.00510035607398129 0.543651790465587]),
  SVector{5}([0.223073748693225 0.223073748693225 0.00510035607398129 0.543651790465587 0.00510035607398129]),
  SVector{5}([0.223073748693225 0.223073748693225 0.543651790465587 0.00510035607398129 0.00510035607398129]),
  SVector{5}([0.223073748693225 0.543651790465587 0.00510035607398129 0.00510035607398129 0.223073748693225]),
  SVector{5}([0.223073748693225 0.543651790465587 0.00510035607398129 0.223073748693225 0.00510035607398129]),
  SVector{5}([0.223073748693225 0.543651790465587 0.223073748693225 0.00510035607398129 0.00510035607398129]),
  SVector{5}([0.543651790465587 0.00510035607398129 0.00510035607398129 0.223073748693225 0.223073748693225]),
  SVector{5}([0.543651790465587 0.00510035607398129 0.223073748693225 0.00510035607398129 0.223073748693225]),
  SVector{5}([0.543651790465587 0.00510035607398129 0.223073748693225 0.223073748693225 0.00510035607398129]),
  SVector{5}([0.543651790465587 0.223073748693225 0.00510035607398129 0.00510035607398129 0.223073748693225]),
  SVector{5}([0.543651790465587 0.223073748693225 0.00510035607398129 0.223073748693225 0.00510035607398129]),
  SVector{5}([0.543651790465587 0.223073748693225 0.223073748693225 0.00510035607398129 0.00510035607398129])
]

const pentatope70wn = [
  0.03901739660684611,
  0.03901739660684611,
  0.03901739660684611,
  0.03901739660684611,
  0.03901739660684611,
  0.055608078377922185,
  0.055608078377922185,
  0.055608078377922185,
  0.055608078377922185,
  0.055608078377922185,
  0.030645207859279433,
  0.030645207859279433,
  0.030645207859279433,
  0.030645207859279433,
  0.030645207859279433,
  0.030645207859279433,
  0.030645207859279433,
  0.030645207859279433,
  0.030645207859279433,
  0.030645207859279433,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0025320024042530293,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186,
  0.0056593499466101186
]

#x = [0.031002102547268737, 0.02871807987242527, 0.3047795257048887, 0.1632598842559579, 0.469494323021712, 0.19097656767256832, 0.037148751846917294, 0.02704627725162003, 0.1392370880634157, 0.036647909737317054, 0.3809416954893026]
#w = [0.04654868655484498, 0.0014672641241909648, 0.0029906128790727765, 0.009907301984077122, 0.011320943419643296, 0.0046932055563252825, 0.010698055349285926]
const pentatope126an = [
  SVector{5}([0.2000000000000000 0.2000000000000000 0.2000000000000000 0.2000000000000000 0.2000000000000000]),
  SVector{5}([0.0310021025472687 0.0310021025472687 0.0310021025472687 0.0310021025472687 0.875991589810925] ),
  SVector{5}([0.0310021025472687 0.0310021025472687 0.0310021025472687 0.875991589810925 0.0310021025472687] ),
  SVector{5}([0.0310021025472687 0.0310021025472687 0.875991589810925 0.0310021025472687 0.0310021025472687] ),
  SVector{5}([0.0310021025472687 0.875991589810925 0.0310021025472687 0.0310021025472687 0.0310021025472687] ),
  SVector{5}([0.875991589810925 0.0310021025472687 0.0310021025472687 0.0310021025472687 0.0310021025472687] ),
  SVector{5}([0.0287180798724253 0.0287180798724253 0.0287180798724253 0.304779525704889 0.609066234677835]  ),
  SVector{5}([0.0287180798724253 0.0287180798724253 0.0287180798724253 0.609066234677835 0.304779525704889]  ),
  SVector{5}([0.0287180798724253 0.0287180798724253 0.304779525704889 0.0287180798724253 0.609066234677835]  ),
  SVector{5}([0.0287180798724253 0.0287180798724253 0.304779525704889 0.609066234677835 0.0287180798724253]  ),
  SVector{5}([0.0287180798724253 0.0287180798724253 0.609066234677835 0.0287180798724253 0.304779525704889]  ),
  SVector{5}([0.0287180798724253 0.0287180798724253 0.609066234677835 0.304779525704889 0.0287180798724253]  ),
  SVector{5}([0.0287180798724253 0.304779525704889 0.0287180798724253 0.0287180798724253 0.609066234677835]  ),
  SVector{5}([0.0287180798724253 0.304779525704889 0.0287180798724253 0.609066234677835 0.0287180798724253]  ),
  SVector{5}([0.0287180798724253 0.304779525704889 0.609066234677835 0.0287180798724253 0.0287180798724253]  ),
  SVector{5}([0.0287180798724253 0.609066234677835 0.0287180798724253 0.0287180798724253 0.304779525704889]  ),
  SVector{5}([0.0287180798724253 0.609066234677835 0.0287180798724253 0.304779525704889 0.0287180798724253]  ),
  SVector{5}([0.0287180798724253 0.609066234677835 0.304779525704889 0.0287180798724253 0.0287180798724253]  ),
  SVector{5}([0.304779525704889 0.0287180798724253 0.0287180798724253 0.0287180798724253 0.609066234677835]  ),
  SVector{5}([0.304779525704889 0.0287180798724253 0.0287180798724253 0.609066234677835 0.0287180798724253]  ),
  SVector{5}([0.304779525704889 0.0287180798724253 0.609066234677835 0.0287180798724253 0.0287180798724253]  ),
  SVector{5}([0.304779525704889 0.609066234677835 0.0287180798724253 0.0287180798724253 0.0287180798724253]  ),
  SVector{5}([0.609066234677835 0.0287180798724253 0.0287180798724253 0.0287180798724253 0.304779525704889]  ),
  SVector{5}([0.609066234677835 0.0287180798724253 0.0287180798724253 0.304779525704889 0.0287180798724253]  ),
  SVector{5}([0.609066234677835 0.0287180798724253 0.304779525704889 0.0287180798724253 0.0287180798724253]  ),
  SVector{5}([0.609066234677835 0.304779525704889 0.0287180798724253 0.0287180798724253 0.0287180798724253]  ),
  SVector{5}([0.163259884255958 0.163259884255958 0.163259884255958 0.469494323021712 0.0407260242104143]    ),
  SVector{5}([0.163259884255958 0.163259884255958 0.163259884255958 0.0407260242104143 0.469494323021712]    ),
  SVector{5}([0.163259884255958 0.163259884255958 0.469494323021712 0.163259884255958 0.0407260242104143]    ),
  SVector{5}([0.163259884255958 0.163259884255958 0.469494323021712 0.0407260242104143 0.163259884255958]    ),
  SVector{5}([0.163259884255958 0.163259884255958 0.0407260242104143 0.163259884255958 0.469494323021712]    ),
  SVector{5}([0.163259884255958 0.163259884255958 0.0407260242104143 0.469494323021712 0.163259884255958]    ),
  SVector{5}([0.163259884255958 0.469494323021712 0.163259884255958 0.163259884255958 0.0407260242104143]    ),
  SVector{5}([0.163259884255958 0.469494323021712 0.163259884255958 0.0407260242104143 0.163259884255958]    ),
  SVector{5}([0.163259884255958 0.469494323021712 0.0407260242104143 0.163259884255958 0.163259884255958]    ),
  SVector{5}([0.163259884255958 0.0407260242104143 0.163259884255958 0.163259884255958 0.469494323021712]    ),
  SVector{5}([0.163259884255958 0.0407260242104143 0.163259884255958 0.469494323021712 0.163259884255958]    ),
  SVector{5}([0.163259884255958 0.0407260242104143 0.469494323021712 0.163259884255958 0.163259884255958]    ),
  SVector{5}([0.469494323021712 0.163259884255958 0.163259884255958 0.163259884255958 0.0407260242104143]    ),
  SVector{5}([0.469494323021712 0.163259884255958 0.163259884255958 0.0407260242104143 0.163259884255958]    ),
  SVector{5}([0.469494323021712 0.163259884255958 0.0407260242104143 0.163259884255958 0.163259884255958]    ),
  SVector{5}([0.469494323021712 0.0407260242104143 0.163259884255958 0.163259884255958 0.163259884255958]    ),
  SVector{5}([0.0407260242104143 0.163259884255958 0.163259884255958 0.163259884255958 0.469494323021712]    ),
  SVector{5}([0.0407260242104143 0.163259884255958 0.163259884255958 0.469494323021712 0.163259884255958]    ),
  SVector{5}([0.0407260242104143 0.163259884255958 0.469494323021712 0.163259884255958 0.163259884255958]    ),
  SVector{5}([0.0407260242104143 0.469494323021712 0.163259884255958 0.163259884255958 0.163259884255958]    ),
  SVector{5}([0.190976567672568 0.190976567672568 0.190976567672568 0.0371487518469173 0.389921545135378]    ),
  SVector{5}([0.190976567672568 0.190976567672568 0.190976567672568 0.389921545135378 0.0371487518469173]    ),
  SVector{5}([0.190976567672568 0.190976567672568 0.0371487518469173 0.190976567672568 0.389921545135378]    ),
  SVector{5}([0.190976567672568 0.190976567672568 0.0371487518469173 0.389921545135378 0.190976567672568]    ),
  SVector{5}([0.190976567672568 0.190976567672568 0.389921545135378 0.190976567672568 0.0371487518469173]    ),
  SVector{5}([0.190976567672568 0.190976567672568 0.389921545135378 0.0371487518469173 0.190976567672568]    ),
  SVector{5}([0.190976567672568 0.0371487518469173 0.190976567672568 0.190976567672568 0.389921545135378]    ),
  SVector{5}([0.190976567672568 0.0371487518469173 0.190976567672568 0.389921545135378 0.190976567672568]    ),
  SVector{5}([0.190976567672568 0.0371487518469173 0.389921545135378 0.190976567672568 0.190976567672568]    ),
  SVector{5}([0.190976567672568 0.389921545135378 0.190976567672568 0.190976567672568 0.0371487518469173]    ),
  SVector{5}([0.190976567672568 0.389921545135378 0.190976567672568 0.0371487518469173 0.190976567672568]    ),
  SVector{5}([0.190976567672568 0.389921545135378 0.0371487518469173 0.190976567672568 0.190976567672568]    ),
  SVector{5}([0.0371487518469173 0.190976567672568 0.190976567672568 0.190976567672568 0.389921545135378]    ),
  SVector{5}([0.0371487518469173 0.190976567672568 0.190976567672568 0.389921545135378 0.190976567672568]    ),
  SVector{5}([0.0371487518469173 0.190976567672568 0.389921545135378 0.190976567672568 0.190976567672568]    ),
  SVector{5}([0.0371487518469173 0.389921545135378 0.190976567672568 0.190976567672568 0.190976567672568]    ),
  SVector{5}([0.389921545135378 0.190976567672568 0.190976567672568 0.190976567672568 0.0371487518469173]    ),
  SVector{5}([0.389921545135378 0.190976567672568 0.190976567672568 0.0371487518469173 0.190976567672568]    ),
  SVector{5}([0.389921545135378 0.190976567672568 0.0371487518469173 0.190976567672568 0.190976567672568]    ),
  SVector{5}([0.389921545135378 0.0371487518469173 0.190976567672568 0.190976567672568 0.190976567672568]    ),
  SVector{5}([0.0270462772516200 0.0270462772516200 0.139237088063416 0.139237088063416 0.667433269369929]   ),
  SVector{5}([0.0270462772516200 0.0270462772516200 0.139237088063416 0.667433269369929 0.139237088063416]   ),
  SVector{5}([0.0270462772516200 0.0270462772516200 0.667433269369929 0.139237088063416 0.139237088063416]   ),
  SVector{5}([0.0270462772516200 0.139237088063416 0.0270462772516200 0.139237088063416 0.667433269369929]   ),
  SVector{5}([0.0270462772516200 0.139237088063416 0.0270462772516200 0.667433269369929 0.139237088063416]   ),
  SVector{5}([0.0270462772516200 0.139237088063416 0.139237088063416 0.0270462772516200 0.667433269369929]   ),
  SVector{5}([0.0270462772516200 0.139237088063416 0.139237088063416 0.667433269369929 0.0270462772516200]   ),
  SVector{5}([0.0270462772516200 0.139237088063416 0.667433269369929 0.0270462772516200 0.139237088063416]   ),
  SVector{5}([0.0270462772516200 0.139237088063416 0.667433269369929 0.139237088063416 0.0270462772516200]   ),
  SVector{5}([0.0270462772516200 0.667433269369929 0.0270462772516200 0.139237088063416 0.139237088063416]   ),
  SVector{5}([0.0270462772516200 0.667433269369929 0.139237088063416 0.0270462772516200 0.139237088063416]   ),
  SVector{5}([0.0270462772516200 0.667433269369929 0.139237088063416 0.139237088063416 0.0270462772516200]   ),
  SVector{5}([0.139237088063416 0.0270462772516200 0.0270462772516200 0.139237088063416 0.667433269369929]   ),
  SVector{5}([0.139237088063416 0.0270462772516200 0.0270462772516200 0.667433269369929 0.139237088063416]   ),
  SVector{5}([0.139237088063416 0.0270462772516200 0.139237088063416 0.0270462772516200 0.667433269369929]   ),
  SVector{5}([0.139237088063416 0.0270462772516200 0.139237088063416 0.667433269369929 0.0270462772516200]   ),
  SVector{5}([0.139237088063416 0.0270462772516200 0.667433269369929 0.0270462772516200 0.139237088063416]   ),
  SVector{5}([0.139237088063416 0.0270462772516200 0.667433269369929 0.139237088063416 0.0270462772516200]   ),
  SVector{5}([0.139237088063416 0.139237088063416 0.0270462772516200 0.0270462772516200 0.667433269369929]   ),
  SVector{5}([0.139237088063416 0.139237088063416 0.0270462772516200 0.667433269369929 0.0270462772516200]   ),
  SVector{5}([0.139237088063416 0.139237088063416 0.667433269369929 0.0270462772516200 0.0270462772516200]   ),
  SVector{5}([0.139237088063416 0.667433269369929 0.0270462772516200 0.0270462772516200 0.139237088063416]   ),
  SVector{5}([0.139237088063416 0.667433269369929 0.0270462772516200 0.139237088063416 0.0270462772516200]   ),
  SVector{5}([0.139237088063416 0.667433269369929 0.139237088063416 0.0270462772516200 0.0270462772516200]   ),
  SVector{5}([0.667433269369929 0.0270462772516200 0.0270462772516200 0.139237088063416 0.139237088063416]   ),
  SVector{5}([0.667433269369929 0.0270462772516200 0.139237088063416 0.0270462772516200 0.139237088063416]   ),
  SVector{5}([0.667433269369929 0.0270462772516200 0.139237088063416 0.139237088063416 0.0270462772516200]   ),
  SVector{5}([0.667433269369929 0.139237088063416 0.0270462772516200 0.0270462772516200 0.139237088063416]   ),
  SVector{5}([0.667433269369929 0.139237088063416 0.0270462772516200 0.139237088063416 0.0270462772516200]   ),
  SVector{5}([0.667433269369929 0.139237088063416 0.139237088063416 0.0270462772516200 0.0270462772516200]   ),
  SVector{5}([0.0366479097373171 0.0366479097373171 0.380941695489303 0.380941695489303 0.164820789546761]   ),
  SVector{5}([0.0366479097373171 0.0366479097373171 0.380941695489303 0.164820789546761 0.380941695489303]   ),
  SVector{5}([0.0366479097373171 0.0366479097373171 0.164820789546761 0.380941695489303 0.380941695489303]   ),
  SVector{5}([0.0366479097373171 0.380941695489303 0.0366479097373171 0.380941695489303 0.164820789546761]   ),
  SVector{5}([0.0366479097373171 0.380941695489303 0.0366479097373171 0.164820789546761 0.380941695489303]   ),
  SVector{5}([0.0366479097373171 0.380941695489303 0.380941695489303 0.0366479097373171 0.164820789546761]   ),
  SVector{5}([0.0366479097373171 0.380941695489303 0.380941695489303 0.164820789546761 0.0366479097373171]   ),
  SVector{5}([0.0366479097373171 0.380941695489303 0.164820789546761 0.0366479097373171 0.380941695489303]   ),
  SVector{5}([0.0366479097373171 0.380941695489303 0.164820789546761 0.380941695489303 0.0366479097373171]   ),
  SVector{5}([0.0366479097373171 0.164820789546761 0.0366479097373171 0.380941695489303 0.380941695489303]   ),
  SVector{5}([0.0366479097373171 0.164820789546761 0.380941695489303 0.0366479097373171 0.380941695489303]   ),
  SVector{5}([0.0366479097373171 0.164820789546761 0.380941695489303 0.380941695489303 0.0366479097373171]   ),
  SVector{5}([0.380941695489303 0.0366479097373171 0.0366479097373171 0.380941695489303 0.164820789546761]   ),
  SVector{5}([0.380941695489303 0.0366479097373171 0.0366479097373171 0.164820789546761 0.380941695489303]   ),
  SVector{5}([0.380941695489303 0.0366479097373171 0.380941695489303 0.0366479097373171 0.164820789546761]   ),
  SVector{5}([0.380941695489303 0.0366479097373171 0.380941695489303 0.164820789546761 0.0366479097373171]   ),
  SVector{5}([0.380941695489303 0.0366479097373171 0.164820789546761 0.0366479097373171 0.380941695489303]   ),
  SVector{5}([0.380941695489303 0.0366479097373171 0.164820789546761 0.380941695489303 0.0366479097373171]   ),
  SVector{5}([0.380941695489303 0.380941695489303 0.0366479097373171 0.0366479097373171 0.164820789546761]   ),
  SVector{5}([0.380941695489303 0.380941695489303 0.0366479097373171 0.164820789546761 0.0366479097373171]   ),
  SVector{5}([0.380941695489303 0.380941695489303 0.164820789546761 0.0366479097373171 0.0366479097373171]   ),
  SVector{5}([0.380941695489303 0.164820789546761 0.0366479097373171 0.0366479097373171 0.380941695489303]   ),
  SVector{5}([0.380941695489303 0.164820789546761 0.0366479097373171 0.380941695489303 0.0366479097373171]   ),
  SVector{5}([0.380941695489303 0.164820789546761 0.380941695489303 0.0366479097373171 0.0366479097373171]   ),
  SVector{5}([0.164820789546761 0.0366479097373171 0.0366479097373171 0.380941695489303 0.380941695489303]   ),
  SVector{5}([0.164820789546761 0.0366479097373171 0.380941695489303 0.0366479097373171 0.380941695489303]   ),
  SVector{5}([0.164820789546761 0.0366479097373171 0.380941695489303 0.380941695489303 0.0366479097373171]   ),
  SVector{5}([0.164820789546761 0.380941695489303 0.0366479097373171 0.0366479097373171 0.380941695489303]   ),
  SVector{5}([0.164820789546761 0.380941695489303 0.0366479097373171 0.380941695489303 0.0366479097373171]   ),
  SVector{5}([0.164820789546761 0.380941695489303 0.380941695489303 0.0366479097373171 0.0366479097373171]   )
]

const pentatope126wn = [
  0.04654868655484498,
  0.0014672641241909648,
  0.0014672641241909648,
  0.0014672641241909648,
  0.0014672641241909648,
  0.0014672641241909648,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.0029906128790727765,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.009907301984077122,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.011320943419643296,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.0046932055563252825,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926,
  0.010698055349285926
]

const pentatopeAn = [
  pentatope1an,
  pentatope5an,
  pentatope15an,
  pentatope35an,
  pentatope70an,
  pentatope126an
]

const pentatopeWn = [
  pentatope1wn,
  pentatope5wn,
  pentatope15wn,
  pentatope35wn,
  pentatope70wn,
  pentatope126wn
]


text:
export CalculateSampling,
Convmtx,
DotTest,
InnerProduct,
Integrate,
MeasureSNR,
Pad5D,
PadFirstAxis,
PhaseShift,
TriangleFilter
include("CalculateSampling.jl")
include("Convmtx.jl")
include("DotTest.jl")
include("InnerProduct.jl")
include("Integrate.jl")
include("MeasureSNR.jl")
include("Pad5D.jl")
include("PadFirstAxis.jl")
include("PhaseShift.jl")
include("TriangleFilter.jl")

ext: jl

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
module TravelingSalesmanHeuristics
using Graphs

include("simulated_annealing.jl")

export solve_tsp, lowerbound, nearest_neighbor, cheapest_insertion, simulated_annealing, two_opt

"""
.. solve_tsp(distmat) ..

One-line interface to approximately solving a TSP by specifying a distance matrix.
This method provides fairly quick solutions but no extra control. For more fine-grained
control over the heuristics used, try nearest_neighbor or cheapest_insertion.
"""
function solve_tsp{T<:Real}(distmat::Matrix{T})
	p1, c1 = nearest_neighbor(distmat)
	p2, c2 = cheapest_insertion(distmat)
	if c1 < c2
		return p1, c1
	else
		return p2, c2
	end
end

###
# path generation heuristics
###

"""
Approximately solve a TSP using the nearest neighbor heuristic. You must pass a square matrix
distmat where distmat[i,j] represents the distance from city i to city j. The matrix needn't be
symmetric and possibly could contain negative values, though nonpositive values have not been tested.

Optional arguments:

firstCity: specifies the city to begin the path on. An empty Nullable corresponds to random selection.
	This argument is ignored if repetitive = true. Defaults to an empty Nullable{Int}

repetitive: boolean for whether to try starting from all possible cities, keeping the best. Defaults to false.

closepath: boolean for whether to include the arc from the last city visited back to the first city in
	cost calculations. If true, the first city appears first and last in the path. Defaults to true.

do2opt: whether to refine the path found by 2-opt switches (corresponds to removing path crossings in
	the planar Euclidean case). Defaults to true.

returns a tuple (path, pathcost) where path is a Vector{Int} corresponding to the order of cities visited
"""
function nearest_neighbor{T<:Real}(distmat::Matrix{T};
							   firstcity::Nullable{Int} = Nullable{Int}(),
							   repetitive = false,
							   closepath = true,
							   do2opt = true)
	# must have a square matrix
	check_square(distmat, "Must pass a square distance matrix to nearest_neighbor")

	numCities = size(distmat, 1)

	# if repetitive, we do all possible cities, and pick the best
	if repetitive
		function nnHelper(i)
			nearest_neighbor(distmat,
						  firstcity = Nullable(i),
						  closepath = closepath,
						  do2opt = do2opt,
						  repetitive = false)
		end
		# do nn for each startin city
		results = map(nnHelper, collect(1:numCities))
		# pick out lowest cost
		_, bestInd = findmin(map(res -> res[2], results))
		return results[bestInd]
	end

	# if not repetitive, we actually perform the heuristic for one starting city

	# put first city on path
	path = Vector{Int}()
	if isnull(firstcity)
		firstcity = rand(1:numCities)
	else
		firstcity = get(firstcity)
	end
	push!(path, firstcity)

	# cities to visit
	citiesToVisit = collect(1:(firstcity - 1))
	append!(citiesToVisit, collect((firstcity + 1):numCities))

	# nearest neighbor loop
	while !isempty(citiesToVisit)
		curCity = path[end]
		dists = distmat[curCity, citiesToVisit]
		_, nextInd = findmin(dists)
		nextCity = citiesToVisit[nextInd]
		push!(path, nextCity)
		deleteat!(citiesToVisit, nextInd)
	end

	# complete cycle? (duplicates first city)
	if closepath
		push!(path, firstcity)
	end

	# do swaps?
	if do2opt
		path, _ = two_opt(distmat, path)
	end

	return (path, pathcost(distmat, path))
end

"""
.. cheapest_insertion(distmat, initpath) ..

Given a distance matrix and an initial path, complete the tour by
repeatedly doing the cheapest insertion. The initial path must have length at least 2, but can be
simply [i, i] for some city index i which corresponds to starting with a self-loop at city i.
Insertions are always in the interior of the current path so this heuristic can also be used for
non-closed TSP paths.
Currently the implementation is a naive n^3 algorithm.
"""
function cheapest_insertion{T<:Real}(distmat::Matrix{T}, initpath::Vector{Int})
	check_square(distmat, "Distance matrix passed to cheapest_insertion must be square.")

	n = size(distmat, 1)
	path = copy(initpath)

	# collect cities to visited
	visitus = setdiff(collect(1:n), initpath)

	# helper for insertion cost
	# tour cost change for inserting node k after the node at index after in the path
	function inscost(k, after)
		return distmat[path[after], k] +
			  distmat[k, path[after + 1]] -
			  distmat[path[after], path[after + 1]]
	end

	while !isempty(visitus)
		bestCost = Inf
		bestInsertion = (-1, -1)
		for k in visitus
			for after in 1:(length(path) - 1) # can't insert after end of path
				c = inscost(k, after)
				if c < bestCost
					bestCost = c
					bestInsertion = (k, after)
				end
			end
		end
		# bestInsertion now holds (k, after)
		# insert into path, remove from to-do list
		k, after = bestInsertion
		insert!(path, after + 1, k)
		visitus = setdiff(visitus, k)
	end

	return (path, pathcost(distmat, path))
end
"""
.. cheapest_insertion(distmat; ...) ..

Cheapest insertion with a self-loop as the initial path.
distmat must be a square real matrix. Non-symmetric distance matrices should work. Negative
distances have not been tested but may also work.

Optional arguments:

firstcity: Specifies which city should have a self loop as the initial path. Nullable(i) for
city i or an empty Nullable{Int} for random selection. This argument is ignored if
repetitive = true. Defaults to random selection.

repetitive: boolean for whether to try starting from all possible cities, keeping the best. Defaults to false.

do2opt: boolean for whether to improve the paths found by 2-opt swaps. Defaults to true.
"""
function cheapest_insertion{T<:Real}(distmat::Matrix{T};
								 firstcity::Nullable{Int} = Nullable{Int}(),
							     repetitive::Bool = false,
								 do2opt::Bool = true)
	#infer size
	n = size(distmat, 1)

	# if repetitive, we do all possible cities, and pick the best
	if repetitive
		function ciHelper(i)
			cheapest_insertion(distmat,
						  firstcity = Nullable(i),
						  do2opt = do2opt,
						  repetitive = false)
		end
		# do cheapest insertion for each starting city
		results = map(ciHelper, collect(1:n))
		# pick out lowest cost
		_, bestInd = findmin(map(res -> res[2], results))
		return results[bestInd]
	end

	# okay, we're not repetitive. Do we pick first city at random?
	firstcity = isnull(firstcity) ? rand(1:n) : get(firstcity)
	path, cost = cheapest_insertion(distmat, [firstcity, firstcity])

	# user may have asked for 2-opt refinement
	if do2opt
		path, cost = two_opt(distmat, path)
	end

	return path, cost
end

###
# helpers
###

# make sure a passed distance matrix is a square
function check_square(m, msg)
	if size(m, 1) != size(m, 2)
		error(msg)
	end
end

# helper for readable one-line path costs
# optionally specify the bounds for the subpath we want the cost of
# defaults to the whole path
# but when calculating reversed path costs can help to have subpath costs
function pathcost{T<:Real}(distmat::Matrix{T}, path::Vector{Int}, lb::Int = 1, ub::Int = length(path))
	cost = zero(T)
	for i in lb:(ub - 1)
		@inbounds cost += distmat[path[i], path[i+1]]
	end
	return cost
end
# calculate the cost of reversing part of a path
# cost of walking along the entire path specified but reversing the sequence from revLow to revHigh, inclusive
function pathcost_rev{T<:Real}(distmat::Matrix{T}, path::Vector{Int}, revLow::Int, revHigh::Int)
	cost = zero(T)
	# if there's an initial unreversed section
	if revLow > 1
		for i in 1:(revLow - 2)
			@inbounds cost += distmat[path[i], path[i+1]]
		end
		# from end of unreversed section to beginning of reversed section
		@inbounds cost += distmat[path[revLow - 1], path[revHigh]]
	end
	# main reverse section
	for i in revHigh:-1:(revLow + 1)
		@inbounds cost += distmat[path[i], path[i-1]]
	end
	# if there's an unreversed section after the reversed bit
	n = length(path)
	if revHigh < length(path)
		# from end of reversed section back to regular
		@inbounds cost += distmat[path[revLow], path[revHigh + 1]]
		for i in (revHigh + 1):(n-1)
			@inbounds cost += distmat[path[i], path[i+1]]
		end
	end
	return cost
end

# helper to get min spanning trees
# returns a (n-1) long Vector of Tuple{Int, Int} where each tuple is an edge in the MST
# and the total weight of the tree
function minspantree(distmat)
	if size(distmat, 1) != size(distmat, 2)
		error("Distance matrix passed to minspantree must be square.")
	end
	n = size(distmat, 1)
	numEdge = convert(Int, n * (n-1)) # not / 2 since directed graph

	# construct graph and list of edge weights
	g = simple_complete_graph(n)
	eWeights = zeros(numEdge)

	for i in eachindex(edges(g))
		e = edges(g)[i]
		eWeights[i] = distmat[source(e), target(e)]
	end

	(treeEdges, treeWeights) = kruskal_minimum_spantree(g, eWeights)

	return (map(e -> (source(e), target(e)), treeEdges),
		    sum(treeWeights))
end

###
# path improvement heuristics
###

"perform 2-opt reversals until doing so does not improve the path cost

First argument is the distance matrix, second is the path to be improved."
function two_opt{T<:Real}(distmat::Matrix{T}, path::Vector{Int})
	# size checks
	n = length(path)
	if size(distmat, 1) != size(distmat, 2)
		error("Distance matrix passed to two_opt must be square.")
	end

	# don't modify input
	path = copy(path)

	# main loop
	# check every possible switch until no 2-swaps reduce objective
	# if the path passed in is a loop (first/last nodes are the same)
	# then we must keep these the endpoints of the path the same
	# ie just keep it a loop, and therefore it doesn't matter which node is at the end
	# if the path is not a cycle, we can do any reversing we like
	isCycle = path[1] == path[end]
	switchLow = isCycle ? 2 : 1
	switchHigh = isCycle ? n - 1 : n
	prevCost = Inf
	curCost = pathcost(distmat, path)
	while prevCost > pathcost(distmat, path)
		prevCost = curCost
		# we can't change the first
		for i in switchLow:(switchHigh-1)
			for j in (i+1):switchHigh
				altCost = pathcost_rev(distmat, path, i, j)
				if altCost < curCost
					curCost = altCost
					reverse!(path, i, j)
				end
			end
		end
	end
	return path, pathcost(distmat, path)
end

###
# bounds
###

# I implement several simple but super loose bounds
# These only apply to closed tours

# the cost of a tour must be >= the sum over all vertices of
# the cost of the cheapest edge leaving that vertex
# likewise for the cheapest edge entering that vertex
# since we must go to and leave each vertex
function vertwise_bound(distmat)
	# the simple code below would tend to pick out the 0 costs on the diagonal
	# so make a doctored copy of the distance matrix with high costs on the diagonal
	m = maximum(distmat)
	distmat_nodiag = distmat + m * eye(distmat)
	leaving = sum(minimum(distmat_nodiag, 2))
	entering = sum(minimum(distmat_nodiag, 1))
	return maximum([leaving, entering])
end

# a simplified/looser version of Held-Karp bounds
# any tour is a spanning tree on (n-1) verts plus two edges
# from the left out vert
# so the maximum over all verts (as the left out vert) of
# MST cost on remaining vertices plus 2 cheapest edges from
# the left out vert is a lower bound
function hkinspired_bound(distmat)
	n = size(distmat, 1)
	if size(distmat, 2) != n
		error("Must pass square distance matrix to hkinspired_bound")
	end
	function del_vert(v)
		keep = [1:(v-1) ; (v+1):n]
		return distmat[keep, keep]
	end

	# make sure min(distmat[v,:]) doesn't pick diagonal elements
	m = maximum(distmat)
	distmat_nodiag = distmat + m * eye(distmat)
	function cost_leave_out(v)
		dmprime = del_vert(v)
		_, c = minspantree(dmprime)
		c += minimum(distmat_nodiag[v,:])
		c += minimum(distmat_nodiag[:,v])
		return c
	end

	return maximum(map(cost_leave_out, 1:n))
end

# best lower bound we have
"""
Lower bound the cost of the optimal TSP tour. At present, the bounds considered
are a simple bound based on the minimum cost of entering and exiting each city and
a slightly better bound inspired by the Held-Karp bounds; note that the implementation
here is simpler and less tight than proper HK bounds.
"""
function lowerbound(distmat)
	return maximum([
			vertwise_bound(distmat),
			hkinspired_bound(distmat)
			])
end

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PTSCOTCH_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PTSCOTCH")
JLLWrappers.@generate_main_file("PTSCOTCH", UUID("b3ec0f5a-9838-5c9b-9e77-5f2c6a4b089f"))
end  # module PTSCOTCH_jll


text:
include("../include.jl")
include("create_report.jl")
overlapping_results = get_CUTEst_results()

##
## FUNCTION VALUES
##
problem_list = collect(keys(first(overlapping_results)[2]))
method_list = collect(keys(overlapping_results))

CON_TOL = 1e-6
F_TOL = 0.5

winners = Dict()
winners["tie"] = 0
winners["fail"] = 0
begin
  for method_name in method_list
      print(pd(""), pd(method_name))
      winners[method_name] = 0
  end
  print("\n")

  for problem_name in problem_list
      print(pd(problem_name))
      best = "fail"
      both_kkt = true
      fbest = 1e14
      for method_name in method_list
        info = overlapping_results[method_name][problem_name]
        fval = info.fval
        con_vio = info.con_vio
        print(rd(con_vio))
        print(rd(fval))

        if info.status == :Optimal
          if con_vio > CON_TOL
            both_kkt = false
            println("CON VIOLATION!!")
          else
            ftol_scaled = F_TOL * (1.0 + abs(fbest))
            if fval <= fbest - ftol_scaled
              fbest = fval
              best = method_name
            elseif fbest - ftol_scaled < fval && fval <= fbest + ftol_scaled
              fbest = fval
              best = "tie"
            end
          end
        else
          both_kkt = false
        end
      end
      if both_kkt
        winners[best] += 1
      end
      print(pd(best))
      print("\n")
  end
end


text:
module MixtureGibbs
using Distributions
using Plots
using StatsBase

include("gibbs.jl")

function update_suffstats!(state)
    state.suffstats = sufficient_statistics(state.data, state.assignment)
    return nothing
end

function log_assignment_score(data_id, cluster_id, state)
    """log p(z_i=k ,|, cdot)
    """
    x = state.data[data_id]
    θ = state.θ[cluster_id]
    var = state.v
    log_pi = log.(state.π[cluster_id])
    return log_pi + logpdf(Normal(θ, var), x)
end

function assigment_probs(data_id, state)
    """p(z_i=cid ,|, cdot) for cid in cluster_ids
    """
    scores = map(cid -> log_assignment_score(data_id, cid, state), state.cluster_ids)
    scores = exp.(scores)
    return scores / sum(scores)
end

function sample_assignment(data_id, state)
    """Sample cluster assignment for data_id given current state

    cf Step 1 of Algorithm 2.1 in Sudderth 2006
    """
    p = assigment_probs(data_id, state)
    return state.cluster_ids[rand(Categorical(p))]
end

function update_assignment!(state)
    """Update cluster assignment for each data point given current state

    cf Step 1 of Algorithm 2.1 in Sudderth 2006
    """
    for (data_id, x) in enumerate(state.data)
        state.assignment[data_id] = sample_assignment(data_id, state)
        update_suffstats!(state)
    end
end

function sample_mixture_weights(state::GibbsState)
    """Sample new mixture weights from current state according to
    a Dirichlet distribution

    cf Step 2 of Algorithm 2.1 in Sudderth 2006
    """
    ss = state.suffstats
    n_clusters = length(state.cluster_ids)
    alpha = [ss[cid].count + state.hp.α / n_clusters for cid in state.cluster_ids]
    return rand(Dirichlet(alpha))
end

function update_mixture_weights!(state::GibbsState)
    """Update state with new mixture weights from current state
    sampled according to a Dirichlet distribution

    cf Step 2 of Algorithm 2.1 in Sudderth 2006
    """
    state.π = sample_mixture_weights(state)
    return nothing
end

function sample_cluster_mean(cluster_id, state)
    cluster_var = state.v
    hp_mean = state.hp.λ1
    hp_var = state.hp.λ2
    ss = state.suffstats[cluster_id]

    numerator = hp_mean / hp_var + ss.μ * ss.count / cluster_var
    denominator = (1.0 / hp_var + ss.count / cluster_var)
    posterior_mu = numerator / denominator
    posterior_var = 1.0 / denominator

    return randn() * sqrt(posterior_var) + posterior_mu
end

function update_cluster_means!(state::GibbsState)
    state.θ = [sample_cluster_mean(cid, state) for cid in state.cluster_ids]
    return nothing
end

function gibb_step!(state::GibbsState)
    update_assignment!(state)
    update_mixture_weights!(state)
    update_cluster_means!(state)
    return nothing
end

function plot_cluster!(state::GibbsState; opts...)
    data = state.data
    assignment = state.assignment
    histogram!(data; group=assignment, alpha=0.5, opts...)
    ylims!(0, 150)
    xlims!(-8, 8)
    return nothing
end

function run()
    n_clusters = 3
    data = generate(1000)
    state = GibbsState(data, n_clusters)

    n_steps = 5
    fig = plot(; legend=false, layout=(n_steps, 1), size=(400, 600))
    for i in 1:n_steps
        plot_cluster!(state; subplot=i)
        annotate!(-5, 120, "step = $((i-1))"; subplot=i)
        gibb_step!(state)
    end
    display(fig)

    return nothing
end

function animation()
    n_clusters = 3
    data = generate(1000)
    state = GibbsState(data, n_clusters)

    anim = @animate for i in 1:50
        plot(; legend=false)
        plot_cluster!(state)
        annotate!(-5, 120, "step = $((i-1))")
        gibb_step!(state)
    end

    return gif(anim, "assets/mixture_gibbs.gif"; fps=1)
end

end


text:
fitting3K = DataFrame(fittingType = ["PIPE",
"ELL-THD-STD",
"ELL-THD-LR",
"ELL-FLG-RD1",
"ELL-FLG-RD2",
"ELL-FLG-RD4",
"ELL-FLG-RD6",
"ELL-MIT-1",
"ELL-MIT-2",
"EL-MIT-3",
"EL45-THD-STD",
"EL45-LR",
"EL45-MIT-1",
"EL45-MIT-2",
"EL180-THD-1",
"EL180-FLG-RD1",
"EL180-FLG-RD15",
"TEE-THD-RD1",
"TEE-THD-RD15",
"TEE-FLG",
"TEE-STUB",
"TEE-THD-THRU",
"TEE-FLG-THRU",
"TEE-STUB-THRU",
"VALVE-ANGLE-45-FULL",
"VALVE-ANGLE-90-FULL",
"VALVE-GLOBE-STD-B1",
"VALVE-PLUG-BRANCH",
"VALVE-PLUG-THRU",
"VALVE-PLUG-3WAY-THRU",
"VALVE-GATE-STD-B1",
"VALVE-BALL-STD-B1",
"VALVE-DIAPH-DAM",
"VALVE-CHECK-SWING",
"VALVE-CHECK-LIFT",
"ENTRANCE-NON-RECOVER",
"EXIT-NON-RECOVER"],
Description = ["PipePipePipe",
"Elbow - 90 DegreeThreaded-Standardr/D = 1",
"Elbow - 90 DegreeThreaded Long Radiusr/D = 1.5",
"Elbow - 90 DegreeFlanged Welded Bendr/D = 1",
"Elbow - 90 DegreeFlanged Welded Bendr/D = 2",
"Elbow - 90 DegreeFlanged Welded Bendr/D = 4",
"Elbow - 90 DegreeFlanged Welded Bendr/D = 6",
"Elbow - 90 DegreeMitered 1 Weld 90 Degree",
"Elbow - 90 DegreeMitered 2 Welds 45 Degree",
"Elbow - 90 DegreeMitered 3 Welds 30 Degree",
"Elbow - 45 DegreeThreaded Standardr/D = 1",
"Elbow - 45 DegreeLong Radiusr/D = 1.5",
"Elbow - 45 DegreeMitered 1 Weld 45 Degree",
"Elbow - 45 DegreeMitered 2 Welds 22.5 Degree",
"Elbows - 180 DegreeThreaded Close Return Bendr/D = 1",
"Elbows - 180 DegreeFlangedr/D = 1",
"Elbows - 180 DegreeAllr/D = 1.5",
"TEEThreadedr/D = 1",
"TEEThreadedr/D = 1.5",
"TEEFlangedr/D = 1",
"TEEStub-in Branch",
"TEERun Through Threadedr/D = 1",
"TEEFlangedr/D = 1",
"TEEStub-in Branch",
"ValveAngle - 45 DegreeFull Line Size b = 1",
"ValveAngle - 90 DegreeFull Line Size b = 1",
"ValveGlobeStandard b = 1",
"ValvePlugBranch Flow",
"ValvePlugStraight Through",
"ValvePlugThree-Way (flow through)",
"ValveGateStandard b = 1",
"ValveBallStandard b = 1",
"ValveDiaphragmDam-Type",
"ValveSwing CheckVmin = 35r-1/2",
"ValveLift CheckVmin = 40r-1/2",
"InletNon-recoverable",
"ExitNon-recoverable"],
Fitting = ["Pipe",
"Elbow - 90 Degree",
"Elbow - 90 Degree",
"Elbow - 90 Degree",
"Elbow - 90 Degree",
"Elbow - 90 Degree",
"Elbow - 90 Degree",
"Elbow - 90 Degree",
"Elbow - 90 Degree",
"Elbow - 90 Degree",
"Elbow - 45 Degree",
"Elbow - 45 Degree",
"Elbow - 45 Degree",
"Elbow - 45 Degree",
"Elbows - 180 Degree",
"Elbows - 180 Degree",
"Elbows - 180 Degree",
"TEE",
"TEE",
"TEE",
"TEE",
"TEE",
"TEE",
"TEE",
"Valve",
"Valve",
"Valve",
"Valve",
"Valve",
"Valve",
"Valve",
"Valve",
"Valve",
"Valve",
"Valve",
"Inlet",
"Exit"],
Type = ["Pipe",
"Threaded-Standard",
"Threaded Long Radius",
"Flanged Welded Bend",
"Flanged Welded Bend",
"Flanged Welded Bend",
"Flanged Welded Bend",
"Mitered 1 Weld 90 Degree",
"Mitered 2 Welds 45 Degree",
"Mitered 3 Welds 30 Degree",
"Threaded Standard",
"Long Radius",
"Mitered 1 Weld 45 Degree",
"Mitered 2 Welds 22.5 Degree",
"Threaded Close Return Bend",
"Flanged",
"All",
"Threaded",
"Threaded",
"Flanged",
"Stub-in Branch",
"Run Through Threaded",
"Flanged",
"Stub-in Branch",
"Angle - 45 Degree",
"Angle - 90 Degree",
"Globe",
"Plug",
"Plug",
"Plug",
"Gate",
"Ball",
"Diaphragm",
"Swing Check",
"Lift Check",
"",
""],
Geometry = ["Pipe",
"r/D = 1",
"r/D = 1.5",
"r/D = 1",
"r/D = 2",
"r/D = 4",
"r/D = 6",
"",
"",
"",
"r/D = 1",
"r/D = 1.5",
"",
"",
"r/D = 1",
"r/D = 1",
"r/D = 1.5",
"r/D = 1",
"r/D = 1.5",
"r/D = 1",
"",
"r/D = 1",
"r/D = 1",
"",
"Full Line Size b = 1",
"Full Line Size b = 1",
"Standard b = 1",
"Branch Flow",
"Straight Through",
"Three-Way (flow through)",
"Standard b = 1",
"Standard b = 1",
"Dam-Type",
"Vmin = 35r-1/2",
"Vmin = 40r-1/2",
"Non-recoverable",
"Non-recoverable"],
K1 = [0,
800,
800,
800,
800,
800,
800,
1000,
800,
800,
500,
500,
500,
500,
1000,
1000,
1000,
500,
800,
800,
1000,
200,
150,
100,
950,
1000,
1500,
500,
300,
300,
300,
300,
1000,
1500,
2000,
0,
0],
Kinf = [0,
0.14,
0.071,
0.091,
0.056,
0.066,
0.075,
0.27,
0.068,
0.035,
0.071,
0.052,
0.086,
0.052,
0.23,
0.12,
0.1,
0.274,
0.14,
0.28,
0.34,
0.091,
0.017,
0,
0.25,
0.69,
1.7,
0.41,
0.084,
0.14,
0.037,
0.017,
0.69,
0.46,
2.85,
1.5,
0],
Kd = [0,
4,
4.2,
4,
3.9,
3.9,
4.2,
4,
4.1,
4.2,
4.2,
4,
4,
4,
4,
4,
4,
4,
4,
4,
4,
4,
4,
0,
4,
4,
3.6,
4,
3.9,
4,
3.9,
4,
4.9,
4,
3.8,
0,
0],
Kp = [1,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0] )


text:
## This turns out to be exactly the same as the symplectit Euler method
# L(xck,vck) -> 1/2 Δt (Ld(xdk,(xdk+1-xdk)/Δt) + Ld(xdk+1,(xdk+1-xdk)/Δt))
# L(qck,ωck) -> 1/2 Δt (Ld(qdk,2 V qdk† (qdk+1-qdk)/Δt) + Ld(qdk+1,2 V qdk† (qdk+1-qdk)/Δt)
# ωckw = sqrt((2/Δt)^2 - ωckᵀωck) - 2/Δt
# Fckᵀxck -> 1/2 (Fdkᵀxdk + Fdk+1ᵀxdk)

METHODORDER = 2
getGlobalOrder() = (global METHODORDER; return METHODORDER)

# Convenience functions
@inline getx3(state::State, Δt) = state.xk[2] + state.vsol[2]*Δt
@inline getq3(state::State, Δt) = state.qk[2] * ωbar(state.ωsol[2],Δt) * Δt / 2

@inline posargsc(state::State) = (state.xc, state.qc)
@inline fullargsc(state::State) = (state.xc, state.vc, state.qc, state.ωc)
@inline posargsk(state::State; k=2) = (state.xk[k], state.qk[k])
@inline posargssol(state::State) = (state.xsol[2], state.qsol[2])
@inline fullargssol(state::State) = (state.xsol[2], state.vsol[2], state.qsol[2], state.ωsol[2])
@inline posargsnext(state::State, Δt) = (getx3(state, Δt), getq3(state, Δt))

@inline function derivωbar(ω::SVector{3}, Δt)
    msq = -sqrt(4 / Δt^2 - dot(ω, ω))
    return [ω' / msq; I]
end

@inline function ωbar(ω, Δt)
    return UnitQuaternion(sqrt(4 / Δt^2 - dot(ω, ω)), ω, false)
end

@inline function setForce!(state::State, F, τ)
    state.Fk[2] = F
    state.τk[2] = τ
    return
end

@inline function discretizestate!(body::Body{T}, Δt) where T
    state = body.state
    xc = state.xc
    qc = state.qc
    vc = state.vc
    ωc = state.ωc

    state.xk[1] = xc
    state.xk[2] = xc + vc*Δt
    state.qk[1] = qc
    state.qk[2] = qc * ωbar(ωc,Δt) * Δt / 2

    state.Fk[1] = szeros(T,3)
    state.Fk[2] = szeros(T,3)
    state.τk[1] = szeros(T,3)
    state.τk[2] = szeros(T,3)

    return
end

@inline function currentasknot!(body::Body)
    state = body.state

    state.xk[2] = state.xc
    state.qk[2] = state.qc

    return
end

@inline function updatestate!(body::Body{T}, Δt) where T
    state = body.state

    state.xc = state.xsol[2]
    state.qc = state.qsol[2]
    state.vc = state.vsol[2]
    state.ωc = state.ωsol[2]

    state.xk[1] = state.xk[2]
    state.xk[2] = state.xk[2] + state.vsol[2]*Δt
    state.qk[1] = state.qk[2]
    state.qk[2] = state.qk[2] * ωbar(state.ωsol[2],Δt) * Δt / 2

    state.xsol[2] = state.xk[2]
    state.qsol[2] = state.qk[2]

    state.Fk[1] = state.Fk[2]
    state.Fk[2] = szeros(T,3)
    state.τk[1] = state.τk[2]
    state.τk[2] = szeros(T,3)
    return
end

@inline function setsolution!(body::Body)
    state = body.state
    state.xsol[2] = state.xk[2]
    state.qsol[2] = state.qk[2]
    state.vsol[1] = state.vc
    state.vsol[2] = state.vc
    state.ωsol[1] = state.ωc
    state.ωsol[2] = state.ωc
    return
end

@inline function settempvars!(body::Body{T}, x, v, F, q, ω, τ, d) where T
    state = body.state
    stateold = deepcopy(state)

    state.xc = x
    state.qc = q
    state.vc = v
    state.ωc = ω
    state.Fk[1] = F
    state.Fk[2] = F
    state.τk[1] = τ
    state.τk[2] = τ
    state.d = d

    return stateold
end


text:
@testset "generic" begin
    k = SqExponentialKernel()
    @test length(k) == 1
    @test iterate(k) == (k,nothing)
    @test iterate(k,1) == nothing
end


text:
using BenchmarkTools
using Random
using RANSAC

const SUITE = BenchmarkGroup()

SUITE["RANSAC"] = BenchmarkGroup()

function setupme(iterations)
    vs, ns, norms4Plot, shape_s = examplepc3()
    pcr = RANSACCloud(vs, ns, 32)
    #p = RANSACParameters{Float64}()
    # plane
    # p = RANSACParameters(p, ϵ_plane=0.3, α_plane=deg2rad(5))
    p_plane = (ϵ=0.3, α=deg2rad(5),)
    # cylidner
    # p = RANSACParameters(p, ϵ_cylinder=0.3, α_cylinder=deg2rad(5))
    p_cyl = (ϵ=0.3, α=deg2rad(5),)
    # sphere
    # p = RANSACParameters(p, ϵ_sphere=0.3, α_sphere=deg2rad(5))
    p_sp = (ϵ=0.3, α=deg2rad(5),)
    # number of minimal subsets drawed in one iteration
	# probability that we found shapes
	# minimum shape size
	# maximum number of iteration
    p_it = (minsubsetN=15, prob_det=0.9, τ=900, itermax=iterations, drawN=3)
    p = ransacparameters(;plane=p_plane, cylinder=p_cyl, sphere=p_sp, iteration=p_it)
    return (pcr, p, true)
end

Random.seed!(1234);
sharp = setupme(100_000)
SUITE["RANSAC"]["smallideal"] = @benchmarkable ransac($sharp..., reset_rand=true)


text:
using JITrench
using Random
using Printf
using Plots

Random.seed!(10)

function generate_datset(N)
    x = rand(N) * 10
    y = 2x .+ 3 .+ (rand(N) .* 5)
    return x, y
end

function predict(x, W, b)
    y = JITrench.matmul(x, W) .+ b
    return y
end

function mse(y_true, y_pred)
    diff = y_true .- y_pred
    return JITrench.sum(diff .^ 2) / length(diff.values)
end

function train(x, y, W_init, b_init, n_iter; lr=0.01)
    x = Variable(x)
    y = Variable(y)
    W = Variable(W_init)
    b = Variable(b_init)
    history = []
    for iter in 1:n_iter
        y_pred = predict(x, W, b)
        loss = mse(y, y_pred)
        JITrench.cleargrad!(W)
        JITrench.cleargrad!(b)
        backward!(loss)
        W.values -= W.grad.values * lr
        b.values -= b.grad.values * lr
        push!(history, loss.values)
        @printf "iters %4i [loss] %.2f [W] %.2f [b] %.2f\n" iter loss.values W.values b.values
    end
    return W, b, history
end

x, y = generate_datset(100)
W = rand()
b = rand()
n_iters = 100

W_trained, b_trained, history = train(x, y, W, b, n_iters)

println("finish train.")
@show W_trained
@show b_trained

plt = scatter(x, y, label=["data"])
plot!(plt, 1:0.1:10, predict(Variable(collect(1:0.1:10)), W_trained, b_trained).values, label=["predict"])
savefig("predict.png")

plot(history, title="learning curve", yaxis=:log10)
savefig("history.png")


text:
using Revise
	using SparseArrays, LinearAlgebra, DiffEqOperators, Setfield, Parameters
	using BifurcationKit
	using Plots
	const BK = BifurcationKit
################################################################################
# case of the SH equation
Nx = 200; Lx = 6.;
X = -Lx .+ 2Lx/Nx*(0:Nx-1) |> collect
hx = X[2]-X[1]

norminf(x) = norm(x, Inf64)
const _weight = rand(Nx)
normweighted(x) = norm(_weight .* x)

# Q = Neumann0BC(hx)
Q = Dirichlet0BC(hx |> typeof)
Dxx = sparse(CenteredDifference(2, 2, hx, Nx) * Q)[1]
Lsh = -(I + Dxx)^2

function R_SH(u, par)
	@unpack l, ν, L1 = par
	out = similar(u)
	out .= L1 * u .+ l .* u .+ ν .* u.^3 - u.^5
end

Jac_sp(u, par) = par.L1 + spdiagm(0 => par.l .+ 3 .* par.ν .* u.^2 .- 5 .* u.^4)
d2R(u,p,dx1,dx2) = @. p.ν * 6u*dx1*dx2 - 5*4u^3*dx1*dx2
d3R(u,p,dx1,dx2,dx3) = @. p.ν * 6dx3*dx1*dx2 - 5*4*3u^2*dx1*dx2*dx3
jet = (R_SH, Jac_sp, d2R, d3R)

parSH = (l = -0.1, ν = 2., L1 = Lsh)
####################################################################################################
sol0 = 1.1cos.(X) .* exp.(-0X.^2/(2*5^2))
	optnew = NewtonPar(verbose = false, tol = 1e-12)
	# allocations 357, 0.8ms
	sol1, = @time newton(
	R_SH, Jac_sp,
	sol0, (@set parSH.l = -.1), optnew)
	Plots.plot(X, sol1)

opts = ContinuationPar(dsmin = 0.0001, dsmax = 0.01, ds = 0.01,
		newtonOptions = setproperties(optnew; maxIter = 30, tol = 1e-8), pMax = 1.,
		maxSteps = 300, plotEveryStep = 40, detectBifurcation = 3, nInversion = 4, tolBisectionEigenvalue = 1e-17, dsminBisection = 1e-7)

function cb(state; kwargs...)
	fromNewton = get(kwargs, :fromNewton, false)
	_x = get(kwargs, :z0, nothing)
	if ~fromNewton
		return norm(_x.u - state.x) < 20.5 && abs(_x.p - state.p)<0.05
	end
	true
end

kwargsC = (verbosity = 3,
	plot = true,
	# tangentAlgo = BorderedPred(),
	linearAlgo  = MatrixBLS(),
	recordFromSolution = (x, p) -> (n2 = norm(x), nw = normweighted(x), s = sum(x), s2 = x[end ÷ 2], s4 = x[end ÷ 4], s5 = x[end ÷ 5]),
	plotSolution = (x, p;kwargs...)->(plot!(X, x; ylabel="solution", label="", kwargs...)),
	callbackN = cb
	)

brflat, = @time continuation(
	R_SH, Jac_sp, sol1, (@set parSH.l = -0.1), (@lens _.l), opts;
	kwargsC...)

plot(brflat, putspecialptlegend = false)
####################################################################################################
# branch switching
function optrec(x, p, l; opt = opts)
	level =  l
	if level <= 2
		return setproperties(opt; maxSteps = 300, detectBifurcation = 3, nev = Nx, detectLoop = false)
	else
		return setproperties(opt; maxSteps = 250, detectBifurcation = 3, nev = Nx, detectLoop = true)
	end
end

diagram = @time bifurcationdiagram(jet..., sol1, (@set parSH.l = -0.1), (@lens _.l), 2, optrec; kwargsC..., halfbranch = true, verbosity = 0, usedeflation = true)

code = ()
	vars = (:param, :n2)
	plot(diagram; code = code, plotfold = false,  markersize = 2, putspecialptlegend = false, vars = vars)
	plot!(brflat, putspecialptlegend = false, vars = vars)
	title!("#branches = $(size(diagram, code))")

diagram2 = bifurcationdiagram!(jet..., BK.getBranch(diagram, (2,)), 3, optrec; kwargsC..., halfbranch = true)

####################################################################################################
deflationOp = DeflationOperator(2, 1.0, [sol1])

br, = @time continuation(
	R_SH, Jac_sp,
	(@set parSH.l = -1.), (@lens _.l),
	setproperties(opts; ds = 0.001, maxSteps = 20000, pMax = 0.25, pMin = -1., newtonOptions = setproperties(optnew; tol = 1e-9, maxIter = 15, verbose = false), saveSolEveryStep = 0, detectBifurcation = 0),
	deflationOp;
	verbosity = 1,
	maxBranches = 150,
	perturbSolution = (sol, p, id) -> sol .+ 0.02 .* rand(length(sol)),
	normN = x -> norm(x, Inf64),
	recordFromSolution = (x, p) -> (s5 = x[end ÷ 5], n2 = norm(x), nw = normweighted(x), s = sum(x), s2 = x[end ÷ 2], s4 = x[end ÷ 4],),
	# tangentAlgo = SecantPred(),
	)

plot(br..., legend=false, linewidth=1, vars = (:param, :n2))


text:
#! /usr/bin/env julia
println(sum(x -> parse(Int, x), readlines()))

max_stars_repo_name: paolobueno/advent-of-julia

max_issues_repo_name: paolobueno/advent-of-julia

max_forks_repo_name: paolobueno/advent-of-julia


text:
module Dtree

using Base.Threads

export DtreeScheduler, nnodes, nodeid, initwork, getwork, runtree, sync, cpu_pause

const fan_out = 2048
const drain_rate = 0.4

const libdtree = joinpath(Pkg.dir("Dtree"), "deps", "Dtree",
        "libdtree.$(Libdl.dlext)")

enter_gc_safepoint() = ccall(:jl_gc_safe_enter, Int8, ())
leave_gc_safepoint(gs) = ccall(:jl_gc_safe_leave, Void, (Int8,), gs)

function __init__()
    ccall((:dtree_init, libdtree), Cint, (Cint, Ptr{Ptr{UInt8}}),
          length(ARGS), ARGS)
    global const nnodes =
	    Int(ccall((:dtree_nnodes, libdtree), Cint, ()))
    global const nodeid =
	    Int(ccall((:dtree_nodeid, libdtree), Cint, ())+1)
    atexit() do
        ccall((:dtree_shutdown, libdtree), Cint, ())
    end
end

type DtreeScheduler
    handle::Array{Ptr{Void}}

    function DtreeScheduler(fan_out::Int, num_work_items::Int64,
            can_parent::Bool, node_mul::Float64,
            first::Float64, rest::Float64, min_dist::Int)
        parents_work = nthreads()>1 ? 1 : 0
        cthrid = cfunction(threadid, Int64, ())
        d = new([0])
        p = [ 0 ]
        r = ccall((:dtree_create, libdtree), Cint, (Cint, Cint, Cint, Cint,
                  Cdouble, Cint, Ptr{Void}, Cdouble, Cdouble, Cshort,
                  Ptr{Void}, Ptr{Int64}), fan_out, num_work_items,
                  can_parent, parents_work, node_mul, nthreads(), cthrid,
                  first, rest, min_dist, pointer(d.handle), pointer(p, 1))
        if r != 0
            error("construction failure")
        end
        finalizer(d, (x -> ccall((:dtree_destroy, libdtree),
                                 Void, (Ptr{Void},), d.handle[1])))
        d, Bool(p[1])
    end
end

DtreeScheduler(num_work_items::Int64, first::Float64) =
    DtreeScheduler(fan_out, num_work_items, true, 1.0, first, drain_rate, 1)
DtreeScheduler(num_work_items::Int64, first::Float64, min_distrib::Int) =
    DtreeScheduler(fan_out, num_work_items, true, 1.0, first, drain_rate, min_distrib)

function initwork(dt::DtreeScheduler)
    w = [ 1, 1 ]::Array{Int64}
    wp1 = pointer(w, 1)
    wp2 = pointer(w, 2)
    gs = enter_gc_safepoint()
    r = ccall((:dtree_initwork, libdtree), Cint,
            (Ptr{Void}, Ptr{Int64}, Ptr{Int64}), dt.handle[1], wp1, wp2)
    leave_gc_safepoint(gs)
    return r, (w[1]+1, w[2])
end

function getwork(dt::DtreeScheduler)
    w = [ 1, 1 ]::Array{Int64}
    wp1 = pointer(w, 1)
    wp2 = pointer(w, 2)
    gs = enter_gc_safepoint()
    r = ccall((:dtree_getwork, libdtree), Cint,
            (Ptr{Void}, Ptr{Int64}, Ptr{Int64}), dt.handle[1], wp1, wp2)
    leave_gc_safepoint(gs)
    return r, (w[1]+1, w[2])
end

function runtree(dt::DtreeScheduler)
    r = 0
    gs = enter_gc_safepoint()
    r = ccall((:dtree_run, libdtree), Cint, (Ptr{Void},), dt.handle[1])
    leave_gc_safepoint(gs)
    Bool(r > 0)
end

@inline sync() = ccall((:dtree_sync, libdtree), Void, ())
@inline cpu_pause() = ccall((:cpu_pause, libdtree), Void, ())
@inline rdtsc() = ccall((:rdtsc, libdtree), Culonglong, ())

end # module


text:
using Unitful: unit

@testset "autoperiod self-consistency" begin
    t, y, dy, params = make_data()
    durations = params.duration .+ range(-0.1, 0.1, length=3)

    period = @inferred autoperiod(t, durations)
    results1 = @inferred BLS(t, y, dy; duration=durations, periods=period)
    results2 = @inferred BLS(t, y, dy; duration=durations)
    @test BoxLeastSquares.power(results1) ≈ BoxLeastSquares.power(results2)

end

@testset "model correctness ($obj)" for obj in [:likelihood, :snr]
    t, y, dy, params = make_data(1000)

    periods = exp.(range(log(params.period) - 0.1, log(params.period) + 0.1, length=1000))

    results = @inferred BLS(t, y, dy; params.duration, periods, objective=obj)
    best_params = BoxLeastSquares.params(results)

    @test best_params.period ≈ params.period atol = 0.01
    @test best_params.t0 ≈ params.t0 atol = 0.01
    @test best_params.duration ≈ params.duration atol = 0.01
    @test best_params.depth ≈ params.depth atol = best_params.depth_err
end

@testset "transit model" begin
    t, y, dy, params = make_data()

    # compute model using linear regression
    A = zeros(length(t), 2)
    dt = @. abs((t - params.t0 + 0.5 * params.period) % params.period - 0.5 * params.period)
    intransit = @. dt < 0.5 * params.duration
    A[.!intransit, 1] .= 1
    A[intransit, 2] .= 1
    w = (A' * (A ./ dy .^2)) \ (A' * (y ./ dy .^2))
    model_true = A * w

    model = @inferred BoxLeastSquares.model(t, y, dy; params...)
    @test model ≈ model_true

    # test interface
    result = BLS(t, y, dy; params.duration)
    model_direct = @inferred BoxLeastSquares.model(result; params...)
    @test model ≈ model_direct
end

@testset "Unitful.jl integration" begin
    @testset "time units" begin
        t, y, dy, params = make_data()

        tu = t * u"d"
        duration = 0.16u"d" |> u"hr"

        results = BLS(tu, y, dy; duration)
        @test unit(eltype(results.periods)) == u"d"
        @test unit(eltype(results.t0)) == u"d"
        @test unit(eltype(results.duration)) == u"hr"
        @test unit(eltype(results.duration_in)) == u"hr"
    end
    @testset "flux units" begin
        t, y, dy, params = make_data()

        yu = y * u"J"
        dyu = dy * u"J"

        results = BLS(t, yu, dyu; params.duration)
        @test unit(eltype(results.depth)) == u"J"
    end
    @testset "flux and time units" begin
        t, y, dy, params = make_data()

        tu = t * u"d"
        yu = y * u"J"
        dyu = dy * u"J"
        duration = 0.16u"d" |> u"hr"

        results = BLS(tu, yu, dyu; duration)
        @test unit(eltype(results.periods)) == u"d"
        @test unit(eltype(results.t0)) == u"d"
        @test unit(eltype(results.duration)) == u"hr"
        @test unit(eltype(results.duration_in)) == u"hr"
        @test unit(eltype(results.depth)) == u"J"
    end
end


text:
############################################################################################
"""
$(TYPEDEF)

Default arguments for printing summary statistics after sampling.

# Fields
$(TYPEDFIELDS)
"""
struct PrintDefault
    "Number of Digits for rounding summary statistics."
    Ndigits::Int64
    "Quantiles for summary statistics."
    quantiles::Vector{Float64}
    function PrintDefault(; Ndigits=3, quantiles=[0.025, 0.25, 0.50, 0.75, 0.975])
        ArgCheck.@argcheck 0 < Ndigits
        ArgCheck.@argcheck all(
            0.0 <= quantiles[iter] <= 1.0 for iter in eachindex(quantiles)
        ) "Quantiles should be between 0.0 and 1.0."
        return new(Ndigits, quantiles)
    end
end

############################################################################################
# Export
export PrintDefault


text:
# Autogenerated wrapper script for LCIO_jll for x86_64-apple-darwin
export liblcio, libsio

using Zlib_jll
JLLWrappers.@generate_wrapper_header("LCIO")
JLLWrappers.@declare_library_product(liblcio, "@rpath/liblcio.2.17.dylib")
JLLWrappers.@declare_library_product(libsio, "@rpath/libsio.0.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        liblcio,
        "lib/liblcio.2.17.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsio,
        "lib/libsio.0.0.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export Chebyshev

function Chebyshev(nlp_stop, PData::PDataFact, dₙ::Vector, gt::Vector)
    # printstyled("on est dans Chebyshev  \n", color = :yellow)
    T = eltype(gt)
    # @show T
    nlp_at_x = nlp_stop.current_state
    x = nlp_at_x.x
    # xdemi = NaN * rand(length(dₙ))
    # @show x
    # @show ∇f³xuv
    # @show typeof(nlp_stop.pb)
    Tdndn = ∇f³xuv(nlp_stop.pb, x, dₙ, dₙ)
    # printstyled("on a Tdndn  \n", color = :yellow)

    ϵ2 = sqrt(eps(eltype(gt)))
    Γ = max.(abs.(PData.Δ), ϵ2)

    d̃ = PData.L \ Tdndn[PData.pp]
    d̂ = PData.L' \ (PData.Q * (PData.Q' * d̃ ./ Γ))
    d = -d̂[invperm(PData.pp)]

    dC = dₙ + T(0.5) * d
    return dC
end


text:
using Test, DataValueInterfaces

@testset "DataValueInterfaces" begin

    @testset "nondatavaluetype" begin

        @test DataValueInterfaces.nondatavaluetype(Int64) == Int64
        @test DataValueInterfaces.nondatavaluetype(Union{}) == Union{}

    end

    @testset "datavaluetype" begin

        @test DataValueInterfaces.datavaluetype(Int64) == Int64
        @test DataValueInterfaces.datavaluetype(Union{}) == Union{}

    end

end # @testset "DataValueInterfaces"


text:
using Distributions
##################
#basically plumbing to generate a multinomial problem with the setup given below
#predictor matrix, p>n
num=200
n=100
X=zeros(n,num)
X[:]=rand(n*num)

#intercept
X[:,1]=1

#true beta values, first 10 rows are nonzero
#four different response categories
beta=zeros(num,4)
beta[1:10,:]=(rand(4*10)-1/2)*3
Xbeta=X*beta
probs = BlockCoordinateDescent.multProbPrecompute(Xbeta)

#simulate response
z=zeros(n,4)
for i=1:n
  t1=Multinomial(1,vec(probs[i,:]))
  z[i,:]=rand(t1,1)
end

#define grouping: columns with same group number/label are considered grouped
groups=zeros(num)
for i=1:(div(num,2)-1)
  groups[2*i+1]=i
  groups[2*i+2]=i
end

groups[1]=-1
groups[2]=0

#the intercept is unpenalized, every other group is penalized equally
penaltyFactor=ones(length(unique(groups)))
penaltyFactor[1]=0
#####################

@time multinomialBCD(X,z,groups,penaltyFactor);


text:
using LightXML
using Dates

mutable struct Params
  lattice_file::String
  hoppings::String # nearest neighbor
  Nhoppings::String # next nearest neighbor
  NNhoppings::String # next next nearest neighbor
  beta::Float64
  delta_tau::Float64
  L::Int
  slices::Int
  safe_mult::Int
  opdim::Int # order parameter dimension

  thermalization::Int # no measurements, no saving
  measurements::Int # save and (maybe) measure

  hsfield::Array{Float64, 3} # dim 1: op component, dim 2: linearized spatial lattice, dim 3: imag time
  boson_action::Float64

  mu1::Float64 # x flavor
  mu2::Float64 # y flavor
  lambda::Float64
  r::Float64
  c::Float64
  u::Float64
  flv::Int # flavors: GF matl_sites size flv*l.sites x flv*l.sites

  box::Float64
  box_global::Float64

  global_updates::Bool
  global_rate::Int64

  sparsity_limit::Float64 # upper sparsity limit for chkr folding heuristic
  chkr::Bool # checkerboard
  Bfield::Bool # artificial magnetic field to reduce finite size effects

  output_file::String
  write_every_nth::Int # measure after 'write_every_nth' up-down sweeps(!)

  all_checks::Bool # true: check for propagation instblts. and S_b consistency

  seed::Int

  resume::Bool
  prethermalized::Int # how many thermal ud-sweeps are captured by custom start configuration 'thermal_init/conf'

  edrun::Bool # if true, only the mass term of the bosonic action is considered
  walltimelimit::Dates.DateTime

  obs::Set{Symbol} # what observables to measure during DQMC (excl. configurations)

  function Params()
    p = new()
    p.global_updates = true
    p.chkr = true
    p.Bfield = false
    p.box = 0.5
    p.box_global = p.box
    p.global_rate = 5
    p.write_every_nth = 1
    p.all_checks = true
    p.opdim = 3
    p.flv = 4
    p.seed = 4729339882041979125
    p.resume = false
    p.prethermalized = 0
    p.sparsity_limit = 0.01
    p.hoppings = "none"
    p.Nhoppings = "none"
    p.NNhoppings = "none"
    p.edrun = false
    p.walltimelimit = Dates.DateTime("2099", "YYYY") # effective infinity
    # p.obs = Set{Symbol}()
    p.obs = Set{Symbol}()
    return p
  end
end

"""
    set_parameters(p::Params, params::Dict)

params::Dict -> p::Params (mandatory and optional).
Also triggers `deduce_remaining_parameters()` to e.g. calculate p.beta and set global DataTypes.
"""
function set_parameters(p::Params, params::Dict)
  ### PARSE MANDATORY PARAMS
  if haskey(params, "WARMUP")
    p.thermalization = parse(Int, params["WARMUP"])/2
  else
    p.thermalization = parse(Int, params["THERMALIZATION"])
  end
  if haskey(params, "SWEEPS")
    p.measurements = parse(Int, params["SWEEPS"])/2
  else
    p.measurements = parse(Int, params["MEASUREMENTS"])
  end
  p.delta_tau = parse(Float64, params["DELTA_TAU"])
  big_delta_tau = parse(BigFloat, params["DELTA_TAU"])
  if haskey(params, "BETA")
    p.slices = Int(parse(BigFloat, params["BETA"])/big_delta_tau)
    println("SLICES = ", p.slices)
  elseif haskey(params, "T")
    p.slices = Int((big"1"/parse(BigFloat, params["T"]))/big_delta_tau)
    println("SLICES = ", p.slices)
  else
    p.slices = parse(Int, params["SLICES"])
  end
  p.safe_mult = parse(Int, params["SAFE_MULT"])
  p.lattice_file = ""
  if haskey(params, "LATTICE_FILE")
    p.lattice_file = params["LATTICE_FILE"]
  else
    p.L = parse(Int, params["L"])
  end

  if haskey(params, "HOPPINGS")
    s = params["HOPPINGS"]
    sum(abs.(parse.(Float64, split(s, ",")))) > 0 && (p.hoppings = params["HOPPINGS"])
  end
  if haskey(params, "N-HOPPINGS")
    s = params["N-HOPPINGS"]
    sum(abs.(parse.(Float64, split(s, ",")))) > 0 && (p.Nhoppings = params["N-HOPPINGS"])
  end
  if haskey(params, "NN-HOPPINGS")
    s = params["NN-HOPPINGS"]
    sum(abs.(parse.(Float64, split(s, ",")))) > 0 && (p.NNhoppings = params["NN-HOPPINGS"])
  end
  if haskey(params, "MU1") && haskey(params, "MU2")
    p.mu1 = parse(Float64, params["MU1"])
    p.mu2 = parse(Float64, params["MU2"])
  elseif haskey(params, "MU1")
    p.mu1 = parse(Float64, params["MU1"])
    p.mu2 = parse(Float64, params["MU1"])
  else
    p.mu1 = parse(Float64, params["MU"])
    p.mu2 = parse(Float64, params["MU"])
  end

  p.lambda = parse(Float64, params["LAMBDA"])
  p.r = parse(Float64, params["R"])
  p.c = parse(Float64, params["C"])
  p.u = parse(Float64, params["U"])
  ###

  ### PARSE OPTIONAL PARAMS
  if haskey(params, "OPDIM")
    p.opdim = parse(Int, params["OPDIM"]);
    if p.opdim == 1 || p.opdim == 2
      p.flv = 2
    elseif p.opdim == 3
      p.flv = 4
    else
      error("OPDIM must be 1, 2 or 3!")
    end
  end
  haskey(params, "SEED") && (p.seed = parse(Int, params["SEED"]))
  haskey(params,"GLOBAL_UPDATES") && (p.global_updates = parse(Bool, lowercase(params["GLOBAL_UPDATES"])))
  haskey(params,"CHECKERBOARD") && (p.chkr = parse(Bool, lowercase(params["CHECKERBOARD"])))
  haskey(params,"BFIELD") && (p.Bfield = parse(Bool, lowercase(params["BFIELD"])))
  if haskey(params,"BOX_HALF_LENGTH")
    p.box = parse(Float64, params["BOX_HALF_LENGTH"])
    p.box_global = p.box
  end
  if haskey(params,"BOX_GLOBAL_HALF_LENGTH")
    p.box_global = parse(Float64, params["BOX_GLOBAL_HALF_LENGTH"])
  end
  haskey(params,"GLOBAL_RATE") && (p.global_rate = parse(Int64, params["GLOBAL_RATE"]))
  haskey(params,"WRITE_EVERY_NTH") && (p.write_every_nth = parse(Int64, params["WRITE_EVERY_NTH"]))
  haskey(params, "SPARSITY_LIMIT") && (p.sparsity_limit = parse(Float64, params["SPARSITY_LIMIT"]))
  haskey(params, "EDRUN") && (p.edrun = parse(Bool, lowercase(params["EDRUN"])))

  # observables
  for obs in (:greens, :boson_action)
    obs_str = uppercase(string(obs))
    if haskey(params, obs_str)
      if parse(Bool, lowercase(params[obs_str]))
        push!(p.obs, obs)
      elseif obs in p.obs
        pop!(p.obs, obs)
      end
    end
  end

  deduce_remaining_parameters(p)

  nothing
end

"""
    deduce_remaining_parameters(p::Params)

Assumes that `p` has been loaded from XML or HDF5 and sets remaining (dependent) fields in `p`.
"""
function deduce_remaining_parameters(p::Params)
  p.beta = p.slices * p.delta_tau
  p.hsfield = zeros(p.opdim, 1, p.slices) # just to initialize it somehow

  if p.lattice_file != ""
    p.L = parse(Int, p.lattice_file[findlast(isequal('_'), collect(p.lattice_file))+1:end-4])
  else
    hn = lowercase(gethostname())
    lat = ""
    if p.Nhoppings == "none" && p.NNhoppings == "none"
      lat = "square_L_$(p.L)_W_$(p.L).xml"
    else
      lat = "NNsquare_L_$(p.L)_W_$(p.L).xml"
    end

    if "LATTICES" in keys(ENV)
      p.lattice_file = joinpath(ENV["LATTICES"], lat)
    else
      # try to deduce lattice path from hostname
      if occursin("cheops", hn)
        p.lattice_file = "/projects/ag-trebst/bauer/lattices/"*lat
      elseif occursin("fz-juelich", hn) || occursin("juwels", hn)
        p.lattice_file = "/gpfs/homea/hku27/hku273/lattices/"*lat
      elseif occursin("thp", hn)
        p.lattice_file = "/home/bauer/lattices/"*lat
      elseif occursin("thinkable", hn)
        p.lattice_file = "C:/Users/carsten/Desktop/sciebo/lattices/"*lat
      elseif occursin("travis", hn)
        p.lattice_file = joinpath(dirname(dirname(@__FILE__)), "test/lattices/", lat)
      else
        @warn "Unrecognized host. Can't deduce lattice file path. Trying to use lattices in test/lattices. (likely to fail)"
        p.lattice_file = joinpath(dirname(dirname(@__FILE__)), "test/lattices/", lat)
      end
    end
  end
  nothing
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
module Sequences

using Interpolations
using Statistics

# export
#     MeasurementSequence,
#     getOverlap,
#     ssd,
#     displace

struct MeasurementSequence{A1 <: AbstractArray, A2 <: AbstractArray, T}
    x::A1
    y::A2
    _interpolate::T
end

function MeasurementSequence(x::AbstractVector{<:Real},y::AbstractVector{<:Real})
    x_ = x |> deepcopy
    y_ = y |> deepcopy
    # flip both x and y if x is decreasing
    # @show Statistics.median(diff(x))
    if Statistics.median(diff(x)) < 0
        reverse!(x_)
        reverse!(y_)
    end

    MeasurementSequence(x, y, LinearInterpolation(x_, y_))
end

(ms::MeasurementSequence)(xval) = ms._interpolate(xval)

function getOverlap(s1::MeasurementSequence, s2::MeasurementSequence)
    (max(minimum(s1.x),minimum(s2.x)), min(maximum(s1.x),maximum(s2.x)))
end

# FIXME BUT NOT ON FIRE (EMBERS), make `res` something on order of incoming sequence grids (and for parent)
function ssd(s1::MeasurementSequence, s2::MeasurementSequence;res::Float64=0.1)
    cr = getOverlap(s1, s2)
    # correlation will be calculated on new grid specific to each overlap
    x = cr[1]:res:cr[2]
    # interpolate both scalar sequences to new grid

    # @info length(x) any(isnan.(x)) any(isnan.(s1.(x))) any(isnan.(s2.(x))) ((1.0/length(x)) * sum((s1.(x) - s2.(x)).^2))
    return (1.0/length(x)) * sum((s1.(x) - s2.(x)).^2)
end

function displace(s::MeasurementSequence, d::Float64)
    MeasurementSequence(s.x .+d, s.y)
end

# _ssdCorr(s1,s2,x) = ssd(s1,displace(s2,x))

# southerly
# a_o = 20
# b_o = 15

# overlap_gt = (15,10) # relative
# overlap_gt = (10,15) # absolute

end


text:
using Test
using QuantLib

mesher1 = Concentrating1dMesher(0.1, 2.5, 1001)
@test mesher1.locations[1] == range(0.1, stop=2.5, length=1001)[1]
mesher2 = Concentrating1dMesher(0.1, 2.5, 1001, (1.1, 0.25), true)
@test mesher2.locations[2] ≈ 0.10333775170976389 atol=1.0e-8


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Base.SecretBuffer()

An [`IOBuffer`](@ref)-like object where the contents will be securely wiped when garbage collected.

It is considered best practice to wipe the buffer using `Base.shred!(::SecretBuffer)` as
soon as the secure data are no longer required. When initializing with existing data, the
`SecretBuffer!` method is highly recommended to securely zero the passed argument. Avoid
initializing with and converting to `String`s as they are unable to be securely zeroed.

# Examples
```jldoctest
julia> s = Base.SecretBuffer()
SecretBuffer("*******")

julia> write(s, 's', 'e', 'c', 'r', 'e', 't')
6

julia> seek(s, 0); Char(read(s, UInt8))
's': ASCII/Unicode U+0073 (category Ll: Letter, lowercase)

julia> Base.shred!(s)
SecretBuffer("*******")

julia> eof(s)
true
```
"""
mutable struct SecretBuffer <: IO
    data::Vector{UInt8}
    size::Int
    ptr::Int

    function SecretBuffer(; sizehint=128)
        s = new(Vector{UInt8}(undef, sizehint), 0, 1)
        finalizer(final_shred!, s)
        return s
    end
end

"""
    SecretBuffer(str::AbstractString)

A convenience constructor to initialize a `SecretBuffer` from a non-secret string.

Strings are bad at keeping secrets because they are unable to be securely
zeroed or destroyed. Therefore, avoid using this constructor with secret data.
Instead of starting with a string, either construct the `SecretBuffer`
incrementally with `SecretBuffer()` and [`write`](@ref), or use a `Vector{UInt8}` with
the `Base.SecretBuffer!(::Vector{UInt8})` constructor.
"""
SecretBuffer(str::AbstractString) = SecretBuffer(String(str))
function SecretBuffer(str::String)
    buf = codeunits(str)
    s = SecretBuffer(sizehint=length(buf))
    for c in buf
        write(s, c)
    end
    seek(s, 0)
    s
end
convert(::Type{SecretBuffer}, s::AbstractString) = SecretBuffer(String(s))

"""
    SecretBuffer!(data::Vector{UInt8})

Initialize a new `SecretBuffer` from `data`, securely zeroing `data` afterwards.
"""
function SecretBuffer!(d::Vector{UInt8})
    len = length(d)
    s = SecretBuffer(sizehint=len)
    for i in 1:len
        write(s, d[i])
    end
    seek(s, 0)
    securezero!(d)
    s
end

function unsafe_SecretBuffer!(s::Cstring)
    if s == C_NULL
        throw(ArgumentError("cannot convert NULL to SecretBuffer"))
    end
    len = Int(ccall(:strlen, Csize_t, (Cstring,), s))
    unsafe_SecretBuffer!(convert(Ptr{UInt8}, s), len)
end
function unsafe_SecretBuffer!(p::Ptr{UInt8}, len=1)
    if p == C_NULL
        throw(ArgumentError("cannot convert NULL to SecretBuffer"))
    end
    s = SecretBuffer(sizehint=len)
    for i in 1:len
        write(s, unsafe_load(p, i))
    end
    seek(s, 0)
    unsafe_securezero!(p, len)
    s
end

show(io::IO, s::SecretBuffer) = print(io, "SecretBuffer(\"*******\")")

# Unlike other IO objects, equality is computed by value for convenience
==(s1::SecretBuffer, s2::SecretBuffer) = (s1.ptr == s2.ptr) && (s1.size == s2.size) && (UInt8(0) == _bufcmp(s1.data, s2.data, min(s1.size, s2.size)))
# Also attempt a constant time buffer comparison algorithm — the length of the secret might be
# inferred by a timing attack, but not its values.
@noinline function _bufcmp(data1::Vector{UInt8}, data2::Vector{UInt8}, sz::Int)
    res = UInt8(0)
    for i = 1:sz
        res |= xor(data1[i], data2[i])
    end
    return res
end
# All SecretBuffers hash the same to avoid leaking information or breaking consistency with ==
const _sb_hash = UInt === UInt32 ? 0x111c0925 : 0xb06061e370557428
hash(s::SecretBuffer, h::UInt) = hash(_sb_hash, h)

function write(io::SecretBuffer, b::UInt8)
    if io.ptr > length(io.data)
        # We need to resize! the array: do this manually to ensure no copies are left behind
        newdata = Vector{UInt8}(undef, (io.size+16)*2)
        copyto!(newdata, io.data)
        securezero!(io.data)
        io.data = newdata
    end
    io.size == io.ptr-1 && (io.size += 1)
    io.data[io.ptr] = b
    io.ptr += 1
    return 1
end

function write(io::IO, s::SecretBuffer)
    nb = 0
    for i in 1:s.size
        nb += write(io, s.data[i])
    end
    return nb
end

cconvert(::Type{Cstring}, s::SecretBuffer) = unsafe_convert(Cstring, s)
function unsafe_convert(::Type{Cstring}, s::SecretBuffer)
    # Ensure that no nuls appear in the valid region
    if any(==(0x00), s.data[i] for i in 1:s.size)
        throw(ArgumentError("`SecretBuffers` containing nul bytes cannot be converted to C strings"))
    end
    # Add a hidden nul byte just past the end of the valid region
    p = s.ptr
    s.ptr = s.size + 1
    write(s, '\0')
    s.ptr = p
    s.size -= 1
    return Cstring(unsafe_convert(Ptr{Cchar}, s.data))
end

seek(io::SecretBuffer, n::Integer) = (io.ptr = max(min(n+1, io.size+1), 1); io)
seekend(io::SecretBuffer) = seek(io, io.size+1)
skip(io::SecretBuffer, n::Integer) = seek(io, position(io) + n)

bytesavailable(io::SecretBuffer) = io.size - io.ptr + 1
position(io::SecretBuffer) = io.ptr-1
eof(io::SecretBuffer) = io.ptr > io.size
isempty(io::SecretBuffer) = io.size == 0
function peek(io::SecretBuffer, ::Type{UInt8})
    eof(io) && throw(EOFError())
    return io.data[io.ptr]
end
function read(io::SecretBuffer, ::Type{UInt8})
    eof(io) && throw(EOFError())
    byte = io.data[io.ptr]
    io.ptr += 1
    return byte
end

function final_shred!(s::SecretBuffer)
    !isshredded(s) && @async @warn("a SecretBuffer was `shred!`ed by the GC; use `shred!` manually after use to minimize exposure.")
    shred!(s)
end

function shred!(s::SecretBuffer)
    securezero!(s.data)
    s.ptr = 1
    s.size = 0
    return s
end

isshredded(s::SecretBuffer) = all(iszero, s.data)

function shred!(f::Function, x)
    try
        f(x)
    finally
        shred!(x)
    end
end

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: TechPenguineer/julia


text:
@inline function un_kernel(nbit::Int, cbits::NTuple{C, Int}, cvals::NTuple{C, Int}, U0::AbstractMatrix, locs::NTuple{M, Int}) where {C, M}
    U = (all(TupleTools.diff(locs).>0) ? U0 : reorder(U0, collect(locs)|>sortperm)) |> staticize
    MM = size(U0, 1)
    locked_bits = [cbits..., locs...]
    locked_vals = [cvals..., zeros(Int, M)...]
    locs_raw = [i+1 for i in itercontrol(nbit, setdiff(1:nbit, locs), zeros(Int, nbit-M))] |> staticize
    configs = itercontrol(nbit, locked_bits, locked_vals)

    length(configs), @inline function kernel(state, inds)
        x = @inbounds configs[inds[1]]
        unrows!(piecewise(state, inds), x .+ locs_raw, U)
    end
end

@inline function u1_kernel(nbit::Int, U1::AbstractMatrix, ibit::Int)
    a, c, b, d = U1
    step = 1<<(ibit-1)
    configs = itercontrol(nbit, [ibit], [0])

    length(configs), @inline function kernel(state, inds)
        i = @inbounds configs[inds[1]]+1
        u1rows!(piecewise(state, inds), i, i+step, a, b, c, d)
    end
end
u1_kernel(nbit::Int, U1::SDSparseMatrixCSC, ibit::Int) = u1_kernel(nbit, U1|>Matrix, ibit)

@inline function u1_kernel(nbit::Int, U1::SDPermMatrix, ibit::Int)
    U1.perm[1] == 1 && return u1_kernel(nbit, Diagonal(U1.vals), ibit)

    mask = bmask(ibit)
    b, c = U1.vals[1], U1.vals[2]
    step = 1<<(ibit-1)
    configs = itercontrol(nbit, [ibit], [0])

    length(configs), @inline function kernel(state, inds)
        x = @inbounds configs[inds[1]] + 1
        swaprows!(piecewise(state, inds), x, x+step, c, b)
    end
end

@inline function u1_kernel(nbit::Int, U1::SDDiagonal, ibit::Int)
    a, d = U1.diag
    u1diag_kernel(nbit,ibit, a, d)
end

@inline function u1diag_kernel(nbit::Int, ibit::Int, a, d)
    mask = bmask(ibit)
    1<<nbit, @inline function kernel(state, inds)
        i = inds[1]
        piecewise(state, inds)[i] *= anyone(i-1, mask) ? d : a
    end
end

################ Specific kernels ##################
x_kernel(nbit::Int, bits::Ints) = cx_kernel(nbit::Int, (), (), bits::Ints)
cx_kernel(nbit::Int, cbits, cvals, loc::Int) = cx_kernel(nbit::Int, cbits, cvals, (loc,))
@inline function cx_kernel(nbit::Int, cbits, cvals, bits::Ints)
    configs = itercontrol(nbit, [cbits..., bits[1]], [cvals..., 0])
    mask = bmask(bits...)
    length(configs), @inline function kernel(state, inds)
        b = @inbounds configs[inds[1]]
        swaprows!(piecewise(state, inds), b+1, flip(b, mask) + 1)
        return
    end
end

@inline function pswap_kernel(nbit::Int, m::Int, n::Int, theta::Real)
    mask1 = bmask(m)
    mask2 = bmask(n)
    mask12 = mask1|mask2
    a, c, b_, d = mat(Rx(theta))
    e = exp(-im/2*theta)
    configs = itercontrol(nbit, [m,n], [0,0])
    length(configs), @inline function kernel(state, inds)
        @inbounds x = configs[inds[1]]
        piecewise(state, inds)[x+1] *= e
        piecewise(state, inds)[x⊻mask12+1] *= e
        y = x ⊻ mask2
        u1rows!(piecewise(state, inds), y+1, y⊻mask12+1, a, b_, c, d)
    end
end

@inline function y_kernel(nbit::Int, bits::Ints)
    mask = bmask(Int, bits...)
    configs = itercontrol(nbit, [bits[1]], [0])
    bit_parity = length(bits)%2 == 0 ? 1 : -1
    factor = (-im)^length(bits)
    length(configs), @inline function kernel(state, inds)
        b = @inbounds configs[inds[1]]
        i_ = flip(b, mask) + 1
        factor1 = count_ones(b&mask)%2 == 1 ? -factor : factor
        factor2 = factor1*bit_parity
        swaprows!(piecewise(state, inds), b+1, i_, factor2, factor1)
    end
end

@inline function z_kernel(nbit::Int,bits::Ints)
    mask = bmask(bits...)
    1<<nbit,@inline function kernel(state, inds)
        i = inds[1]
        piecewise(state, inds)[i] *= count_ones((i-1)&mask)%2==1 ? -1 : 1
        return
    end
end

@inline function zlike_kernel(nbit::Int,bits::Ints, d::Union{ComplexF32, ComplexF64, Float64, Float32})
    mask = bmask(Int32, bits...)
    1<<nbit,@inline function kernel(state, inds)
        i = inds[1]
        piecewise(state, inds)[i] *= CUDAnative.pow(d, bit_count(Int32(i-1)&mask))
        return
    end
end

@inline function cdg_kernel(nbit::Int, cbits, cvals, ibit, a, d)
    ctrl = controller((cbits..., ibit), (cvals..., 1))

    1<<nbit, @inline function kernel(state, inds)
        i = inds[1]
        piecewise(state, inds)[i] *= ctrl(i-1) ? d : a
    end
end

@inline function cy_kernel(nbit, cbits, cvals, bit::Int)
    configs = itercontrol(nbit, [cbits..., bit], [cvals..., 0])
    mask = bmask(bit)
    length(configs), @inline function kernel(state, inds)
        b = @inbounds configs[inds[1]]
        swaprows!(piecewise(state, inds), b+1, flip(b, mask) + 1, im, -im)
    end
end

for (G, FACTOR) in zip([:z, :s, :t, :sdag, :tdag], [:(-one(Int32)), :(1f0im), :($(exp(im*π/4))), :(-1f0im), :($(exp(-im*π/4)))])
    KERNEL = Symbol(G, :_kernel)
    CKERNEL = Symbol(:c, KERNEL)
    if G != :z
        @eval $KERNEL(nbit::Int, bits::Ints) = zlike_kernel(nbit, bits, $FACTOR)
    end
    @eval $KERNEL(nbit::Int, bit::Int) = u1diag_kernel(nbit::Int, bit, one($FACTOR), $FACTOR)
    @eval $CKERNEL(nbit::Int, cbits, cvals, ibit::Int) = cdg_kernel(nbit::Int, cbits, cvals, ibit, one($FACTOR), $FACTOR)
end

"""
build a simple kernel function from functions like f(::VecOrMat, inds).
"""
@generated function simple_kernel(func, state::AbstractArray{T, N}) where {T, N}
    if N == 1
        ex = :(inds = (blockIdx().x-1) * blockDim().x + threadIdx().x)
    else
        ex = :(inds = ((blockIdx().x-1) * blockDim().x + threadIdx().x,
                       (blockIdx().y-1) * blockDim().y + threadIdx().y); inds[2]<=size(state, 2) || return nothing)
    end
    ex = :($ex; func(state, inds); nothing)
    ex
end


text:
module SWScale
  include(joinpath(dirname(@__FILE__),"..","..","..","init.jl"))
  w(f) = joinpath(swscale_dir, f)

  using AVUtil

  include(w("LIBSWSCALE.jl"))
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO