text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -13.266290625856113 <= q <= 11.094199160852137)

                     add_NL_constraint(m, :(gelu(-0.2977454197977769 + 0.7722469146312885*gelu(0.9100824262656433 + 0.0704650714228956*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.8966302258455006*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.7890296052704904*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.5815923034237733*gelu(-0.15598180224250147 + 0.7627892034632455*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.6806349200893864*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.163059274878961*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + -0.731486602090937*gelu(-0.9577813343303001 + 0.660774759467007*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.16269598559087406*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.9308521240523335*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + gelu(0.23088458854123672 + 0.9025271253252471*gelu(0.9100824262656433 + 0.0704650714228956*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.8966302258455006*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.7890296052704904*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.9612479731086991*gelu(-0.15598180224250147 + 0.7627892034632455*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.6806349200893864*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.163059274878961*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.9633061422344626*gelu(-0.9577813343303001 + 0.660774759467007*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.16269598559087406*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.9308521240523335*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) + gelu(-0.8299198997441475 + 0.7709117851978275*gelu(0.9100824262656433 + 0.0704650714228956*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.8966302258455006*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.7890296052704904*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.5338511667423695*gelu(-0.15598180224250147 + 0.7627892034632455*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.6806349200893864*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.163059274878961*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4])))) + 0.3695524647159538*gelu(-0.9577813343303001 + 0.660774759467007*gelu(-0.7110372864748489 + -0.3189000958940382*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.039301419137926175*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.19790129907631693*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + -0.16269598559087406*gelu(0.7299768289984225 + 0.7045346620386517*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.8068588624409387*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + 0.8888276792068996*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))) + 0.9308521240523335*gelu(-0.5099367345236407 + -0.3810722083004521*gelu(0.8606897861111369 + 0.08970944375635792*$(x[1]) + 0.5451334799438907*$(x[2]) + -0.8811882246809297*$(x[3]) + -0.4223772070829366*$(x[4])) + 0.4850891682061631*gelu(0.47662806300654736 + 0.5011612847792533*$(x[1]) + 0.5521868073085425*$(x[2]) + 0.45455279409713745*$(x[3]) + -0.33071846313642217*$(x[4])) + -0.158929561309892*gelu(-0.45543581836352764 + -0.399986773274851*$(x[1]) + -0.07747330631488092*$(x[2]) + -0.8376058727182878*$(x[3]) + 0.18120897903061906*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Fermi
using LinearAlgebra

Nt = Threads.nthreads()
Fermi.tblis_set_num_threads(Nt)
BLAS.set_num_threads(Nt)
output("NUMBER OF THREADS: {}", Threads.nthreads())

function get_hc(N)

    molstring = ""
    f = false

    for l = eachline(joinpath(@__DIR__, "alkanes.xyz"))
        if occursin("C$N", l)
            f = true
            continue
        end

        if f && occursin("--", l)
            break
        end

        if f
            molstring *= l*"\n"
        end
    end

    Fermi.Options.set("molstring", molstring)
end

MAX_CHAIN_LENGTH = 10

# Cold run
@set printstyle none
@energy CCSD

@reset
@set {
    basis cc-pvdz
    printstyle file
}

for N = 1:MAX_CHAIN_LENGTH

    get_hc(N)
    fc = 2*N
    Fermi.Options.set("drop_occ", fc)

    run(`echo TBLIS N = $N Threads = $Nt`)
    @set tblis true
    Fermi.Options.set("output", "TBLIS_C"*"$N"*"_T"*"$Nt"*".out")
    output("NUMBER OF THREADS: {}", Threads.nthreads())
    @energy ccsd

    run(`echo BLAS N = $N Threads = $Nt`)
    @set tblis false
    output("NUMBER OF THREADS: {}", Threads.nthreads())
    Fermi.Options.set("output", "BLAS_C"*"$N"*"_T"*"$Nt"*".out")
    @energy ccsd
end

ext: jl


text:
module TraitDispatch

using InteractiveUtils
using Printf
using MacroTools
using MacroTools: postwalk, @q, combinedef

export
    AbstractTrait,
    NullTrait,
    @define_trait,
    @implement_trait,
    @define_traitfn,
    @implement_traitfn,
    istrait,
    isleaftrait,
    subtraits,
    hastrait,
    hastraits,
    implementstrait,
    implementstraits,
    notraits,
    traits

abstract type AbstractTrait{T} end
abstract type NullTrait{T} end
traitdispatch(::Type{Trait}) where {T,Trait<:AbstractTrait{T}} = NullTrait{T}

include("trait_functions.jl")
include("trait_definition.jl")
include("trait_utils.jl")

end # module


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

if !isdefined(@__MODULE__, :DefsChannelName)
    const DefsChannelName = String
else
    @warn("Skipping redefinition of DefsChannelName to String")
end


text:
using KnowledgeGraphNote
using Test
include("notes.jl")

@testset "KnowledgeGraphNote.jl" begin

    ######################### testing library api
    # print_concepts(math_notes)
    @test Set(get_missing_concepts(math_notes)) == Set(["affinely independent", "measure", "limit point", "empty set", "power set", "neighborhood", "set closure", "compact set"])
    @test get_duplicate_concepts(duplicate_notes) == Dict("a"=>2)
    @test length(get_duplicate_concepts(math_notes)) == 0
    @test Set(get_known_concepts(duplicate_notes)) == Set(["a", "b"])

    mathkg = init_knowledge_graph(math_notes)
    @test length(find_cycles(mathkg)) == 0
    cyclickg = init_knowledge_graph(cyclic_notes)
    @test [["a", "b", "c"]] == find_cycles(cyclickg)
    # println(generate_learning_plan(mathkg, math_notes,  "Borel measure"))

    target = "borel sigma-algebra"
    @test mathkg.idtoname[generate_learning_plan(mathkg, target)] == ["set", "sigma-algebra", "metric", "metric space", "open set", "closed set", "topological space", "borel sigma-algebra"]

    #### the dot file format is described here: https://en.wikipedia.org/wiki/DOT_(graph_description_language). you can find a viewer from there.
    #### it's not a good idea to write files in unit tests. uncommment these lines locally to see what they do.
    # export_knowledge_graph(mathkg, "./math.dot")
    # export_knowledge_graph_towards_target(mathkg, target, "./target.dot")
    # write_notes_for_target_concept_to_md_file(mathkg, math_notes, target, "./target.md")

    @test check_notes(duplicate_notes) == false
    @test check_notes(math_notes) == false
    @test check_notes(cyclic_notes) == false
    @test getfield.(fuzzy_search(math_notes, "rado measure"), :name) == ["radon measure"]

    ######################### testing library internal functions
    dfs_postordering_test1_kg = KnowledgeGraphNote.init_knowledge_graph(dfs_postordering_test1)
    # KnowledgeGraphNote.export_knowledge_graph(dfs_postordering_test1_kg, "./dfs_postordering_test1.dot")
    postorder = KnowledgeGraphNote.dfs_postordering(dfs_postordering_test1_kg.graph, 1)
    @test ["d", "b", "e", "c", "t", "a"] == dfs_postordering_test1_kg.idtoname[postorder]

    dfs_postordering_test2_kg = KnowledgeGraphNote.init_knowledge_graph(dfs_postordering_test2)
    # KnowledgeGraphNote.export_knowledge_graph(dfs_postordering_test2_kg, "./dfs_postordering_test2.dot")
    postorder = KnowledgeGraphNote.dfs_postordering(dfs_postordering_test2_kg.graph, 1)
    @test ["d", "b", "c", "a"] == dfs_postordering_test2_kg.idtoname[postorder]

    dfs_postordering_test3_kg = KnowledgeGraphNote.init_knowledge_graph(dfs_postordering_test3)
    # KnowledgeGraphNote.export_knowledge_graph(dfs_postordering_test3_kg, "./dfs_postordering_test3.dot")
    postorder = KnowledgeGraphNote.dfs_postordering(dfs_postordering_test3_kg.graph, 1)
    @test ["f", "e", "b", "d", "c", "a"] == dfs_postordering_test3_kg.idtoname[postorder]

end

max_stars_repo_name: bluepine/KnowledgeGraphNote

max_issues_repo_name: bluepine/KnowledgeGraphNote

max_forks_repo_name: bluepine/KnowledgeGraphNote


text:
using MarketData

facts("Volatility") do

    context("bollinger_bands") do
        @fact bollingerbands(cl)["up"].values[1]   --> roughly(117.3251, atol=.01)      # TTR default uses sample=FALSE and value 117.3251
        @fact bollingerbands(cl)["down"].values[1] --> roughly(89.39392, atol=.01)      # TTR default uses sample=FALSE and value 89.39392
        @fact bollingerbands(cl)["mean"].values[1] --> roughly(103.3595, atol=.01)      # TTR 103.3595
        @fact bollingerbands(cl).timestamp[1]      --> Date(2000,1,31)
        @fact bollingerbands(cl).timestamp[end]    --> Date(2001,12,31)
    end

    context("truerange") do
        @fact truerange(ohlc).values[end]    --> roughly(0.83, atol=.01)       # TTR 0.83
        @fact truerange(ohlc).timestamp[end] --> Date(2001,12,31)
    end

    context("atr") do
        @fact atr(ohlc).values[1]    --> roughly(8.343571428571428)   # TTR value 8.343571
        @fact atr(ohlc).values[end]  --> roughly(0.9664561242651976)  # TTR value 0.9664561
        @fact atr(ohlc).timestamp[1] --> Date(2000,1,24)
    end

   #  context("keltner_bands") do
   #      @fact keltnerbands(ohlc)["kma"].values[end] --> roughly(21.42)  # needs confirmation
   #      @fact keltnerbandb(ohlc)["kup"].values[end] --> roughly(22.32)  # needs confirmation
   #      @fact keltnerbands(ohlc)["kdn"].values[end] --> roughly(20.52)  # needs confirmation
   #      @fact keltnerbands(ohlc).timestamp[1]       --> Date(2000,1,14)
   #  end

   #  context("chaikin_volatility") do
   #     @fact chk --> 17.0466
   #  end
end


text:
seed!(12345)

function testMultinomial(n::Int64, p::Vector{Float64})
  counts::Vector{Int64} = sampleMultinomial(n, p)
  return testGOFMultinomial(p, counts)
end

@test testMultinomial(100, [0.5; 0.1; 0.25; 0.15]) > 0.01
@test testMultinomial(10000, [0.5; 0.1; 0.25; 0.15]) > 0.01


text:
using Documenter, Lifting

import Lifting:
    Lifting_Aux, Lifting_Exercise_Names, Lifting_Progressions, Lifting_Programmes

makedocs(;
    modules = [Lifting],
    format = Documenter.HTML(),
    pages = [
        "Home" => "index.md",
        "IO" => "io.md",
        "Types" => "types.md",
        "Functions & Methods" => "functions.md",
        "Defaults" => "dictionaries.md",
    ],
    repo = "https://github.com/dcelisgarza/Lifting.jl/blob/{commit}{path}#L{line}",
    sitename = "Lifting.jl",
    authors = "Daniel Celis Garza",
    # assets=String[],
)

deploydocs(; repo = "github.com/dcelisgarza/Lifting.jl.git")

max_stars_repo_name: dcelisgarza/Lifting


text:
@testset "GS" begin

m = 20
n = 10
A = randn(m,n);

Aold  = copy(A)
Qc,Rc = cgs(A)
@test all(Aold.==A) # make sure input is unchanged
Qt,Rt = cgs!(A)
@test all(Qt.==A) # check that input IS overwritten
@test maximum(abs.(Qt-Qc)) < 1e-15
@test norm(Aold-Qc*Rc,Inf)/norm(A,Inf) < 1e-14
@test norm(Qc'*Qc-I,Inf) < 1e-14

Aold  = copy(A)
Qm,Rm = mgs(A)
@test all(Aold.==A) # make sure input is unchanged
Qt,Rt = mgs!(A)
@test all(Qt.==A) # check that input IS overwritten
@test maximum(abs.(Qt-Qm)) < 1e-15
@test norm(Aold-Qm*Rm,Inf)/norm(A) < 1e-14
@test norm(Qm'*Qm-I,Inf) < 1e-14

end


text:
export LogHomotopy

"""
    LogHomotopy(H::AbstractHomotopy)

Reparametrizes the homotopy ``H(x,t)`` as ``Ĥ(x,s) := H(x,e^{-s})``.
"""
struct LogHomotopy{H<:AbstractHomotopy} <: AbstractHomotopy
    homotopy::H
end

Base.size(H::LogHomotopy) = size(H.homotopy)
basehomotopy(H::LogHomotopy) = basehomotopy(H.homotopy)

struct LogHomotopyCache{HC} <: AbstractHomotopyCache
    cache::HC
end

function cache(H::LogHomotopy, x, s)
    t = exp(-s)
    c = cache(H.homotopy, x, t)
    LogHomotopyCache(c)
end

@propagate_inbounds function evaluate!(u, H::LogHomotopy, x, s, c::LogHomotopyCache)
    evaluate!(u, H.homotopy, x, exp(-s), c.cache)
end
@propagate_inbounds function jacobian!(U, H::LogHomotopy, x, s, c::LogHomotopyCache)
    jacobian!(U, H.homotopy, x, exp(-s), c.cache)
end
@propagate_inbounds function dt!(u, H::LogHomotopy, x, s, c::LogHomotopyCache)
    t = exp(-s)
    dt!(u, H.homotopy, x, t, c.cache)
    LinearAlgebra.rmul!(u, -t)
    u
end

@propagate_inbounds function evaluate_and_jacobian!(
    u,
    U,
    H::LogHomotopy,
    x,
    s,
    c::LogHomotopyCache,
)
    evaluate_and_jacobian!(u, U, H.homotopy, x, exp(-s), c.cache)
end

@propagate_inbounds function jacobian_and_dt!(
    U,
    u,
    H::LogHomotopy,
    x,
    s,
    c::LogHomotopyCache,
)
    t = exp(-s)
    jacobian_and_dt!(U, u, H.homotopy, x, t, c.cache)
    LinearAlgebra.rmul!(u, -t)

    nothing
end

function evaluate(H::LogHomotopy, x, s, c::LogHomotopyCache)
    evaluate(H.homotopy, x, exp(-s), c.cache)
end
function jacobian(H::LogHomotopy, x, s, c::LogHomotopyCache)
    jacobian(H.homotopy, x, exp(-s), c.cache)
end
function dt(H::LogHomotopy, x, s, c::LogHomotopyCache)
    t = exp(-s)
    u = dt(H.homotopy, x, t, c.cache)
    LinearAlgebra.rmul!(u, -t)
    u
end


text:
using DynamicGrids, Test, Colors, ColorSchemes, FieldDefaults
using DynamicGrids: normaliseframe, frametoimage, @Image, @Graphic, @Output
using ColorSchemes: leonardo

init = [8.0 10.0;
        0.0  5.0]

# Define a simple image output
@Image @Graphic @Output struct TestImageOutput{} <: ImageOutput{T} end
processor = Nothing
output = TestImageOutput(init, minval=0.0, maxval=10.0)

ruleset = Ruleset(Life())

# Test level normalisation
normed = normaliseframe(output, output[1])
@test normed == [0.8 1.0
                 0.0 0.5]

# Test greyscale Image conversion
@test frametoimage(ColorProcessor(zerocolor=RGB24(1.0,0.0,0.0)), output, ruleset, init, 1) ==
    [RGB24(0.8, 0.8, 0.8) RGB24(1.0, 1.0, 1.0)
     RGB24(1.0, 0.0, 0.0) RGB24(0.5, 0.5, 0.5)]

l0 = RGB24(get(leonardo, 0))
l05 = RGB24(get(leonardo, 0.5))
l08 = RGB24(get(leonardo, 0.8))
l1 = RGB24(get(leonardo, 1))
@test frametoimage(ColorProcessor(;scheme=leonardo), output, ruleset, init, 1) == [l08 l1
                                                                                l0 l05]
z0 = RGB24(1, 0, 0)
@test frametoimage(ColorProcessor(;scheme=leonardo, zerocolor=z0), output, ruleset, init, 1) == [l08 l1
                                                                                z0 l05]


text:
function MCPP_solve(data, graph)
    """Given data, solve a MRCPP problem with integer programming."""

    ### Constants linked to the problem ###

    Er = graph["req_edges"]
    Enr = graph["noreq_edges"]
    Ar = graph["req_arcs"]
    Anr = graph["noreq_arcs"]

    depot = data[1, :depots][1]

    # Return acs for edges
    Er_swap = copy(Er)
    Er_swap[:start_node], Er_swap[:end_node] = Er_swap[:end_node], Er_swap[:start_node]
    Enr_swap = copy(Enr)
    Enr_swap[:start_node], Enr_swap[:end_node] = Enr_swap[:end_node], Enr_swap[:start_node]

    # Edges as and return arcs
    E = vcat(Er, Enr, Er_swap, Enr_swap)

    # All arcs
    A = vcat(E, Ar, Anr)

    # All requested arcs
    R = vcat(Er, Er_swap, Ar)

    # Number of edges/arcs by category
    card_Er = nrow(Er)
    card_Enr = nrow(Enr)
    card_Er_swap = nrow(Er_swap)
    card_Enr_swap = nrow(Enr_swap)
    card_E = nrow(E)
    card_Ar = nrow(Ar)
    card_Anr = nrow(Anr)
    card_A = nrow(A)
    card_R = nrow(R)

    # Traversing cost
    trav_cost = A[:, :trav_cost]

    # Service cost
    serv_cost = R[:, :serv_cost]

    # Demand
    demand = R[:, :demand]

    ### Model definition ###
    model = Model(with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=0))

    ### Variables ###
    X_Enr = [0 for arc in 1:card_Enr]
    X_Enr_swap = [0 for arc in 1:card_Enr]
    X_Anr = [0 for arc in 1:card_Anr]
    X_Ar = [1 for arc in 1:card_Ar]
    @variable(model, X_Er[1:card_Er], Bin)
    @variable(model, X_Er_swap[1:card_Er_swap], Bin)
    X_R = vcat(X_Er, X_Er_swap,  X_Ar)
    X_tot = vcat(X_Er, X_Enr, X_Er_swap, X_Enr_swap, X_Ar, X_Anr)
    @variable(model, Y_Ar[1:card_Ar], Int)
    @variable(model, Y_Anr[1:card_Anr], Int)
    @variable(model, Y_Er[1:card_Er], Int)
    @variable(model, Y_Enr[1:card_Enr], Int)
    @variable(model, Y_Er_swap[1:card_Er_swap], Int)
    @variable(model, Y_Enr_swap[1:card_Enr_swap], Int)
    Y_tot = vcat(Y_Er, Y_Enr, Y_Er_swap, Y_Enr_swap, Y_Ar, Y_Anr)

    ### Objective ###
    @objective(model, Min, sum(X_R.*serv_cost)+sum(Y_tot.*trav_cost))

    ### Constraint ###
    # Positivity constraint
    @constraint(model, positivity, Y_tot .>= 0)
    # Service constraint on edges
    @constraint(model, [arc in 1:card_Er], X_R[arc]+X_R[arc+card_Er] == 1)

    # Connectivity constraint
    @constraint(model, connectivity[i=1:data[1,:nodes]],
        sum([X_tot[arc]+Y_tot[arc] for arc in 1:card_A if A[arc,:start_node]==i]) ==
        sum([X_tot[arc]+Y_tot[arc] for arc in 1:card_A if A[arc,:end_node]==i])
        )

    # Depot contraint
    @constraint(model, depot_constraint, sum([X_tot[arc]+Y_tot[arc] for arc in 1:card_A if starts_in(arc, depot, A)])>=1)

    ### Solve problem ###

    start = now()
    optimize!(model)
    stop = now()

    solving_time = stop-start

    ### Format output ###

    tours = []

    cost = sum(value.(X_R).*serv_cost) + sum(value.(Y_tot).*trav_cost)

    X = Dict(
        "req_arcs" => [X_Ar[arc] for arc in 1:nrow(Ar)],
        "req_edges" => [value.(X_Er[arc]) for arc in 1:nrow(Er)],
        "req_edges_swap" => [value.(X_Er_swap[arc]) for arc in 1:nrow(Er_swap)]
        )

    Y = Dict(
        "req_arcs" => [value.(Y_Ar[arc]) for arc in 1:nrow(Ar)],
        "noreq_arcs" => [value.(Y_Anr[arc]) for arc in 1:nrow(Anr)],
        "req_edges" => [value.(Y_Er[arc]) for arc in 1:nrow(Er)],
        "noreq_edges" => [value.(Y_Enr[arc]) for arc in 1:nrow(Enr)],
        "req_edges_swap" => [value.(Y_Er_swap[arc]) for arc in 1:nrow(Er)],
        "noreq_edges_swap" => [value.(Y_Enr_swap[arc]) for arc in 1:nrow(Enr)]
        )

    tour = Dict(
        "depot" => depot,
        "cost" => cost,
        "capacity" => sum(value.(X_R).*demand),
        "X" => X,
        "Y" => Y
        )

    push!(tours, tour)

    return(tours, solving_time)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
module JungleHelperMossyWall

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/MossyWall" MossyWall(x::Integer, y::Integer, height::Integer=8, left::Bool=false, spriteDirectory::String="JungleHelper/Moss")

const placements = Ahorn.PlacementDict(
    "Mossy Wall (Right) (Jungle Helper)" => Ahorn.EntityPlacement(
        MossyWall,
        "rectangle",
        Dict{String, Any}(
            "left" => true
        )
    ),
    "Mossy Wall (Left) (Jungle Helper)" => Ahorn.EntityPlacement(
        MossyWall,
        "rectangle",
        Dict{String, Any}(
            "left" => false
        )
    )
)

mossColor = (51, 193, 17, 1) ./ (255, 255, 255, 1)

Ahorn.minimumSize(entity::MossyWall) = 0, 8
Ahorn.resizable(entity::MossyWall) = false, true

function Ahorn.selection(entity::MossyWall)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, 8, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::MossyWall, room::Maple.Room)
    left = get(entity.data, "left", false)

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    height = Int(get(entity.data, "height", 8))
    tileHeight = div(height, 8)

    spriteDirectory = get(entity.data, "spriteDirectory", "JungleHelper/Moss")

    if left
        for i in 2:tileHeight - 1
            Ahorn.drawImage(ctx, "$(spriteDirectory)/moss_mid1", 8, (i - 1) * 8, tint=mossColor)
        end

        Ahorn.drawImage(ctx, "$(spriteDirectory)/moss_top", 8, 0, tint=mossColor)
        Ahorn.drawImage(ctx, "$(spriteDirectory)/moss_bottom", 8, (tileHeight - 1) * 8, tint=mossColor)

    else
        Ahorn.Cairo.save(ctx)
        Ahorn.scale(ctx, -1, 1)

        for i in 2:tileHeight - 1
            Ahorn.drawImage(ctx, "$(spriteDirectory)/moss_mid1", 0, (i - 1) * 8, tint=mossColor)
        end

        Ahorn.drawImage(ctx, "$(spriteDirectory)/moss_top", 0, 0, tint=mossColor)
        Ahorn.drawImage(ctx, "$(spriteDirectory)/moss_bottom", 0, (tileHeight - 1) * 8, tint=mossColor)

        Ahorn.restore(ctx)
    end
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using PyCall

pyscf = pyimport("pyscf")
scipy = pyimport("scipy")
numpy = pyimport("numpy")
jax = pyimport("jax")
jnp = jax.numpy

jax.config.update("jax_enable_x64", true)

pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)
adscf = pyimport("adscf")

function cayley(mol)
    # overlap matrix
    S = mol.intor_symmetric("int1e_ovlp")
    # kinetic energy
    T = mol.intor_symmetric("int1e_kin")
    # potential energy
    V = mol.intor_symmetric("int1e_nuc")
    Hcore = T + V
    # two electron integrals
    JK = pyscf.ao2mo.restore(1, mol.intor("int2e"), mol.nao_nr())
    # the number of electrons
    Ne = mol.nelectron

    calcEnergy, gradEnergy = adscf.calcEnergy_create(S, Hcore, JK, Ne)

    # parameters
    tau = 1.0
    tau_m = 1e-10
    tau_M = 1e10
    rho = 1e-4
    delta = 0.1
    eta = 0.5
    epsilon = 1e-3
    max_iter = 5000

    # 1. initialize X0
    S64 = numpy.asarray(S, dtype=numpy.float64)
    X_np = scipy.linalg.inv(scipy.linalg.sqrtm(S64))
    X = jnp.asarray(X_np)

    # 2. set C=f(X0) and Q0=1
    C = calcEnergy(X)
    Q = 1.0

    # 4. calculate G0 and A0
    G = gradEnergy(X)
    A = jnp.matmul(jnp.matmul(G, X.T), S) - jnp.matmul(jnp.matmul(S, X), G.T)

    # main loop
    for it=1:40
        Y_tau = adscf.Y_tau_create(S)
        Y = Y_tau(tau, X, A)
        A_norm = jnp.linalg.norm(A, "fro")
        X_old, Q_old, G_old = X, Q, G

        # 5
        while numpy.greater(calcEnergy(Y), C - rho * tau * A_norm^2.0)
            tau *= delta    # 6
            Y = Y_tau(tau, X, A)
        end

        # 8
        X_new = Y
        Q_new = eta * Q + 1.0
        C = (eta * Q * C + calcEnergy(X_new)) / Q_new
        G_new = gradEnergy(X_new)
        Sk = X_new - X
        Yk = G_new - G
        # 9
        if it % 2 == 0
            tau_k = jnp.trace(jnp.matmul(Sk.T, Sk)) / abs(jnp.trace(jnp.matmul(Sk.T, Yk)))
        else
            tau_k = abs(jnp.trace(jnp.matmul(Sk.T, Yk))) / jnp.trace(jnp.matmul(Yk.T, Yk))
        end
        tau = max(min(tau_k, tau_M), tau_m)

        # Check loop condition
        X, Q, G = X_new, Q_new, G_new
        A = jnp.matmul(jnp.matmul(G, X.T), S) - jnp.matmul(jnp.matmul(S, X), G.T)
        c = jnp.linalg.norm(jnp.matmul(A, X))
        if numpy.less(c, epsilon)
            println(it)
            break
        end
    end
    e = calcEnergy(X)+mol.energy_nuc()
    return convert(Float64, e)
end

function main()
    mol = pyscf.gto.Mole()
    mol.build(atom = "H 0.0 0.0 0.0; H 0.0 0.0 1.4", basis="STO-3G", unit="Angstrom")
    e = @time cayley(mol)
    println("Energy: ", e)
end

main()


text:
my_f(x,y) = 2x+3y

# a closure is used -> see 18337 notes (18.337J Parallel Computing)
derivative_of_my_f(x,y) = ForwardDiff.derivative(x->my_f(x,y),x)


text:
#################### Empirical Move ####################

#################### Types and Constructors ####################

mutable struct EmpiricalTune <: SamplerTune
    samplerfun::Union{Function, Missing}
    width::Int64
    replacement::Bool

    EmpiricalTune() = new()

    EmpiricalTune(samplerfun::Union{Function, Missing}, width::Int64, replacement::Bool=false)=
        new(samplerfun, width, false)
end

EmpiricalTune(x::Vector, width::Int64, replacement::Bool=false)=
    EmpiricalTune(missing, width, replacement)

const EmpiricalVariate = SamplerVariate{EmpiricalTune}

#################### Sampler Constructor ####################
"""

function Empirical(params::ElementOrVector{Symbol}, width::Int64; args...)

This Sampler generates Samples from a distribution by selecting a random value
from this distribution. The move will always be accepted.
"""
function Empirical(params::ElementOrVector{Symbol},
                width::Int64;
                args...)
    samplerfx = function(model::Model, block::Integer)
        block = SamplingBlock(model, block)
        v = SamplerVariate(block, width, args...)
        sample!(v, x -> rand!(block, x))
        relist(block, v)
    end
    Sampler(params, samplerfx, EmpiricalTune())
end

sample!(v::EmpiricalVariate) = sample(v, v.tune.samplerfun)
function sample!(v::EmpiricalVariate, sampler_func::Function)
    width = v.tune.width
    sam = sampler_func(width)
    if width > 1 && v.tune.replacement
        while length(Set(sam)) != width
            sam = sampler_func(width)
        end
    end
    v[:] = sam
    v
end


text:
using Comrade, ComradeAdaptMCMC, Distributions
using Test

load_ehtim()
include(joinpath(@__DIR__, "../../../test/test_util.jl"))

@testset "ComradeAdaptMCMC.jl" begin

    m, vis, amp, lcamp, cphase = load_data()
    prior = test_prior()
    lklhd = RadioLikelihood(lcamp, cphase)
    post = Posterior(lklhd, prior, test_model)
    a1 = AdaptMCMC(ntemp=5)
    a2 = AdaptMCMC(ntemp=5, all_levels=true)

    x0 = prior_sample(post)

    chain, stats = sample(post, a1, 100_000; thin=5)
    chain, stats = sample(post, a2, 100_000; thin=5)
    chain, stats = sample(post, a1, 100_000; thin=5, init_params=x0)

    cpost = ascube(post)
    l0 = logdensityof(cpost, Comrade.HypercubeTransform.inverse(cpost, x0))
    @test l0 < mean(stats.logl[1])
end


text:
using Base.Test
# weighting.jl
import WaveletScattering: EqualWeighting, LoudnessWeighting, weight_frequencies

# EqualWeighting
weighting = EqualWeighting()
@test weight_frequencies(weighting, [0.4, 0.2, 0.1]) ≈ [1.0, 1.0, 1.0]

# LoudnessWeighting
weighting = LoudnessWeighting(44100)
@test weight_frequencies(weighting, 0.0) ≈ 0.0
@test weight_frequencies(LoudnessWeighting(12000.0), [1e48]) ≈ [0.0]


text:
module fug

using Plots
using JuMP
using Ipopt

#geometric mean of relative volitilies, useful for shortcut column
"""
geoMean(kLKd::Float64, kLKb::Float64, kHKd::Float64, kHKb::Float64)

returns the geometric mean of k values

```math
\\sqrt{\\frac{LK_d \\cdot LK_b}{HK_d \\cdot HK_b}}
```

# Arguments

*'kLKd'*: k value of light key in distilate

*'kLKb'*: k value of light key in bottoms

*'kHKd'*: k value of heavy key in distilate

*'kHKb'*: k value of heavy key in bottoms
"""
function geoMean(kLKd, kLKb, kHKd, kHKb)
    return(sqrt((kLKd*kLKb)/(kHKd*kHKb)))
end

#Fenske equation used in step 3 of column sizing to get the
"""
fenske(xLKd, kLKd, xHKd, kHKd, xLKb, kLKb, xHKb, kHKb)

returns the theoretical minimun number of trays for the distillation column

```math
n = \\frac{log_{10}(\\frac{xLKd \\cdot xHKb}{xLKb \\cdot xHKd})}{log_{10}(\\alpha_{avg})}
```

# Arguments
*'xLKd'*: x value of light key in distilate

*'kLKd'*: k value of light key in distilate

*'xHKd'*: x value of heavy key in distilate

*'kHKd'*: k value of heavy key in distilate

*'xLKb'*: x value of light key in bottoms

*'kLKb'*: k value of light key in bottoms

*'xHKb'*: x value of heavy key in bottoms

*'kHKb'*: k value of heavy key in bottoms
"""
function fenske(xLKd, kLKd, xHKd, kHKd, xLKb, kLKb, xHKb, kHKb)
    α_avg = geoMean(kLKd, kLKb, kHKd, kHKb)
    n = log10(xLKd*xHKb/(xLKb*xHKd)) / log10(α_avg)
    return(n) #return number of trays
end

#nonKeyComp finds the compistions of the non-key components in both the distillate and bottoms
"""
nonKeyComp(alpha, nmin, F, D, B, xHKd, xHKb, xfi)

returns the molar fraction of a non-key component in both the bottoms and the distillate

(xBi, xDi)

by solving the 2 equations at the same time:
```math
\\frac{x_{di}}{x_{bi}} = \\alpha^{n_{min}} \\cdot \\frac{x_{HK,d}}{x_{HK,b}}
```

```math
F_i = D_i + B_i
```

done by using the following optimization:
```math
Min[(\\frac{\\alpha^{N_{min}} \\cdot x_{HK,d}}{x_{HK,b}} - \\frac{x_{i,d} \\cdot B}{F \\cdot x_{i,f} - D \\cdot x_{i,d}})^2]
```

# Arguments
*'alpha'*: relative volitility of component with the heavy key

*'nmin'*: the theoretical minimum number of stages for the column, get from Fenske equation

*'F'*: feed molar flow rate

*'D'*: distillate molar flow rate

*'B'*: bottoms molar flow rate

*'xHKd*': x value of the heavy key in the distillate

*'xHKb'*: x value of the heavy key in the bottoms

*'xfi'*: the x value of the non-key component in the feed
"""
function nonKeyComp(α, nmin, F, D, B, xHKd, xHKb, xfi)
    model = Model(with_optimizer(Ipopt.Optimizer))
    set_silent(model)
    @variable(model, xdi)
    @NLconstraint(model, xdi >= 1e-9)
    @NLconstraint(model, xdi <= xfi*F/D)
    @NLconstraint(model, (xfi*F-D*xdi)/B <= 1)
    @NLobjective(model, Min, (α^nmin*xHKd/xHKb - xdi*B/(F*xfi-D*xdi))^2)
    optimize!(model)
    xDi = value(xdi)
    xBi = (F*xfi-D*value(xdi))/B
    empty!(model)
    return(xBi, xDi)
end

#Underwood equation for finding theta
"""
underwood(kDist, kBot, kHKd, kHKb, xFeed, qual, xDist)

The Underwood equation for finding the minimum reflux ratio by solving the first equation subject to the second condition:

```math
\\sum_{i=1}^n \\frac{\\alpha_i \\cdot x_{i,f}}{\\alpha_i - \\Theta} = 1-q
```
```math
R_{min} = \\sum_{i=1}^n \\frac{\\alpha_i \\cdot x_{i,d}}{\\alpha_i - \\Theta} - 1
```

# Arguments
*'kDist'*: vector of k values for all components in the distillate

*'kBot'*: vector of k values for all components in the bottoms

*'kHKd'*: k value of heavy key in distillate

*'kHKb'*: k value of heavy key in bottoms

*'xFeed'*: vector of x values for the feed

*'qual'*: quality of the feed stream (percentage of vapor if in VLE)

*'xDist'*: vector of x values for the distillate
"""
function underwood(kDist, kBot, kHKd, kHKb, xFeed, qual, xDist)
    l = length(kDist)
    α = zeros(l)
    for i = 1:l
       α[i] = geoMean(kDist[i], kBot[i], kHKd, kHKb)
    end

    model = Model(with_optimizer(Ipopt.Optimizer))
    set_silent(model)
    @variable(model, Θ)
    @NLconstraint(model, 1.01 <= Θ <= maximum(α))
    @NLobjective(model, Min, (sum((α[i]*xFeed[i])/(α[i]-Θ) for i in 1:l) + qual -1)^2)
    optimize!(model)

    rmin = sum(((α .* xDist) ./ (α .- value(Θ)))) - 1
    empty!(model)
    return(rmin)
end

#Gilliand correlation for FUG method
"""
gilliand(ract, rmin, nmin)

The Gilliand correlation for the actual number of stages. Solve the following equation for N

```math
\\frac{N-N_{min}}{N+1} = 0.75 \\cdot [1-(\\frac{R-R_{min}}{R+1})^{0.566}]
```

# Arguments
*'ract'*: The actual reflux ratio, typically 1.1-1.4 times the minimum reflux ratio

*'rmin'*: Minimum reflux ratio (Underwood equation)

*'nmin'*: Theoretical minimum number of stages (Fenske equation)
"""
function gilliand(ract, rmin, nmin)
    model = Model(with_optimizer(Ipopt.Optimizer))
    set_silent(model)
    @variable(model, Nact)
    @NLconstraint(model, nmin <= Nact)
    @NLobjective(model, Min, (0.75*(1-((ract-rmin)/(ract+1))^0.566) - (Nact-nmin)/(Nact+1))^2)
    optimize!(model)
    nact = value(Nact)
    empty!(model)
    return(nact)
end

#Kirkbride function for finding optimal feed stage, uses ipopt solver
"""
kirkbride(B,D,xHKf,xLKf,xHKd,xLKb,nTot)

The Kirkbride function for finding the optimial feed stage. Solves the first equation subject to the constraint:

```math
ln(\\frac{N_d}{N_b}) = 0.206 \\cdot ln(\\frac{B \\cdot x_{HK,f} \\cdot x_{LK,b}^2}{D \\cdot x_{LK,f} \\cdot x_{HK,d}^2})
```
```math
N_{tot} = N_d + N_b
```

# Arguments
*'B'*: bottoms molar flow rate

*'D'*: distillate molar flow rate

*'xHKf'*: x value of heavy key in feed

*'xLKf'*: x value of light key in feed

*'xHKd'*: x value of heavy key in distillate

*'xLKb'*: x value of light key in bottoms

*'nTot'*: total number of stages in column
"""
function kirkbride(B,D,xHKf,xLKf,xHKd,xLKb,nTot)
    t = ((B/D) * (xHKf/xLKf) * (xLKb/xHKd)^2)^0.206
    model = Model(with_optimizer(Ipopt.Optimizer))
    set_silent(model)
    @variable(model, m) # where m is number of stages above feed
    @NLconstraint(model, 0.2 <= m <= nTot-0.2)
    @NLobjective(model, Min, (m/(nTot-m) - t)^2)
    optimize!(model)
    M = value(m)
    empty!(model)
    return(M, nTot-M) #return number stages above, number of stages below
end

"""
diameter(F_ha, F_f, F_st, C_st, den_l, den_g, f, A_d, A_t, G)

function for finding the diameter of the column based on fluid properities. Solves equation for:

```math
D_T = \\sqrt{\\frac{4G}{fU_f \\pi (1-(A_d/A_T)) \\rho_g}}
```
```math
U_f = C \\cdot \\sqrt{\\frac{\\rho_l-\\rho_g}{\\rho_g}}
```
```math
C = C_{st} \\cdot F_{st} \\cdot F_f \\cdot F_{ha}
```

# Arguments
*'F_ha'*: Hole factor, 1 for valve and bubble cap trays;
"""
function diameter(F_ha, F_f, F_st, C_st, ρ_l, ρ_g, f, A_d, A_t, G)
    C = F_ha * F_f * F_st * C_st
    Uf = C*sqrt((ρ_l-ρ_g)/ρ_g)
    Dt = sqrt(4*G/ (f*Uf*pi*(1-A_d/A_t)*ρ_g))
    return(Dt)
end

### This is the master function for the fug method. The inputs will be:
### comp: List of component names
### xFeed: the component fractions in order from the least to the most volitile
### feedBasis: The molar flow rate of the feed
### kDist: the K values for the distillate
### kBot: the K values for the bottoms
### LKi: the index of the light key
### HKi: in the index of the heavy key
### splitFracLK: The percentage of the light key in the distillate for a 90/10 split this would be 0.9
### splitFracHK: the percentage of the heavy key in the bottoms
### F is molar flow rate of feed
### B is molar flow rate of bottoms
### D is molar flow rate of distillate
### quality is vapor fraction of the feed
### RRHeuristic is what to multiply the minimum reflux ratio by to get the actual reflux ratio

"""
fugk(comp, xFeed, kBot, kDist, HKi, LKi, splitFracHK, splitFracLK, F, B, D, quality, RRHeuristic)

main function for the Fenske-Underwood-Gilliand-Kirkbride function. Please make sure to put your components in the correct order.

# Arguments

*'comp'*: vector of component names in order from lowest to highest volitility

*'xFeed'*: vector of fractions for each component in comp

*'kBot'*: vector of k values for the bottoms of the column from lowest to highest

*'kDist'*: vector of k values for the distillate of the column from lowest to highest

*'HKi'*: Index of the heavy key (ie index in the comp list)

*'LKi'*: Index of the light key (ie index in the comp list)

*'splitFracHK'*: The fraction of the heavy key that will go to the bottoms. If 90\\% of the heavy key goes to the bottoms then 0.9

*'splitFracLK'*: The fraction of the light key that will go to the distilate

*'F'*: Feed molar feed rate

*'B'*: Bottoms molar flow rate

*'D'*: Distilate molar flow rate

*'quality'*: vapor fraction of the feed

*'RRHeuristic'*: reflux ratio heuristic value, typically in range [1.2,1.4]
"""
function fugk(comp, xFeed, kBot, kDist, HKi, LKi, splitFracHK, splitFracLK, F, B, D, quality, RRHeuristic)
    len = length(xFeed)
    feedMolFlow = xFeed * F
    botMolFlow=zeros(len); distMolFlow=zeros(len); xBot=zeros(len); xDist=zeros(len)

    botMolFlow[1:HKi-1]=feedMolFlow[1:HKi-1]
    botMolFlow[HKi]=feedMolFlow[HKi]*splitFracHK
    botMolFlow[LKi]=feedMolFlow[LKi]*(1-splitFracLK)

    distMolFlow[HKi]=feedMolFlow[HKi]*(1-splitFracHK)
    distMolFlow[LKi]=feedMolFlow[LKi]*splitFracLK
    distMolFlow[LKi+1:end]=feedMolFlow[LKi+1:end]

    for i ∈ 1:len
        xDist[i] = distMolFlow[i]/sum(distMolFlow)
        xBot[i] = botMolFlow[i]/sum(botMolFlow)
    end
    print(xDist[HKi], " x val of HK in distilate " , xBot[LKi], " x val of LK in bottoms \n")

    xLKd = xDist[LKi]; kLKd = kDist[LKi]; xLKb = xBot[LKi]; kLKb = kBot[LKi];
    xHKd = xDist[HKi]; kHKd = kDist[HKi]; xHKb = xBot[HKi]; kHKb = kBot[HKi];

    nmin = fenske(xLKd, kLKd, xHKd, kHKd, xLKb, kLKb, xHKb, kHKb)
    print("Theoretical min number of stages is ", nmin, '\n')

    names=Array{Any}(undef,len-2); kNKd=zeros(len-2); kNKb=zeros(len-2); xfi=zeros(len-2);
    names[1:HKi-1] = comp[1:HKi-1]; names[HKi:end] = comp[LKi+1:end];
    kNKd[1:HKi-1] = kDist[1:HKi-1]; kNKd[HKi:end] = kDist[LKi+1:end];
    kNKb[1:HKi-1] = kBot[1:HKi-1]; kNKb[HKi:end] = kBot[LKi+1:end];
    xfi[1:HKi-1] = xFeed[1:HKi-1]; xfi[HKi:end] = xFeed[LKi+1:end];
    xHKd = xDist[HKi]
    xHKb = xBot[HKi]

    alp = zeros(len-2)
    for i = 1:len-2
         alp[i] = geoMean(kNKd[i], kNKb[i], kHKd, kHKb)
    end

    xDist=zeros(len-2); xBot=zeros(len-2);
    for i = 1:len-2
        xBot[i], xDist[i] = nonKeyComp(alp[i], nmin, F, D, B, xHKd, xHKb, xfi[i])
        print(names[i], " bottoms frac ", xBot[i], " distilate frac ", xDist[i], '\n')
    end

    insert!(xBot,HKi,xHKb); insert!(xBot,LKi,xLKb);
    insert!(xDist,HKi,xHKd); insert!(xDist,LKi,xLKd);

    rmin = underwood(kDist, kBot, kHKd, kHKb, xFeed, quality, xDist)
    print("The min reflux ratio is ", rmin, '\n')
    ract = rmin * RRHeuristic

    nact = gilliand(ract, rmin, nmin)
    print("Our actual number of stages will be ", nact, '\n')

    above, below = kirkbride(B,D,xFeed[HKi],xFeed[LKi],xHKd,xLKb,nact)
    print("our optimal feed stage is ", above, " from the top and ", below, " from the bottom \n")

end

end # module


text:
# This file is a part of SimilaritySearch.jl

"""
struct SearchGraphCallbacks
    hints::Union{Nothing,Callback} = DisjointHints()
    hyperparameters::Union{Nothing,Callback} = OptimizeParameters(kind=ParetoRecall())
    logbase::Float32 = 1.5
    starting::Int32 = 8
end

Call insertions and indexing methods with `SearchGraphCallbacks` objects to control how the index structure is adjusted (callbacks are called when ``n > starting`` and ``\\lceil(\\log(logbase, n)\\rceil != \\lceil\\log(logbase, n+1)\\rceil``)
"""
@with_kw struct SearchGraphCallbacks
    hints::Union{Nothing,Callback} = DisjointHints()
    hyperparameters::Union{Nothing,Callback} = OptimizeParameters(kind=ParetoRecall())
    logbase::Float32 = 1.5
    starting::Int32 = 8
end

Base.copy(g::SearchGraph;
    dist=g.dist,
    db=g.db,
    links=g.links,
    locks=g.locks,
    hints=g.hints,
    search_algo=copy(g.search_algo),
    verbose=true
) = SearchGraph(; dist, db, links, locks, hints, search_algo, verbose)

"""
    SearchGraphCallbacks(kind::ErrorFunction;
        hints=DisjointHints(),
        logbase=1.5,
        starting=8,
        initialpopulation=16,
        maxpopulation=initialpopulation,
        maxiters=12,
        bsize=4,
        mutbsize=4bsize,
        crossbsize=2bsize,
        tol=-1.0,
        params=SearchParams(; maxpopulation, bsize, mutbsize, crossbsize, tol, maxiters),
        ksearch=10,
        numqueries=32,
        space::BeamSearchSpace=BeamSearchSpace()
    )

Convenient constructor function to create `SearchGraphCallbacks` structs. See [`SearchGraphCallbacks`](@ref), [`SearchParams`](@ref), and [`BeamSearchSpace`](@ref) for more details.

# Arguments

- `kind`: The kind of error function, e.g. `MinRecall(0.9)`.
- `hints`: How search hints should be computed.
- `logbase`: Controls the periodicity of callback executions.
- `starting`: Controls the minimum size of the index before execute callbacks.
- `initialpopulation`: Optimization argument that determines the initial number of configurations.
- `maxiters`: Optimization argument that determines the number of iterations.
- `bsize`: Optimization argument that determines how many top configurations are allowed to mutate and cross.
- `tol`: Optimization argument that determines the minimal tolerance improvement to stop the optimization.
- `params`: The `SearchParams` arguments (if separated optimization arguments are not enough)
- `ksearch`: The number of neighbors to be retrived by the optimization process.
- `numqueries`: The number of queries to be performed during the optimization process.
- `space`: The cofiguration search space
"""
function SearchGraphCallbacks(kind::ErrorFunction;
    hints=DisjointHints(),
    logbase=1.5,
    starting=8,
    initialpopulation=16,
    maxiters=12,
    bsize=4,
    tol=-1.0,
    params=SearchParams(; maxpopulation=initialpopulation, bsize, mutbsize=4bsize, crossbsize=2bsize, tol, maxiters),
    ksearch=10,
    numqueries=32,
    space::BeamSearchSpace=BeamSearchSpace()
)
hyperparameters = OptimizeParameters(; kind, initialpopulation, params, ksearch, numqueries, space)
SearchGraphCallbacks(; hints, hyperparameters, logbase, starting)
end

"""
    execute_callbacks(index::SearchGraph, n=length(index), m=n+1)

Process all registered callbacks
"""
function execute_callbacks(callbacks::SearchGraphCallbacks, index::SearchGraph, n=length(index), m=n+1; force=false)
    if force || (n >= callbacks.starting && ceil(Int, log(callbacks.logbase, n)) != ceil(Int, log(callbacks.logbase, m)))
        callbacks.hints !== nothing && execute_callback(callbacks.hints, index)
        callbacks.hyperparameters !== nothing && execute_callback(callbacks.hyperparameters, index)
    end
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
using SFML

# This example downloads the html for the webpage http://www.sfml-dev.org/learn.php
http = Http()

set_host(http, "http://www.sfml-dev.org")

request = HttpRequest()
set_uri(request, "learn.php")
set_method(request, SFML.HttpGet)

response = send_request(http, request)

status = get_status(response)
if status == SFML.HttpOk
    println("Success!")
    println("Got: $(get_body(response))")
else
    println("Error")
    println(status)
end


text:
using Gadfly
set_default_plot_size(21cm, 8cm)

p1 = plot(xmin=[1.0, 5.0, 7.0], xmax=[2.0, 6.5, 8.0] , Geom.vband, Theme(default_color="green"));

p2 = plot(ymin=[2.5], ymax=[7.5], Geom.hband, Theme(default_color="red"));

hstack(p1, p2)


text:
Base.one{N,ES}(T::Type{Valid{N,ES}}) = Valid(one(Vnum{N,ES}), one(Vnum{N,ES}))
Base.zero{N,ES}(T::Type{Valid{N,ES}}) = Valid(zero(Vnum{N,ES}), zero(Vnum{N,ES}))

type ∅; end
type ℝ; end
type ℝp; end

Base.convert{N, ES}(T::Type{Valid{N, ES}}, ::Type{∅})  = Valid(zero(     Vnum{N,ES}), maxneg( Vnum{N,ES}))
Base.convert{N, ES}(T::Type{Valid{N, ES}}, ::Type{ℝ})  = Valid(-maxpos( Vnum{N,ES}), maxpos(      Vnum{N,ES}))
Base.convert{N, ES}(T::Type{Valid{N, ES}}, ::Type{ℝp}) = Valid(Vnum{N,ES}(Inf),       maxpos(      Vnum{N,ES}))

type __plusstar; end
Base.:+(::typeof(Base.:*)) = __plusstar
type __minusstar; end
Base.:-(::typeof(Base.:*)) = __minusstar
type __positives; end
type __negatives; end

(::Type{ℝ})(::Type{__plusstar}) =      __plusstar
(::Type{ℝ})(::Type{__minusstar}) =     __minusstar
(::Type{ℝ})(::typeof(+)) =             __positives
(::Type{ℝ})(::typeof(-)) =             __negatives

Base.convert{N,ES}(T::Type{Valid{N,ES}}, ::Type{__plusstar}) =  Valid(minpos(Vnum{N, ES}),  maxpos(Vnum{N, ES}))
Base.convert{N,ES}(T::Type{Valid{N,ES}}, ::Type{__minusstar}) = Valid(-maxpos(Vnum{N, ES}), maxneg(Vnum{N, ES}))
Base.convert{N,ES}(T::Type{Valid{N,ES}}, ::Type{__positives}) = Valid(zero(Vnum{N, ES}),          maxpos(Vnum{N, ES}))
Base.convert{N,ES}(T::Type{Valid{N,ES}}, ::Type{__negatives}) = Valid(-maxpos(Vnum{N, ES}),         zero(Vnum{N, ES}))

export ∅, ℝ, ℝp

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
using DynamicGrids, FreeTypeAbstraction, Test

@testset "Fonts" begin
	@test DynamicGrids.autofont() isa FreeTypeAbstraction.FTFont
	face = DynamicGrids.autofont()
	name = face.family_name
 	@testset "TextConfig accepts font as String" begin
		@test name isa String
		textconfig = TextConfig(; font=name)
		@test textconfig.face isa FreeTypeAbstraction.FTFont
	end
 	@testset "TextConfig accepts font as FTFont" begin
		@test face isa FreeTypeAbstraction.FTFont
		textconfig = TextConfig(; font=face)
		@test textconfig.face === face
	end
	@test_throws ArgumentError TextConfig(; font="not_a_font")
	@test_throws ArgumentError TextConfig(; font=:not_a_string)
end


text:
"""
    OverdeterminedTracker

An `OverdeterminedTracker` allows to solve overdetermined systems by
first squaring them up, then solving them with a `PolyhedralTracker` or `PathTracker`
and finally to find all actual solutions.
"""
struct OverdeterminedTracker{
    T<:AbstractPathTracker,
    AV<:AbstractVector,
    S<:HomotopyWithCache,
    NC<:NewtonCorrectorCache,
} <: AbstractPathTracker
    tracker::T
    # stuff to verify containment
    system::S
    y::AV
    jacobian::JacobianMonitor{Float64}
    newton::NC
end

function construct_tracker(prob::OverdeterminedProblem, start_solutions; kwargs...)
    tracker = construct_tracker(prob.problem, start_solutions)
    y = copy(path_tracker_state(tracker).solution)
    system = HomotopyWithCache(ConstantHomotopy(prob.target_system), y, 0.0)
    J = JacobianMonitor(jacobian(system, y, 0.0))
    newton_corrector = NewtonCorrectorCache(system, y, 0.0)

    OverdeterminedTracker(tracker, system, y, J, newton_corrector)
end

seed(OT::OverdeterminedTracker) = seed(OT.tracker)
function PathResult(OT::OverdeterminedTracker, x, path_number = nothing; kwargs...)
    PathResult(OT.tracker, x, path_number; kwargs...)
end

result_type(OT::OverdeterminedTracker) = result_type(OT.tracker)
prepare!(OT::OverdeterminedTracker, S) = prepare!(OT.tracker, S)
function track!(
    OT::OverdeterminedTracker,
    x;
    accuracy::Union{Nothing,Float64} = nothing,
    max_corrector_iters::Union{Nothing,Int} = nothing,
)
    retcode = track!(
        OT.tracker,
        x;
        accuracy = accuracy, max_corrector_iters = max_corrector_iters,
    )
    is_success(retcode) || return retcode

    # We take the solution and try an overdetermined Newton method to see whether this also
    # converges
    s = solution(OT.tracker)
    result = newton!(
        OT.y,
        OT.system,
        s,
        0.0,
        OT.jacobian,
        norm(OT.tracker),
        OT.newton;
        tol = min_accuracy(OT.tracker),
        # Make sure to evaluate the Jacobian only *once*. Otherwise it can happen at singuar
        # solutions that we bounce away from a good solution.
        full_steps = 1,
        max_iters = 2,
        double_64_evaluation = false,
    )

    # overwrite the solution information used for constructing a `PathResult`
    state = path_tracker_state(OT)
    state.solution_cond = cond!(OT.jacobian, InfNorm())
    original_residual = state.solution_residual
    state.solution_residual = LA.norm(OT.newton.r, InfNorm())
    state.solution_accuracy = result.accuracy

    if is_converged(result)
        state.solution .= OT.y
        state.status = PathTrackerStatus.success
    else
        if state.solution_cond > 1e8
            evaluate!(OT.newton.r, OT.system, s, 0.0)
            state.solution_residual = LA.norm(OT.newton.r, InfNorm())
            if state.solution_residual / original_residual < 1e5
                state.status = PathTrackerStatus.success
            else
                state.status = PathTrackerStatus.excess_solution
            end
        else
            state.status = PathTrackerStatus.excess_solution
        end
    end
end

function track(
    tracker::OverdeterminedTracker,
    x;
    path_number::Union{Int,Nothing} = nothing,
    accuracy::Union{Nothing,Float64} = nothing,
    max_corrector_iters::Union{Nothing,Int} = nothing,
    details::Symbol = :default,
)
    track!(tracker, x; accuracy = accuracy, max_corrector_iters = max_corrector_iters)
    PathResult(tracker, x, path_number; details = details)
end

path_tracker_state(OT::OverdeterminedTracker) = path_tracker_state(OT.tracker)
path_tracker_state(PT::PolyhedralTracker) = path_tracker_state(PT.generic_tracker)
path_tracker_state(T::PathTracker) = T.state


text:
using CxxWrap
using Base.Test

function func1(arr)
  @test arr[1] == 1.0
  @test arr[2] == 2.0
  @test arr[3] == 3.0
end

wrap_modules(CxxWrap._l_containers)
using Containers

#println(Base.uncompressed_ast(first(methods(test_tuple))))

@test test_tuple() == (1,2.0,3.0f0)

cptr = const_ptr()
@test isbits(typeof(cptr))
@test const_ptr_arg(cptr) == (1., 2., 3.)

cv = const_vector()
@test size(cv) == (3,)
@test cv == [1.,2.,3.]

cm = const_matrix()
@test size(cm) == (3,2)
@test cm == [[1.,2.,3.] [4.,5.,6.]]
println("Displaying const matrix")
display(cm)
println()

mm = Containers.mutable_array()
println("Displaying mutable matrix")
display(mm)
println()
mm[:,:] = 1.0
@test Containers.check_mutable_array(mm)
Containers.do_embedding_test()


text:
DBG = false

macro debug(msg)
    if DBG
        :(println($msg))
    end
end

function idem(needle::AbstractString, haystack::AbstractString, idx)
    for c1 in needle
        if done(haystack, idx)
            return false
        end
        c2, idx = next(haystack, idx)
        if c2 != c1
            return false
        end
    end
    return true
end

getchar(state::ParserState) = state.subject[state.index]

function nextchar!(state::ParserState)
    if done(state.subject, state.index)
        return :eof
    end
    (char, state.index) = next(state.subject, state.index)
    @debug ("nextchar: $char")
    char
end

function next_non_space!(state::ParserState)
    while true
        c = nextchar!(state)
        if c != ' ' && c != '\t'
            return c
        end
    end
end

function endlineP(c, state::ParserState) # Where 'P' really is a '?'.
    if c == '\r' && state.subject[state.index] == '\n'
        nextchar!(state)
    end
    if c == '\r' || c == '\n'
        state.line += 1
        return true
    else
        return false
    end
end

function endline!(state::ParserState)
    c = next_non_space!(state);
    if endlineP(c, state)
        return
    elseif c == :eof
        return
    elseif c == '#'
        while (c = nextchar!(state); !endlineP(c, state) && c != :eof)
        end
        return
    else
        _error("Illegal character \'$c\'", state)
    end
end

function next_non_comment!(state::ParserState)
    local in_comment = false
    while true
        c = nextchar!(state)
        if in_comment
            if c == :eof
                return c
            elseif endlineP(c, state)
                in_comment = false
            end
        else
            if c == '#'
                in_comment = true
            elseif !endlineP(c,state) && c != ' ' && c != '\t'
                return c
            end
        end
    end
end

immutable TOMLError
    msg
end

Base.show(io::IO, e::TOMLError) = print(io, repr(e))
Base.repr(e::TOMLError) = "TOMLError: $(e.msg)"

_error(msg, state) = throw(TOMLError("$msg on line $(state.line)."))


text:
using Test, FinancialMonteCarlo, DualNumbers;
@show "Black Scholes Model"
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

underlying_name = "ENJ"

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

display(Model)

port_ = EUData * 2.0 + FwdData - FwdData / 2.0 - FwdData;
port_ = port_ * 2.0;
port_ = 2.0 * port_;
display(port_)
print(port_)
display(FwdData)
display(-FwdData)
print(FwdData)
print(-FwdData)
display(Forward(dual(T, 1.0)))
display(-Forward(dual(T, 1.0)))
print(Forward(dual(T, 1.0)))
print(-Forward(dual(T, 1.0)))

display(Forward(dual(T, 1.0)) + EUData)
display(-Forward(dual(T, 1.0)) + EUData)
print(Forward(dual(T, 1.0)) + EUData)
print(-Forward(dual(T, 1.0)) + EUData)

display(Forward(dual(T, 1.0)) - EUData)
display(-Forward(dual(T, 1.0)) - EUData)
print(Forward(dual(T, 1.0)) - EUData)
print(-Forward(dual(T, 1.0)) - EUData)
fwd = Forward(dual(T, 1.0)) * 1.0
fwd = 1.0 * Forward(dual(T, 1.0))
ok = FwdData - fwd;
display(fwd)
print(fwd)
fwd = 1.0 * fwd
fwd = fwd * 1.0
fwd = fwd / 1.0
fwd = -fwd
fwd = Forward(dual(T, 1.0)) * dual(1.0, 1.0) + EUData * 2.0 + AMData
display(fwd)
print(fwd)
fwd = -Forward(dual(T, 1.0)) * dual(1.0, 1.0) + EUData * 2.0 + AMData
display(fwd)
print(fwd)

mktdataset = underlying_name → Model
portfolio_ = [FwdData; EUData; AMData; BarrierData; AsianFloatingStrikeData; AsianFixedStrikeData; Spot()];
portfolio = underlying_name → FwdData
portfolio += underlying_name → EUData
portfolio += underlying_name → AMData
portfolio += underlying_name → Spot()
portfolio += underlying_name → BarrierData
portfolio += underlying_name → AsianFloatingStrikeData
portfolio += underlying_name → AsianFixedStrikeData

display(portfolio)

price_mkt = pricer(mktdataset, rfCurve, mc, portfolio)
price_old = sum(pricer(Model, rfCurve, mc, portfolio_))

@test abs(price_mkt - price_old) < 1e-8


text:
# Utility functions for common algebraic operations.

"""
    rotateTensorToVectorZ(tensor, vector)

Rotates `tensor` with a rotation matrix that aligns z-axis with `vector`.
"""
function rotateTensorToVectorZ(tensor::AbstractMatrix{T}, v::AbstractVector{T}) where T
   unitz = SVector{3, T}(0.0, 0.0, 1.0)
   vz = v × unitz::SVector{3, T}
   vz ./= hypot(vz[1], vz[2], vz[3])
   angle = acos(v ⋅ unitz / hypot(v[1], v[2], v[3]))
   R = getRotationMatrix(vz, angle)
   # Rotate tensor
   R * tensor * R'
end

"""
    getRotationMatrix(vector, angle)

Creates a rotation matrix that rotates around a unit `vector` by an `angle` in radians.
References: https://en.wikipedia.org/wiki/Rodrigues'_rotation_formula
https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
"""
function getRotationMatrix(v::AbstractVector, θ)
   sinθ, cosθ = sincos(eltype(v)(θ))
   tmp = 1 - cosθ
   m =  @SMatrix [
        cosθ+v[1]^2*tmp         v[1]*v[2]*tmp-v[3]*sinθ v[1]*v[3]*tmp+v[2]*sinθ;
        v[1]*v[2]*tmp+v[3]*sinθ cosθ+v[2]^2*tmp         v[2]*v[3]*tmp-v[1]*sinθ;
        v[1]*v[3]*tmp-v[2]*sinθ v[3]*v[2]*tmp+v[1]*sinθ cosθ+v[3]^2*tmp]
end

"""
    getRotationB(B) -> SMatrix

Obtain a rotation matrix with each column being a unit vector which is parallel (`v3`) and
perpendicular (`v1,v2`) to the magnetic field `B`. The two perpendicular directions are
chosen based on the reference vector of z-axis in the Cartesian coordinates.
"""
function getRotationB(B::AbstractVector{T}) where T
   b = hypot(B[1], B[2], B[3])
   if B[3] / b ≈ -1.0 # B aligned with reference vector
      R = SMatrix{3,3,T}([0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 1.0])
   elseif B[3] / b ≈ 1.0 # B aligned with reference vector
      R = SMatrix{3,3,T}([0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 1.0])
   else
      v0 = SVector{3,T}(0.0, 0.0, 1.0) # reference vector
      v3 = SVector{3,T}(B[1]/b, B[2]/b, B[3]/b) # unit vector along B
      v1 = v0 × v3::SVector{3, T}
      v2 = v3 × v1::SVector{3, T}

      R = @SMatrix [v1[1] v2[1] v3[1]; v1[2] v2[2] v3[2]; v1[3] v2[3] v3[3]]
   end
   R
end

"""
    rotateWithB(T, B) -> Matrix

Rotate the tensor `T` with the 3rd direction aligned with `B`.
See also: [`rotateWithB!`](@ref)
"""
function rotateWithB(T::AbstractMatrix, B::AbstractVector)
   R = getRotationB(B)
   R * T * R'
end

"""
    rotateWithB!(T, B)

Rotate the tensor `T` with the 3rd direction aligned with `B`.
See also: [`rotateWithB`](@ref)
"""
function rotateWithB!(T::AbstractMatrix, B::AbstractVector)
   R = getRotationB(B)
   T[:] = R * T * R'
end


text:
@testset "LDA" begin

    doc1 = "a a a sample text text"
    doc2 = "another example example text text"

    crps = Corpus(Any[StringDocument(doc1), StringDocument(doc2)])

    update_lexicon!(crps)

    dtm = DocumentTermMatrix(crps)
    ϕ, θ = lda(dtm, 2, 25, 0.1, 0.1)
    @test ϕ isa SparseMatrixCSC
    @test θ isa Matrix{Float64}
    @test all(sum(θ, dims=1) .≈ 1)
end


text:
module heatbath_module
using LinearAlgebra

import ..AbstractGaugefields_module:normalize3!,normalizeN!,AbstractGaugefields,evaluate_gaugelinks_evenodd!, map_U!,map_U_sequential!
import Wilsonloop:loops_staple
import ..GaugeAction_module:GaugeAction,evaluate_staple_eachindex!
using InteractiveUtils

struct Heatbath{T}
    _tempotal_gauges::Vector{T}
    β::Float64
    ITERATION_MAX::Int64

    function Heatbath(U::Array{T,1},β;ITERATION_MAX=10^5)  where T <: AbstractGaugefields
        _tempotal_gauges = Array{T,1}(undef,3)
        for i=1:3
            _tempotal_gauges[i] = similar(U[1])
        end
        return new{T}(_tempotal_gauges,β,ITERATION_MAX)
    end

end

struct Heatbath_update{Dim,T}
    _temporary_gaugefields::Vector{T}
    gauge_action::GaugeAction{Dim,T}
    ITERATION_MAX::Int64

    function Heatbath_update(U::Array{T,1},gauge_action;ITERATION_MAX=10^5) where T <: AbstractGaugefields
        _temporary_gaugefields = Array{T,1}(undef,3)
        Dim = length(U)
        for i=1:3
            _temporary_gaugefields[i] = similar(U[1])
        end
        return new{Dim,T}(_temporary_gaugefields,gauge_action,ITERATION_MAX)
    end
end

const heatbath_factor = 2

function heatbath_update_eachsite_SU2!(A,μ,U::Vector{<: AbstractGaugefields{NC,Dim}},h::Heatbath_update{Dim,T},mat_temps,indices...) where {NC,Dim,T}
    @assert NC == 2
    V = mat_temps[1]
    evaluate_staple_eachindex!(V,μ,h.gauge_action,U,view(mat_temps,2:6),indices...)
    SU2update_KP!(A,V,heatbath_factor ,NC,view(mat_temps,7:8),h.ITERATION_MAX)
end

function heatbath!(U::Array{<: AbstractGaugefields{2,Dim},1},h::Heatbath_update{Dim,T}) where {Dim,T}
    NC = 2
    nt = 8
    temps = Vector{Matrix{ComplexF64}}(undef,nt)
    for i=1:nt
        temps[i] = zeros(ComplexF64,NC,NC)
    end

    for μ=1:Dim
        mapfunc!(A,U,indices...) = heatbath_update_eachsite_SU2!(A,μ,U,h,temps,indices...)
        map_U_sequential!(U[μ],mapfunc!,U)
    end

end

function heatbath_update_eachsite_SU3!(A,μ,U::Vector{<: AbstractGaugefields{NC,Dim}},h::Heatbath_update{Dim,T},mat_temps1,mat_temps2,indices...) where {NC,Dim,T}
    @assert NC == 3
    V = mat_temps1[1]
    evaluate_staple_eachindex!(V,μ,h.gauge_action,U,view(mat_temps1,2:6),indices...)
    SU3update_matrix!(A,V,heatbath_factor ,NC,view(mat_temps1,7:11),mat_temps2,h.ITERATION_MAX)
    #SU2update_KP!(A,V,1,NC,view(mat_temps,7:8),h.ITERATION_MAX)
end

function heatbath_update_eachsite_SUN!(A,μ,U::Vector{<: AbstractGaugefields{NC,Dim}},h::Heatbath_update{Dim,T},mat_temps1,mat_temps2,indices...) where {NC,Dim,T}
    V = mat_temps1[1]
    evaluate_staple_eachindex!(V,μ,h.gauge_action,U,view(mat_temps1,2:6),indices...)
    SUNupdate_matrix!(A,V,heatbath_factor ,NC,view(mat_temps1,7:11),mat_temps2,h.ITERATION_MAX)
    #SU2update_KP!(A,V,1,NC,view(mat_temps,7:8),h.ITERATION_MAX)
end

function overrelaxation_update_eachsite_SUN!(A,μ,U::Vector{<: AbstractGaugefields{NC,Dim}},h::Heatbath_update{Dim,T},mat_temps1,mat_temps2,indices...) where {NC,Dim,T}
    V = mat_temps1[1]
    evaluate_staple_eachindex!(V,μ,h.gauge_action,U,view(mat_temps1,2:6),indices...)
    SUN_overrelaxation!(A,V,heatbath_factor ,NC,view(mat_temps1,7:11),mat_temps2,h.ITERATION_MAX)
    #SU2update_KP!(A,V,1,NC,view(mat_temps,7:8),h.ITERATION_MAX)
end

function heatbath!(U::Array{<: AbstractGaugefields{NC,Dim},1},h::Heatbath_update{Dim,T}) where {Dim,T,NC}
    nt = 11
    temps = Vector{Matrix{ComplexF64}}(undef,nt)
    for i=1:nt
        temps[i] = zeros(ComplexF64,NC,NC)
    end

    if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    end

    for μ=1:Dim
        mapfunc!(A,U,indices...) = heatbath_update_eachsite_SUN!(A,μ,U,h,temps,temps3,indices...)
        map_U_sequential!(U[μ],mapfunc!,U)
    end

end

function overrelaxation!(U::Array{<: AbstractGaugefields{NC,Dim},1},h::Heatbath_update{Dim,T}) where {Dim,T,NC}
    nt = 11
    temps = Vector{Matrix{ComplexF64}}(undef,nt)
    for i=1:nt
        temps[i] = zeros(ComplexF64,NC,NC)
    end

    #if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    #end

    for μ=1:Dim
        mapfunc!(A,U,indices...) = heatbath_update_eachsite_SUN!(A,μ,U,h,temps,temps3,indices...)
        map_U_sequential!(U[μ],mapfunc!,U)
    end

end

function heatbath!(U::Array{<: AbstractGaugefields{3,Dim},1},h::Heatbath_update{Dim,T}) where {Dim,T}
    NC = 3
    nt = 11
    temps = Vector{Matrix{ComplexF64}}(undef,nt)
    for i=1:nt
        temps[i] = zeros(ComplexF64,NC,NC)
    end

    if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    end

    for μ=1:Dim
        mapfunc!(A,U,indices...) = heatbath_update_eachsite_SU3!(A,μ,U,h,temps,temps3,indices...)
        map_U_sequential!(U[μ],mapfunc!,U)
    end

end

function heatbath!(U::Array{<: AbstractGaugefields{NC,Dim},1},h::Heatbath) where {Dim,NC}
    heatbath!(U,h._tempotal_gauges,h.β;ITERATION_MAX=h.ITERATION_MAX)
end

function overrelaxation!(U::Array{<: AbstractGaugefields{NC,Dim},1},h::Heatbath) where {Dim,NC}
    overrelaxation!(U,h._tempotal_gauges,h.β;ITERATION_MAX=h.ITERATION_MAX)
end

function heatbath!(U::Array{<: AbstractGaugefields{2,Dim},1},temps,β;ITERATION_MAX=10^5) where {Dim}
    NC = 2
    temp1 = temps[1]
    temp2 = temps[2]
    V = temps[3]

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
    end

    if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    end

    mapfunc!(A,B) = SU2update_KP!(A,B,β,NC,temps2,ITERATION_MAX)

    for μ=1:Dim

        loops = loops_staple[(Dim,μ)]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function heatbath!(U::Array{<: AbstractGaugefields{3,Dim},1},temps,β;ITERATION_MAX=10^5) where {Dim}
    NC = 3
    temp1 = temps[1]
    temp2 = temps[2]
    V = temps[3]

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
    end

    if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    end

    mapfunc!(A,B) = SU3update_matrix!(A,B,β,NC,temps2,temps3,ITERATION_MAX)

    for μ=1:Dim

        loops = loops_staple[(Dim,μ)]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function heatbath!(U::Array{<: AbstractGaugefields{3,Dim},1},temps,β,gauge_action;ITERATION_MAX=10^5) where {Dim} #This function is for debugging
    NC = 3
    temp1 = temps[1]
    temp2 = temps[2]
    V = temps[3]

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
    end

    if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    end

    mapfunc!(A,B) = SU3update_matrix!(A,B,β,NC,temps2,temps3,ITERATION_MAX)

    for μ=1:Dim

        #loops = loops_staple[(Dim,μ)]
        loops = gauge_action.dataset[1].staples[μ]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function heatbath!(U::Array{<: AbstractGaugefields{NC,Dim},1},temps,β;ITERATION_MAX=10^5) where {Dim,NC}

    temp1 = temps[1]
    temp2 = temps[2]
    V = temps[3]

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
    end

    if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    end

    mapfunc!(A,B) = SUNupdate_matrix!(A,B,β,NC,temps2,temps3,ITERATION_MAX)

    for μ=1:Dim

        loops = loops_staple[(Dim,μ)]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function overrelaxation!(U::Array{<: AbstractGaugefields{NC,Dim},1},temps,β;ITERATION_MAX=10^5) where {Dim,NC}

    temp1 = temps[1]
    temp2 = temps[2]
    V = temps[3]

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
    end

    #if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    #end

    mapfunc!(A,B) = SUN_overrelaxation!(A,B,β,NC,temps2,temps3,ITERATION_MAX) #SUNupdate_matrix!(A,B,β,NC,temps2,temps3,ITERATION_MAX)

    for μ=1:Dim

        loops = loops_staple[(Dim,μ)]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function heatbath!(U::Array{<: AbstractGaugefields{NC,Dim},1},h::Heatbath_update) where {Dim,NC}
    heatbath!(U,h._tempotal_gauges,h.gauge_action;ITERATION_MAX=h.ITERATION_MAX)
end

function heatbath!(U::Array{<: AbstractGaugefields{2,Dim},1},temps,S::GaugeAction;ITERATION_MAX=10^5) where {Dim}
    NC = 2

    temps2 = Array{Matrix{ComplexF64},1}(undef,5)
    for i=1:5
        temps2[i] = zeros(ComplexF64,2,2)
    end

    if NC != 2
        temps3 = Array{Matrix{ComplexF64},1}(undef,5)
        for i=1:5
            temps3[i] = zeros(ComplexF64,NC,NC)
        end
    end

    mapfunc!(A,B) = SU2update_KP!(A,B,1,NC,temps2,ITERATION_MAX)

    numterm = length(S.dataset)
    temp1 = S._temp_U[1]
    temp2 = S._temp_U[2]
    temp3 = S._temp_U[3]

    error("error")

    for μ=1:Dim

        loops = loops_staple[(Dim,μ)]
        iseven = true

        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)

        iseven = false
        evaluate_gaugelinks_evenodd!(V,loops,U,[temp1,temp2],iseven)
        map_U!(U[μ],mapfunc!,V,iseven)
    end

end

function SU2update_KP!(Unew,V,beta,NC,temps,ITERATION_MAX = 10^5)
    V0 = temps[1]
    temp = temps[2]

    eps = 0.000000000001

    ρ0 = real(V[1,1]+V[2,2])/2
    ρ1 = -imag(V[1,2]+V[2,1])/2
    #ρ1 = imag(V[1,2]+V[2,1])/2
    ρ2 = real(V[2,1]-V[1,2])/2
    ρ3 = imag(V[2,2]-V[1,1])/2
    ρ = sqrt(ρ0^2+ρ1^2+ρ2^2+ρ3^2)
    #println("R = ",R," ρ ",ρ)
    #println("detV = , ", det(V)," ",ρ0^2+ρ1^2+ρ2^2+ρ3^2)
    detV = V[1,1]*V[2,2] - V[1,2]*V[2,1]
    V0[1,1] = ρ*V[2,2]/detV
    V0[1,2] = -ρ*V[1,2]/detV
    V0[2,1] = -ρ*V[2,1]/detV
    V0[2,2] = ρ*V[1,1]/detV
    #println(V0)
    #V0 = inv(V/ρ)

    #
    #Nc = 2 # Since Ishikawa's book uses 1/g^2 notation.
    #k = (beta/NC)*ρ
    k = 2*(beta/NC)*ρ
    #println("k $k, $ρ")

    #k = (beta/2)*ρ

    R = rand() + eps
    Rp = rand() + eps
    X = -log(R)/k
    Xp = -log(Rp)/k
    Rpp = rand()
    C = cos(2pi*Rpp)^2
    A = X*C
    delta = Xp + A
    Rppp = rand()

    a = zeros(Float64,4)
    while(Rppp^2 > 1-0.5*delta)
        R = rand()
        Rp = rand()
        X = -log(R)/k
        Xp = -log(Rp)/k
        Rpp = rand()
        C = cos(2pi*Rpp)^2
        A = X*C
        delta = Xp + A
        Rppp = rand()
        #println(Rppp^2,"\t",1-0.5*delta)
    end
    a[1] = 1-delta

    rr = sqrt(1.0-a[1]^2)
    ϕ = rand()*pi*2.0 # ϕ = [0,2pi]
    cosθ = (rand()-0.5)*2.0 # -1<cosθ<1
    sinθ = sqrt(1-cosθ^2)

    a[2]=rr*cos(ϕ)*sinθ
    a[3]=rr*sin(ϕ)*sinθ
    a[4]=rr*cosθ
    temp[1,1] = a[1]+im*a[4]
    temp[1,2] = a[3]+im*a[2]
    temp[2,1] = -a[3]+im*a[2]
    temp[2,2] = a[1]-im*a[4]
    mul!(Unew,temp,V0)
    #Unew[:,:] = Unew*V0

    α = Unew[1,1]*0.5 + conj(Unew[2,2])*0.5
    β = Unew[2,1]*0.5 - conj(Unew[1,2])*0.5

    detU = abs(α)^2 + abs(β)^2
    Unew[1,1] = α/detU
    Unew[2,1]  = β/detU
    Unew[1,2] = -conj(β)/detU
    Unew[2,2] = conj(α)/detU

end

function SU2update_KP(V,beta,NC,temps,ITERATION_MAX = 10^5)
    #println("V = ",V)
    Unew = zero(V)
    SU2update_KP!(Unew,V,beta,NC,temps,ITERATION_MAX)
    return Unew
end

function SUNupdate_matrix!(u,V,beta,NC,temps2,temps3,ITERATION_MAX)
    UV = temps3[1]
    A = temps3[2]
    AU = temps3[3]
    temp1 = temps2[1]
    temp2 = temps2[2]
    S = temps2[3]
    K = temps2[4]

    for l=1:NC
        #for l=1:2NC

        mul!(UV,u,V)
        #UV = u[:,:]*V

        n = rand(1:NC-1)#l
        m = rand(n:NC)
        while(n==m)
            m = rand(n:NC)
        end

        #=
        if l < NC
            n = l
            m = l+1
        else
            n = rand(1:NC)#l
            m = rand(1:NC)
            while(n==m)
                m = rand(1:NC)
            end
        end
        =#

        make_submatrix!(S,UV,n,m)
        #S = make_submatrix(UV,n,m)

        #gramschmidt_special!(S)
        project_onto_SU2!(S)

        SU2update_KP!(K,S,beta,NC,[temp1,temp2],ITERATION_MAX)

        make_largematrix!(A,K,n,m,NC)

        mul!(AU,A,u)

        #AU = A*u[:,:]

        u[:,:] .= AU
        #println("det U ",det(AU))

    end

    AU[:,:] .= u #u[:,:]
    normalizeN!(AU)
    u[:,:] .= AU
end

function SU3update_matrix!(u,V,beta,NC,temps2,temps3,ITERATION_MAX)
    UV = temps3[1]
    A = temps3[2]
    AU = temps3[3]
    temp1 = temps2[1]
    temp2 = temps2[2]
    S = temps2[3]
    K = temps2[4]

    #println("#Heatbath for one SU(3) link started")
    for l=1:3
        mul!(UV,u,V)
        #UV = u*V
        #println("UV $UV $V $u")

        if l==1
            n,m = 1,2
        elseif l==2
            n,m = 2,3
        else
            n,m = 1,3

        end

        #S = make_submatrix(UV,n,m)
        make_submatrix!(S,UV,n,m)
        #gramschmidt_special!(S)
        project_onto_SU2!(S)

        #K = SU2update_KP(S,beta,NC,ITERATION_MAX)
        SU2update_KP!(K,S,beta,NC,[temp1,temp2],ITERATION_MAX)

        #A = make_largematrix(K,n,m,NC)
        make_largematrix!(A,K,n,m,NC)

        mul!(AU,A,u)
        #AU = A*u

        u[:,:] .= AU
    end

    AU[:,:] .= u #u[mu][:,:,ix,iy,iz,it]
    normalize3!(AU)
    u[:,:] .= AU
    #u[mu][:,:,ix,iy,iz,it] = AU
end

function SUN_overrelaxation!(u,V,beta,NC,temps2,temps3,ITERATION_MAX)
    UV = temps3[1]
    A = temps3[2]
    AU = temps3[3]
    temp1 = temps2[1]
    temp2 = temps2[2]
    w = temps2[3]
    h = temps2[4]
    #K = temps2[5]

    for l=1:NC
        #for l=1:2NC

        mul!(UV,u,V)
        #UV = u[:,:]*V

        n = rand(1:NC-1)#l
        m = rand(n:NC)
        while(n==m)
            m = rand(n:NC)
        end

        # we emplay DeGrand's textbook notation
        make_submatrix!(w,UV,n,m)
        #S = make_submatrix(UV,n,m)

        #gramschmidt_special!(S)
        project_onto_SU2!(w)

        # following two lines are only difference to HB
        #SU2update_KP!(K,S,beta,NC,[temp1,temp2],ITERATION_MAX)
        for j=1:2
            for i=1:2
                h[i,j] = 0
                for k=1:2
                    h[i,j] += w'[i,k]*w'[k,j]
                end
            end
        end
        normalizeN!(h)

        make_largematrix!(A,h,n,m,NC)

        mul!(AU,A,u)

        #AU = A*u[:,:]

        u[:,:] .= AU
        #println("det U ",det(AU))

    end

    AU[:,:] .= u #u[:,:]
    normalizeN!(AU)
    u[:,:] .= AU
end

function project_onto_SU2!(S) # This project onto SU(2) up to normalization.
    #S2 = zeros(ComplexF64,2,2)
    α = S[1,1]*0.5 + conj(S[2,2])*0.5
    β = S[2,1]*0.5 - conj(S[1,2])*0.5
    S[1,1] = α
    S[2,1] = β
    S[1,2] = -conj(β)
    S[2,2] = conj(α)
    #return S2
end

function make_submatrix(UV,i,j)
    S = zeros(ComplexF64,2,2)
    S[1,1] = UV[i,i]
    S[1,2] = UV[i,j]
    S[2,1] = UV[j,i]
    S[2,2] = UV[j,j]
    return S
end

function make_submatrix!(S,UV,i,j)
    S[1,1] = UV[i,i]
    S[1,2] = UV[i,j]
    S[2,1] = UV[j,i]
    S[2,2] = UV[j,j]
    return
end

function make_largematrix(K,i,j,NC)
    A = zeros(ComplexF64,NC,NC)
    for n=1:NC
        A[n,n] = 1
    end
    #K = project_onto_su2(K)
    A[i,i] = K[1,1]
    A[i,j] = K[1,2]
    A[j,i] = K[2,1]
    A[j,j] = K[2,2]
    return A
end

function make_largematrix!(A,K,i,j,NC)
    #A = zeros(ComplexF64,NC,NC)
    A .= 0
    for n=1:NC
        A[n,n] = 1
    end
    #K = project_onto_su2(K)
    A[i,i] = K[1,1]
    A[i,j] = K[1,2]
    A[j,i] = K[2,1]
    A[j,j] = K[2,2]
    return A
end

const nhit = 6
const rwidth = 0.4

"""
-------------------------------------------------c
 su2-submatrix(c) in su3 matrix(x)
        su2            su3
 k=1         <-    1-2 elements
 k=2         <-    2-3 elements
 k=3         <-    1-3 elements
 k=4          ->   1-2 elements
 k=5          ->   2-3 elements
 k=6          ->   1-3 elements
-------------------------------------------------c
"""
function submat!(x,c,n,k,id)

    if k==1
        for i=1:n
            c[1,i] = real(x[1,1,i]+x[2,2,i])*0.5
            c[2,i] = imag(x[1,2,i]+x[2,1,i])*0.5
            c[3,i] = real(x[1,2,i]-x[2,1,i])*0.5
            c[4,i] = imag(x[1,1,i]-x[2,2,i])*0.5
        end
    elseif k==2
        for i=1:n
            c[1,i] = real(x[2,2,i]+x[3,3,i])*0.5
            c[2,i] = imag(x[3,2,i]+x[2,3,i])*0.5
            c[3,i] = real(x[3,2,i]-x[2,3,i])*0.5
            c[4,i] = imag(x[2,2,i]-x[3,3,i])*0.5
        end

    elseif k==3
        for i=1:n
            c[1,i] = real(x[1,1,i]+x[3,3,i])*0.5
            c[2,i] = imag(x[3,1,i]+x[1,3,i])*0.5
            c[3,i] = real(x[1,3,i]-x[3,1,i])*0.5
            c[4,i] = imag(x[1,1,i]-x[3,3,i])*0.5
        end
    elseif k==4

        for i=1:n
            #println("i = $i")
            #println(c[:,i])
            if id[i] == 1
                x[1,1,i] = c[1,i] + im*c[4,i]
                x[1,2,i] = c[3,i] + im*c[2,i]
                x[1,3,i] = 0
                x[2,1,i] = -c[3,i] + im*c[2,i]
                x[2,2,i] = c[1,i] - im*c[4,i]
                x[2,3,i] = 0
                x[3,1,i] = 0
                x[3,2,i] = 0
                x[3,3,i] = 1

            elseif id[i] == 0
                x[1,1,i] = 1
                x[1,2,i] = 0
                x[1,3,i] = 0
                x[2,1,i] = 0
                x[2,2,i] = 1
                x[2,3,i] = 0
                x[3,1,i] = 0
                x[3,2,i] = 0
                x[3,3,i] = 1
            end
        end
    elseif k==5
        for i=1:n
            if id[i] == 1
                x[1,1,i] = 1
                x[1,2,i] = 0
                x[1,3,i] = 0
                x[2,1,i] = 0
                x[2,2,i] = c[1,i] + im*c[4,i]
                x[2,3,i] = -c[3,i] + im*c[2,i]
                x[3,1,i] = 0
                x[3,2,i] = c[3,i] + im*c[2,i]
                x[3,3,i] = c[1,i] -im*c[4,i]

            elseif id[i] == 0
                x[1,1,i] = 1
                x[1,2,i] = 0
                x[1,3,i] = 0
                x[2,1,i] = 0
                x[2,2,i] = 1
                x[2,3,i] = 0
                x[3,1,i] = 0
                x[3,2,i] = 0
                x[3,3,i] = 1
            end
        end

    elseif k==6
        for i=1:n
            if id[i] == 1
                x[1,1,i] = c[1,i] + im*c[4,i]
                x[1,2,i] = 0
                x[1,3,i] = c[3,i] + im*c[2,i]
                x[2,1,i] = 0
                x[2,2,i] = 1
                x[2,3,i] = 0
                x[3,1,i] = -c[3,i] + im*c[2,i]
                x[3,2,i] = 0
                x[3,3,i] = c[1,i] -im*c[4,i]

            elseif id[i] == 0
                x[1,1,i] = 1
                x[1,2,i] = 0
                x[1,3,i] = 0
                x[2,1,i] = 0
                x[2,2,i] = 1
                x[2,3,i] = 0
                x[3,1,i] = 0
                x[3,2,i] = 0
                x[3,3,i] = 1
            end
        end
    end
end

function rndprd!(ranf,n)
    rn = zeros(Float64,n)
    rndprd!(ranf,rn,n)
    return rn
end

function rndprd!(ranf,rn,n)
    for i=1:n
        rn[i] = ranf()
    end
    return rn
end

function rndprd2!(ranf,n)
    xrn = zeros(Float64,3,n)
    rndprd2!(ranf,xrn,n)
    return xrn
end

function rndprd2!(ranf,xrn,n)
    for j=1:n
        for i=1:3
            xrn[i,j] = ranf()
        end
    end
    return
end

end


text:
using JuMP, GLPK, JuDGE, Gurobi, Random, Test

include("solvers/setup_gurobi.jl")
#include("solvers/setup_cplex.jl")
#include("solvers/setup_coin.jl")
#include("solvers/setup_glpk.jl")

function inventory()
	mytree = narytree(4,4)

	Random.seed!(1000)
	price_array=rand(length(collect(mytree)))
	price=Dict(zip(collect(mytree),price_array))

	function sub_problems(n)
	    sp = JuMP.Model(JuDGE_SP_Solver)

	    @state(sp, -50<=Δstock<=50, state_name=stock, lb=0, ub=200, initial=0)

	    @ongoingcosts(sp, 0.01*stock)

	    @variable(sp, buy>=0)
		@variable(sp, sell>=0)
		#@variable(sp, change, Int)
		#@constraint(sp, Δstock==10*change)
	    @constraint(sp, trade, Δstock == buy - sell)

	    @objective(sp, Min, (buy*1.01)*price[n]-(sell/1.01)*price[n])

	    sp
	end

	model = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver,check=true, risk=JuDGE.Risk(0.1,bound=0.0))

	JuDGE.solve(model)
	JuDGE.print_expansions(model,onlynonzero=true,inttol=10^-5)
	println("\nRe-solved Objective: " * string(resolve_subproblems(model)))
	solution=JuDGE.solution_to_dictionary(model)
	solution[:prices]=price
	JuDGE.visualize_tree(mytree,solution)
	JuDGE.get_objval(model)
end

@test inventory() ≈  -17.708 atol = 5e-3